ロトカ＝ヴォルテラ方程式

非線形微分方程式を考えましょう。

以下は、ロトカ＝ヴォルテラ方程式と呼ばれる非線形微分方程式の例です。

これをPythonで解き、その解をグラフ化します。

$$
[ \frac{dx}{dt} = \alpha x - \beta xy ]
$$

$$
[ \frac{dy}{dt} = \delta xy - \gamma y ]
$$

ここで、$ (x) $と$ (y) $は時間$ (t) $の関数です。

この方程式は、捕食者と被食者の関係をモデル化したもので、$ (x) $が被食者（例えばウサギ）、$ (y) $が捕食者（例えばキツネ）の個体数を表します。

それでは、この方程式を解いて結果をグラフ化するPythonコードを見てみましょう。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# ロトカ＝ヴォルテラ方程式
def predator_prey_equations(xy, t, alpha, beta, delta, gamma):
    x, y = xy
    dxdt = alpha * x - beta * x * y
    dydt = delta * x * y - gamma * y
    return [dxdt, dydt]

# 初期条件とパラメータ
x0 = 10  # 初期被食者の個体数
y0 = 5   # 初期捕食者の個体数
alpha = 1.1
beta = 0.4
delta = 0.1
gamma = 0.4

# 時間
t = np.linspace(0, 30, 1000)

# 微分方程式を解く
solution = odeint(predator_prey_equations, [x0, y0], t, args=(alpha, beta, delta, gamma))
x, y = solution[:, 0], solution[:, 1]

# 結果をグラフ化
plt.figure(figsize=(8, 6))
plt.plot(t, x, label='被食者 (Prey)')
plt.plot(t, y, label='捕食者 (Predator)')
plt.xlabel('時間')
plt.ylabel('個体数')
plt.title('ロトカ＝ヴォルテラ方程式の解')
plt.legend()
plt.grid(True)
plt.show()

このコードでは、Scipyのodeint関数を使用して微分方程式を数値的に解いています。

結果を時系列でグラフ化し、被食者と捕食者の個体数の変化を示しています。

ソースコード解説

このソースコードは、Pythonでロトカ＝ヴォルテラ方程式（捕食者と被食者の間の相互作用をモデル化する非線形微分方程式）を解き、結果をグラフ化するものです。

それぞれの部分について順番に説明します。

ライブラリのインポート

1
2
3

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

この部分では、NumPy（数値計算用のライブラリ）、Scipyのodeint（常微分方程式を数値的に解くための関数）、Matplotlib（グラフ描画ライブラリ）をインポートしています。

ロトカ＝ヴォルテラ方程式の定義

1
2
3
4
5

def predator_prey_equations(xy, t, alpha, beta, delta, gamma):
    x, y = xy
    dxdt = alpha * x - beta * x * y
    dydt = delta * x * y - gamma * y
    return [dxdt, dydt]

この部分では、ロトカ＝ヴォルテラ方程式を関数predator_prey_equationsとして定義しています。

この方程式は被食者と捕食者の個体数の時間変化を示します。

初期条件とパラメータの設定

1
2
3
4
5
6

x0 = 10  # 初期被食者の個体数
y0 = 5   # 初期捕食者の個体数
alpha = 1.1
beta = 0.4
delta = 0.1
gamma = 0.4

初期条件（被食者と捕食者の初期個体数）と方程式のパラメータ（alpha, beta, delta, gamma）を設定しています。

時間の設定

1

t = np.linspace(0, 30, 1000)

$0$から$30$までの時間を$1000分割$した配列を生成しています。

微分方程式の解を求める

1
2

solution = odeint(predator_prey_equations, [x0, y0], t, args=(alpha, beta, delta, gamma))
x, y = solution[:, 0], solution[:, 1]

odeint関数を使用して微分方程式を数値的に解き、被食者と捕食者の個体数の時間変化を求めています。

結果のグラフ化

1
2
3
4
5
6
7
8
9

plt.figure(figsize=(8, 6))
plt.plot(t, x, label='被食者 (Prey)')
plt.plot(t, y, label='捕食者 (Predator)')
plt.xlabel('時間')
plt.ylabel('個体数')
plt.title('ロトカ＝ヴォルテラ方程式の解')
plt.legend()
plt.grid(True)
plt.show()

Matplotlibを使って、被食者と捕食者の個体数の時間変化をグラフ化しています。

横軸は時間、縦軸は個体数を表し、それぞれのラインは被食者と捕食者を表しています。

グラフは凡例付きで、グリッド線も表示されています。

グラフ解説

このグラフは、ロトカ＝ヴォルテラ方程式に基づいてモデル化された捕食者と被食者の個体数の時間変化を表しています。
この方程式は捕食者と被食者の相互作用を記述し、個体数の変動を表す非線形微分方程式です。

横軸は時間を示し、縦軸は被食者（Prey）と捕食者（Predator）の個体数を表します。
グラフでは、被食者（ウサギなど）の個体数が青色で、捕食者（キツネなど）の個体数がオレンジ色で表示されています。

初期条件では、被食者と捕食者の個体数がそれぞれ$ x0 $と$ y0 $で設定されています。
それぞれの個体数は時間の経過とともに変化し、捕食者は被食者に依存して増減します。
被食者が増えると捕食者も増えるが、その後被食者の減少により捕食者も減少します。
この相互作用が繰り返され、周期的な変動が生じます。

グラフの振動や変動は、捕食者と被食者の間の相互作用とその影響を示しています。
このようなモデルは生態学や生物学などで捕食者と被食者の個体数の相互作用を理解するために利用されます。

Watts-Strogatzモデル

Watts-Strogatzモデルを使ってランダムなグラフを作成し、その次数分布をグラフ化してみましょう。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

import networkx as nx
import matplotlib.pyplot as plt

# Watts-Strogatzモデルによるグラフの生成
n = 100  # ノード数
k = 4    # 各ノードの近傍ノード数
p = 0.1  # エッジの再接続確率

G = nx.watts_strogatz_graph(n, k, p)

# 次数分布の取得
degrees = [G.degree(node) for node in G.nodes()]

# 次数分布のグラフ化
plt.hist(degrees, bins='auto', alpha=0.7, color='skyblue', edgecolor='black')
plt.xlabel('Degree')
plt.ylabel('Frequency')
plt.title('Degree Distribution of a Watts-Strogatz Graph')
plt.show()

このコードでは、Watts-Strogatzモデルを使用して100個のノードを持つランダムなグラフを生成しています。

その後、各ノードの次数（接続数）を取得し、次数分布をヒストグラムで描画しています。

次数分布はネットワークの特性を示す重要な指標の1つです。

これにより、ノードがどれだけ他のノードと接続しているかの様子を可視化することができます。

ソースコード解説

このコードは、PythonのNetworkXライブラリとMatplotlibを使用して、Watts-Strogatzモデルに基づくランダムなグラフを生成し、その次数分布を可視化しています。

1. ライブラリのインポート:

1
2

import networkx as nx
import matplotlib.pyplot as plt

• networkxはグラフ理論や複雑ネットワークの作成や操作を行うためのライブラリです。
• matplotlib.pyplotはグラフの描画に使用されます。

2. Watts-Strogatzモデルによるグラフの生成:

1
2
3
4
5

n = 100  # ノード数
k = 4    # 各ノードの近傍ノード数
p = 0.1  # エッジの再接続確率

G = nx.watts_strogatz_graph(n, k, p)

• nx.watts_strogatz_graph(n, k, p)はWatts-Strogatzモデルに基づいてランダムなグラフを生成します。
• nはノード数、kは各ノードの近傍ノード数、pはエッジの再接続確率を表します。

3. 次数分布の取得:

1

degrees = [G.degree(node) for node in G.nodes()]

• G.degree(node)はネットワーク内の各ノードの次数（接続数）を取得します。
• G.nodes()はグラフ内のすべてのノードを返します。
  これを利用して各ノードの次数をリストに収集します。

4. 次数分布のグラフ化:

1
2
3
4
5

plt.hist(degrees, bins='auto', alpha=0.7, color='skyblue', edgecolor='black')
plt.xlabel('Degree')
plt.ylabel('Frequency')
plt.title('Degree Distribution of a Watts-Strogatz Graph')
plt.show()

• plt.hist()は与えられたデータ（degrees）のヒストグラムを作成します。
• bins='auto'は適切なビンの数を自動的に決定します。
• alpha、color、edgecolorはヒストグラムのスタイルを設定します。
• plt.xlabel()、plt.ylabel()、plt.title()はそれぞれx軸のラベル、y軸のラベル、タイトルを設定します。
• plt.show()はグラフを表示します。

これにより、Watts-Strogatzモデルに基づくランダムグラフの次数分布が可視化されます。

次数分布は、ネットワーク内の各ノードが持つ接続の数を表し、ネットワークの構造を理解するのに役立ちます。

最適化問題

PuLPを使用して最適化問題を解いてみましょう。

例として、複数の製品を生産する工場があり、利益を最大化するためにそれぞれの製品の生産量を決定したいとします。

以下は、3つの製品 (A, B, C) の生産にかかる時間と利益、および生産できる最大数量の制約条件を考慮した最適化問題です。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

from pulp import LpMaximize, LpProblem, LpVariable

# 問題の定義
model = LpProblem(name="product_production", sense=LpMaximize)

# 変数の定義
prod_A = LpVariable(name="product_A", lowBound=0, cat="Integer")
prod_B = LpVariable(name="product_B", lowBound=0, cat="Integer")
prod_C = LpVariable(name="product_C", lowBound=0, cat="Integer")

# 利益と制約条件の設定
model += 300 * prod_A + 200 * prod_B + 500 * prod_C, "total_profit"

model += 2 * prod_A + 3 * prod_B + 4 * prod_C <= 100  # 生産時間の制約
model += prod_A + prod_B <= 40  # 製品Aと製品Bの同時生産の制約

# 問題の解決
model.solve()

# 結果の出力
print(f"Product A: {prod_A.varValue}")
print(f"Product B: {prod_B.varValue}")
print(f"Product C: {prod_C.varValue}")
print(f"Total Profit: {model.objective.value()}")

# グラフ化
import matplotlib.pyplot as plt

products = ['Product A', 'Product B', 'Product C']
production_quantities = [prod_A.varValue, prod_B.varValue, prod_C.varValue]

plt.figure(figsize=(8, 5))
plt.bar(products, production_quantities, color='skyblue')
plt.xlabel('Products')
plt.ylabel('Production Quantity')
plt.title('Optimal Production Quantities')
plt.show()

このコードを実行すると、最適な製品の生産量とその利益が表示され、それらを棒グラフで可視化することができます。

このように、PuLPを使用して最適化問題を解き、結果を分かりやすくグラフ化することができます。

[実行結果]

ソースコード解説

このコードは、PuLPとMatplotlibを使用して製品生産最適化問題を解き、結果をグラフ化するものです。

1. PuLPの利用

• from pulp import LpMaximize, LpProblem, LpVariable:
  PuLPの必要なクラスや関数をインポートしています。
• model = LpProblem(name="product_production", sense=LpMaximize):
  問題を定義するための空のモデルを作成しています。
  ここでは、最大化問題（LpMaximize）としてproduct_productionという名前の問題を作成しています。

2. 変数の定義

• prod_A, prod_B, prod_C:
  製品A、B、Cのそれぞれの生産量を表す変数を定義しています。
  LpVariableを使用して整数として設定し、0以上の値を取るように制約を設けています。

3. 利益と制約条件の設定

• model += 300 * prod_A + 200 * prod_B + 500 * prod_C, "total_profit":

目的関数を定義しています。
製品ごとの利益を考慮して、それらの生産量の総利益を最大化するように設定しています。

• model += 2 * prod_A + 3 * prod_B + 4 * prod_C <= 100:
  生産時間の制約条件を設定しています。
  各製品の生産にかかる時間が合計で100以下であることを意味します。
• model += prod_A + prod_B <= 40:
  製品AとBの同時生産量に関する制約条件を設定しています。
  製品AとBの合計生産量が40以下であることを要求しています。

4. 問題の解決

• model.solve():
  定義された最適化問題を解きます。

5. 結果の出力とグラフ化

• 最適解から得られた製品ごとの生産量と総利益を出力しています。
• matplotlibを使用して、製品ごとの最適生産量を棒グラフとして可視化しています。

これにより、PuLPを使って製品生産の最適化問題を定義し、解決し、その結果をグラフ化するプロセスが示されています。

結果解説

この問題の解は以下のようになりました：

[実行結果]

• 製品Aを40個生産し、製品Bは0個、製品Cは5個生産されます。
• これにより、総利益は14500になります。

グラフでは、製品ごとの生産量が棒グラフとして示されています。

製品Aが最大生産量である40で、製品Bは制約条件により生産されず、製品Cは5個生産されています。

この生産量のもとで最大の利益である14500が得られます。

製品Bが生産されない理由は、製品Aと製品Bの同時生産に関する制約条件を考慮したためです。

製品Cは比較的高い利益が得られるため、最適解では少量生産されています。

食糧供給最適化

Google OR-Toolsは最適化問題を解くためのライブラリであり、主にルーティング、ネットワークフロー、線形プログラミングなどの問題を扱います。

具体的な貧困問題として、一定の予算で最大数の家庭に食糧を供給するという問題を考えてみます。

この問題は整数線形プログラミング問題として定式化で、Google OR-ToolsのLinearSolverを用いて解くことができます。

それぞれの家庭が必要とする食糧の量と、運搬コスト（距離や時間に基づいた値）が既知であるとします。

以下に、この問題を解くためのPythonコードを示します。

（ただし、このコードはあくまで一例であり、実際の問題はより多くの要素や複雑な制約を考慮する必要があります。）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

from ortools.linear_solver import pywraplp
import matplotlib.pyplot as plt

# データ
families = [0, 1, 2, 3, 4]  # 表現の一部として家族を表すインデックス
food_required = [3, 7, 2, 5, 4]  # 各家族が必要とする食糧の量
transport_costs = [10, 20, 5, 15, 10]  # 各家族への食糧配達のコスト

# ベクトル全体で最小化するためのスカラー乗数
scalar_multiplier = 10

# 問題を作成
solver = pywraplp.Solver.CreateSolver('SCIP')

# 変数を作成
x = []
for i in families:
    x.append(solver.IntVar(0, 1, 'x[%i]' % i))  # x[i] = 家族iに食糧を配達するかどうか (0 または 1)

# 制約を定義
budget = 30  # 予算
solver.Add(sum(transport_costs[i] * x[i] for i in families) <= budget)

# 目的関数を定義
solver.Maximize(scalar_multiplier * sum(food_required[i] * x[i] for i in families) - sum(transport_costs[i] * x[i] for i in families))

# 問題を解く
solver.Solve()

# 結果を出力
print("Total food delivered = ", sum(food_required[i] * x[i].solution_value() for i in families))
delivered = [x[i].solution_value() for i in families]
plt.figure(figsize=(10,5))
plt.bar(range(len(families)), delivered)
plt.title("Food delivered to each family")
plt.xlabel("Family")
plt.ylabel("Food quantity")
plt.xticks(range(len(families)), families)
plt.show()

[実行結果]

ソースコード解説

このコードは、Google OR-Toolsを使用して最適な食糧配達問題を解決し、結果を可視化するものです。

詳細について見ていきましょう。

1. ライブラリのインポートとデータの定義

1
2
3
4
5
6
7

from ortools.linear_solver import pywraplp
import matplotlib.pyplot as plt

families = [0, 1, 2, 3, 4]  # 家族を表すインデックス
food_required = [3, 7, 2, 5, 4]  # 各家族の食糧必要量
transport_costs = [10, 20, 5, 15, 10]  # 各家族への食糧配達コスト
scalar_multiplier = 10  # 最小化するためのスカラー乗数

まず、Google OR-ToolsとMatplotlibをインポートし、最適化のためのデータを定義しています。
familiesは家族を表すインデックスのリスト、food_requiredは各家族の食糧必要量のリスト、transport_costsは各家族への食糧配達コストのリストです。
また、scalar_multiplierは後で使われるスカラー乗数です。

2. 問題の作成と変数の定義

1
2
3
4
5

solver = pywraplp.Solver.CreateSolver('SCIP')

x = []
for i in families:
   x.append(solver.IntVar(0, 1, 'x[%i]' % i))  # 家族iに食糧を配達するかどうか (0 または 1)

Google OR-Toolsの線形ソルバーを作成し、バイナリ変数 x を定義しています。
x[i] は家族 i に食糧を配達するかどうかを示す変数で、$0$か$1$の値を取ります。

3. 制約の定義

1
2

budget = 30
solver.Add(sum(transport_costs[i] * x[i] for i in families) <= budget)

予算制約を設定しています。
すべての家族への食糧配達コストの合計が budget 以下である必要があります。

4. 目的関数の定義と最適化

1
2

solver.Maximize(scalar_multiplier * sum(food_required[i] * x[i] for i in families) - sum(transport_costs[i] * x[i] for i in families))
solver.Solve()

目的関数は、家族への食糧配達量を最大化し、同時にコストを最小化するように設定されています。
そして、ソルバーを使って問題を解決しています。

5. 結果の出力と可視化

1
2
3
4
5
6
7
8
9

print("Total food delivered = ", sum(food_required[i] * x[i].solution_value() for i in families))
delivered = [x[i].solution_value() for i in families]
plt.figure(figsize=(10,5))
plt.bar(range(len(families)), delivered)
plt.title("Food delivered to each family")
plt.xlabel("Family")
plt.ylabel("Food quantity")
plt.xticks(range(len(families)), families)
plt.show()

最後に、配達された食糧量を表示し、Matplotlibを使用して各家族への食糧配達量を可視化しています。
print文で配達された総食糧量を表示し、棒グラフで各家族の配達量を示しています。

[実行結果]

緯度経度間の距離

地理関連の問題を解決する手法として、緯度経度間の距離を計算し、その結果をグラフ化することができます。

以下はその例です。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

from scipy.spatial.distance import cdist
import matplotlib.pyplot as plt

# 緯度経度のデータを定義
# 例: Tokyo, New York, London, Sydney
cities = {
    'Tokyo': (35.6895, 139.6917),
    'New York': (40.7128, -74.0060),
    'London': (51.5074, -0.1278),
    'Sydney': (-33.8688, 151.2093)
}

# 緯度経度をnumpy配列に変換
locations = list(cities.values())

# 距離を計算（ここではユークリッド距離を使用）
distances = cdist(locations, locations)

# 結果のグラフ化
plt.imshow(distances, cmap='viridis', interpolation='nearest')
plt.colorbar(label='Distance (degrees)')
plt.xticks(ticks=range(len(cities)), labels=cities.keys())
plt.yticks(ticks=range(len(cities)), labels=cities.keys())
plt.title('Distance between Cities')
plt.show()

このコードは、4つの都市（東京、ニューヨーク、ロンドン、シドニー）の緯度経度を使用し、それらの間の距離を計算し、結果をヒートマップとして可視化しています。

各都市間の距離が色で表され、より遠い場所ほど色が濃くなります。

こうした手法を用いて、地理的な位置関係を視覚的に表現することができます。

ソースコード解説

このコードは、PythonのSciPyライブラリを使用して、4つの都市の緯度経度からそれぞれの都市間の距離を計算し、その結果をヒートマップで可視化するものです。

詳細に見ていきましょう。

ライブラリのインポート

1
2

from scipy.spatial.distance import cdist
import matplotlib.pyplot as plt

必要なライブラリをインポートしています。

• cdist: SciPyの距離関数を使用して、配列内の要素間の距離を計算するために使用されます。
• matplotlib.pyplot: グラフ描画のためのライブラリです。

都市の緯度経度の定義

1
2
3
4
5
6

cities = {
    'Tokyo': (35.6895, 139.6917),
    'New York': (40.7128, -74.0060),
    'London': (51.5074, -0.1278),
    'Sydney': (-33.8688, 151.2093)
}

4つの都市（東京、ニューヨーク、ロンドン、シドニー）の緯度経度データを辞書形式で定義しています。

緯度経度のnumpy配列への変換

1

locations = list(cities.values())

cities辞書から値（緯度経度）を取り出し、それをnumpy配列に変換しています。

都市間の距離の計算

1

distances = cdist(locations, locations)

cdist関数を使用して、各都市間の距離を計算しています。
ここではデフォルトでユークリッド距離が使用されています。

結果のグラフ化

1
2
3
4
5
6

plt.imshow(distances, cmap='viridis', interpolation='nearest')
plt.colorbar(label='Distance (degrees)')
plt.xticks(ticks=range(len(cities)), labels=cities.keys())
plt.yticks(ticks=range(len(cities)), labels=cities.keys())
plt.title('Distance between Cities')
plt.show()

imshow関数を使って都市間の距離をヒートマップとして表示しています。
colorbar関数でカラーバーを表示し、xticksおよびyticksでx軸とy軸の目盛りを設定し、都市名をラベルとして表示しています。
最後にtitle関数でグラフのタイトルを設定し、show関数でグラフを表示しています。

Rosenbrock関数

非線形最適化問題の一例として、Rosenbrock関数の最小化問題を考えてみましょう。

Rosenbrock関数は、次のように定義されます：

$
f(x, y) = (a - x)^2 + b * (y - x^2)^2
$

ここで、$ a = 1, b = 100 $とします。

この関数の最小値は$ 0 $で、その位置は$ (a, a^2) = (1, 1) $です。

Pythonでは、Scipyの最適化関数を使用してこの最適化問題を解くことができます。

また、Matplotlibを使用して結果をグラフ化することも可能です。

以下に、PythonでRosenbrock関数の最小化問題を解くコードを示します。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt

# Rosenbrock関数を定義
def rosenbrock(x):
   return (1 - x[0])**2 + 100 * (x[1] - x[0]**2)**2

# Rosenbrock関数の勾配を定義
def rosenbrock_der(x):
   return np.array([-2 * (1 - x[0]) - 400 * x[0] * (x[1] - x[0]**2), 200 * (x[1] - x[0]**2)])

# Rosenbrock関数のHessianを定義
def rosenbrock_hess(x):
   return np.array([[2 + 1200 * x[0]**2 - 400 * x[1], -400 * x[0]],
                   [-400 * x[0], 200]])

# 初期値を設定
x0 = np.array([0.5, 0])

# 最適化問題を解く
res = minimize(rosenbrock, x0, method='trust-constr', jac=rosenbrock_der, hess=rosenbrock_hess)

# 結果を表示
print(res.x)

# 結果をグラフ化
x = np.linspace(-1, 2, 100)
y = np.linspace(-1, 3, 100)
X, Y = np.meshgrid(x, y)
Z = rosenbrock([X, Y])

plt.contour(X, Y, Z, levels=np.logspace(-2, 2, 20))
plt.xlabel('x')
plt.ylabel('y')
plt.title('Rosenbrock function')
plt.show()

このコードは、Rosenbrock関数の最小化問題を解くためのPythonのコードです。

まず、Rosenbrock関数とその勾配 Hessianを定義します。

次に、初期値を設定し、scipy.optimize.minimize関数を使用して最適化問題を解きます。

最後に、結果を表示し、結果をグラフ化します。

ソースコード解説

ソースコードを各部分を説明します。

1. ライブラリのインポート

1
2
3

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt

必要なライブラリをインポートしています。

• numpy: 数値計算用のライブラリ
• scipy.optimize.minimize: 最適化問題を解くための関数
• matplotlib.pyplot: グラフ描画のためのライブラリ

2. Rosenbrock関数の定義

1
2

def rosenbrock(x):
   return (1 - x[0])**2 + 100 * (x[1] - x[0]**2)**2

Rosenbrock関数を定義しています。
この関数は最小化問題でよく使用される非線形関数で、$ (f(x, y) = (1 - x)^2 + 100(y - x^2)^2) $です。

3. Rosenbrock関数の勾配とHessianの定義

1
2
3
4
5
6

def rosenbrock_der(x):
   return np.array([-2 * (1 - x[0]) - 400 * x[0] * (x[1] - x[0]**2), 200 * (x[1] - x[0]**2)])

def rosenbrock_hess(x):
   return np.array([[2 + 1200 * x[0]**2 - 400 * x[1], -400 * x[0]],
                   [-400 * x[0], 200]])

Rosenbrock関数の勾配とHessian（ヘシアン、2階偏微分行列）をそれぞれ定義しています。

4. 初期値の設定と最適化問題の解法

1
2

x0 = np.array([0.5, 0])
res = minimize(rosenbrock, x0, method='trust-constr', jac=rosenbrock_der, hess=rosenbrock_hess)

初期値を設定し、minimize関数を使ってRosenbrock関数の最小化問題を解いています。
method='trust-constr'は信頼領域法を用いて最適化を行います。
jac引数とhess引数はそれぞれ勾配とHessianの関数を指定しています。

5. 結果の表示

1

print(res.x)

最適化結果（最小値に達する$ x $の値）を出力しています。

6. 結果のグラフ化

1
2
3
4
5
6
7
8
9
10

x = np.linspace(-1, 2, 100)
y = np.linspace(-1, 3, 100)
X, Y = np.meshgrid(x, y)
Z = rosenbrock([X, Y])

plt.contour(X, Y, Z, levels=np.logspace(-2, 2, 20))
plt.xlabel('x')
plt.ylabel('y')
plt.title('Rosenbrock function')
plt.show()

Rosenbrock関数の等高線プロットを作成しています。
X軸とY軸の値を用意し、等高線の高さを関数値に対応させて関数を可視化しています。

結果解説

このグラフは、Rosenbrock関数の3次元表示を示しています。
Rosenbrock関数は、非線形最適化問題の一つで、最小値は$ 0 $で、その位置は$ (a, a^2) = (1, 1) $です。

グラフの$x軸$は、$x$の値を表し、$y軸$は$y$の値を表します。
Z軸は、Rosenbrock関数の値を表します。
この関数の最小値は$ 0 $で、その位置は$ (a, a^2) = (1, 1) $です。

グラフの中心部分は、関数の最小値を示しています。
この部分は、$ x=1 $と$ y=1 $の位置で、これが関数の最小値を示しています。

グラフの上部は、関数の値が大きい部分を示しています。
この部分は、$ x $や$ y $の値が大きいと、関数の値も大きくなります。

グラフの下部は、関数の値が小さい部分を示しています。
この部分は、$ x $や$ y $の値が小さいと、関数の値も小さくなります。

このグラフを見ると、Rosenbrock関数の最小値は、$ x=1 $と$ y=1 $の位置で、これが関数の最小値を示しています。

この位置は、$ x $や$ y $の値が最も小さいときに達します。
そのため、この位置は、最適化問題を解くための最良の解となります。

月と地球の運動

月と地球の運動をPythonを使って3Dで表現してみます。

このコードでは、地球が静止していると仮定し、月の動きのみを描画します。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np

# 定数定義
G = 6.67430e-11  # 万有引力定数
m1 = 5.972e24  # 地球の質量(kg)
m2 = 7.35e22  # 月の質量(kg)
r = 384400000  # 地球から月までの距離(m)
v1 = 0  # 地球の初期速度(m/s)
v2 = 1022  # 月の初期速度(m/s)
timestep = 1000  # タイムステップ
num_steps = 10000  # ステップ数

# 初期条件
x1, y1, z1 = 0, 0, 0  # 地球の初期位置
x2, y2, z2 = r, 0, 0  # 月の初期位置

# 速度成分
vx1, vy1, vz1 = v1, 0, 0
vx2, vy2, vz2 = 0, v2, 0

# 位置、速度の更新
x1_values, y1_values, z1_values = [], [], []
x2_values, y2_values, z2_values = [], [], []

for _ in range(num_steps):
    # 地球と月の距離
    dx = x2 - x1
    dy = y2 - y1
    dz = z2 - z1
    r = np.sqrt(dx ** 2 + dy ** 2 + dz ** 2)

    # 万有引力の計算
    F = G * m1 * m2 / r ** 2

    # 重力の影響による加速度
    ax1 = F * dx / r / m1
    ay1 = F * dy / r / m1
    az1 = F * dz / r / m1
    ax2 = -F * dx / r / m2
    ay2 = -F * dy / r / m2
    az2 = -F * dz / r / m2

    # 速度の更新
    vx1 += ax1 * timestep
    vy1 += ay1 * timestep
    vz1 += az1 * timestep
    vx2 += ax2 * timestep
    vy2 += ay2 * timestep
    vz2 += az2 * timestep

    # 位置の更新
    x1 += vx1 * timestep
    y1 += vy1 * timestep
    z1 += vz1 * timestep
    x2 += vx2 * timestep
    y2 += vy2 * timestep
    z2 += vz2 * timestep

    # 位置を保存
    x1_values.append(x1)
    y1_values.append(y1)
    z1_values.append(z1)
    x2_values.append(x2)
    y2_values.append(y2)
    z2_values.append(z2)

# 3D軌道のプロット
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(x2_values, y2_values, z2_values, label='Moon')
ax.plot([0], [0], [0], marker='o', markersize=8, color='blue', label='Earth')
ax.set_xlabel('X coordinate (m)')
ax.set_ylabel('Y coordinate (m)')
ax.set_zlabel('Z coordinate (m)')
ax.set_title('Moon Orbit around Earth')
ax.legend()
plt.show()

このコードは、月の軌道を3Dで描画し、地球を静止して描画します。

月は地球の周りを楕円軌道で回る様子が視覚化されます。

ソースコード解説

コードの各部分を詳しく見ていきましょう：

1. ライブラリのインポート:

• import matplotlib.pyplot as plt: プロット用のMatplotlibライブラリをインポートします。
• from mpl_toolkits.mplot3d import Axes3D: 3Dグラフを作成するためのAxes3Dをインポートします。
• import numpy as np: 数値計算のためのNumPyライブラリをインポートします。

2. 定数の定義:

• G = 6.67430e-11など、物理定数（万有引力定数など）や天体の質量、初期位置、速度などの値を定義します。

3. 初期条件:

• x1, y1, z1など、地球と月の初期位置を設定します。

4. 速度成分:

• vx1, vy1, vz1など、地球と月の初期速度を設定します。

5. 位置、速度の更新:

• forループ内で、地球と月の位置および速度を時間ステップごとに計算して更新します。ニュートンの運動方程式と万有引力の法則に基づいて、加速度を計算し、それをもとに速度と位置を更新します。

6. 位置の保存:

• 各ステップでの地球と月の位置を保存するためのリストを作成し、更新された位置をそれぞれのリストに追加します。

7. 3D軌道のプロット:

• 地球と月の位置を元に、3Dグラフを作成します。
• figとaxを使用して3Dグラフを作成し、plotメソッドで月の軌道を、plotメソッドとmarkerを使用して地球を示します。
• 軸ラベル、タイトル、凡例を設定しています。

最後に、plt.show()で作成したグラフを表示します。

このコードは、地球と月の間の引力による運動をシミュレートし、月の楕円軌道を地球の周りにプロットしています。

結果解説

このグラフは、地球と月の間の万有引力を考慮した際の月の軌道を3Dで表現したものです。

青い点が地球を表し、月の軌道がその周りを回る様子が示されています。
地球は固定された位置にあり、月はその周りを楕円軌道で動いています。

月の軌道は、地球に引かれる重力によって決まります。
月が地球から遠ざかるときには引力が弱まり、近づくときには引力が強くなります。
その結果、月は楕円軌道を描きます。

この3Dグラフは、地球を中心に月が回転しながら地球の周りを周回する様子を視覚的に表現しています。
月の軌道が立体的に示されており、空間内での運動がわかりやすくなっています。

3D散布図／3Dサーフェスプロット

Plotlyを使って3Dグラフを描く方法はいくつかあります。

以下に、3D散布図と3Dサーフェスプロットの例を示します。

まず、必要なライブラリをインポートします。

1
2

import plotly.graph_objects as go
import numpy as np

3D散布図

以下のコードは、3D散布図を描く例です。

ここでは、ヘリックスの方程式を使って3D空間に点を配置します。

1
2
3
4
5
6

# ヘリックスの方程式
t = np.linspace(0, 10, 50)
x, y, z = np.cos(t), np.sin(t), t

fig = go.Figure(data=[go.Scatter3d(x=x, y=y, z=z, mode='markers')])
fig.show()

3Dサーフェスプロット

以下のコードは、3Dサーフェスプロットを描く例です。

ここでは、3D空間にサーフェスを配置します。

1
2
3
4
5
6

x = np.outer(np.linspace(-2, 2, 30), np.ones(30))
y = x.copy().T
z = np.cos(x ** 2 + y ** 2)

fig = go.Figure(data=[go.Surface(x=x, y=y, z=z)])
fig.show()

これらのコードは、Plotlyを使って3D散布図と3Dサーフェスプロットを描きます。

ソースコード解説

ソースコードの詳細を説明します。

まず、必要なライブラリをインポートします。

1
2

import plotly.graph_objects as go
import numpy as np

3D散布図

次に、ヘリックスの方程式を使用して3D空間に点を配置します。

ヘリックスの方程式は、3次元空間でのヘリックスの形状を定義する数学的な式です。

ここでは、時間tに対してxとyの位置を計算し、zは時間t自体とします。

1
2
3

# ヘリックスの方程式
t = np.linspace(0, 10, 50)
x, y, z = np.cos(t), np.sin(t), t

これらの位置を使用して3D散布図を描画します。

go.Scatter3dは、3D空間での散布図を作成するためのPlotlyの関数です。

x=x, y=y, z=zは、散布図の座標を指定し、mode='markers'は、各点をマーカーとして表示することを指定します。

1
2

fig = go.Figure(data=[go.Scatter3d(x=x, y=y, z=z, mode='markers')])
fig.show()

3Dサーフェスプロット

次に、3Dサーフェスプロットを描画します。

3Dサーフェスプロットは、3次元空間でのサーフェスを表示するためのプロットです。

ここでは、3D空間でのサーフェスの形状を定義する数学的な式を使用します。

1
2
3

x = np.outer(np.linspace(-2, 2, 30), np.ones(30))
y = x.copy().T
z = np.cos(x ** 2 + y ** 2)

これらの位置を使用して3Dサーフェスプロットを描画します。

go.Surfaceは、3Dサーフェスプロットを作成するためのPlotlyの関数です。

x=x, y=y, z=zは、サーフェスの座標を指定します。

1
2

fig = go.Figure(data=[go.Surface(x=x, y=y, z=z)])
fig.show()

このコードは、Plotlyを使用して3D散布図と3Dサーフェスプロットを描画します。

カーネルメソッド

カーネルメソッドは、データを高次元空間にマッピングすることで、高次元空間上での計算を低次元空間上で行うことで計算効率を向上させる手法の一つです。

以下に、ガウスカーネルを用いたカーネルメソッドの具体的な例を示します。

まず、必要なライブラリをインポートします。

1
2

import numpy as np
import matplotlib.pyplot as plt

次に、カーネル関数を定義します。
ここでは、ガウスカーネルを使用します。

1
2

def kernel(x, y):
 return np.exp(-0.5*(x - y)**2 / 2)

次に、目的関数を定義します。

1
2

def f(x):
 return 2 * x + 5 * np.cos(x)

次に、データを生成します。

1
2
3
4
5
6
7

np.random.seed(seed=4)
n = 13
m = 100
X = np.linspace(0, 10, m)
y_truth = f(X)
x = 0.5 + 8.5 * np.random.rand(n)
y = f(x)

次に、カーネルマトリックスを計算します。

1

K = np.array([[kernel(x[i], x[j]) for j in range(n)] for i in range(n)])

次に、カーネルメソッドによる予測を行います。

1
2
3
4
5
6
7
8

c = np.linalg.inv(K) @ y
y_pred = []
for i in range(m):
 sum = 0
 for j in range(n):
   sum += c[j] * kernel(x[j], X[i])
 y_pred.append(sum)
y_pred = np.array(y_pred)

最後に、予測結果をグラフに表示します。

1
2
3
4
5
6

fig, ax = plt.subplots()
ax.scatter(x, y, label="data")
ax.plot(X, y_pred, color="red", label="predict")
ax.plot(X, y_truth, label="ground truth")
ax.legend()
plt.show()

このコードは、ガウスカーネルを用いたカーネルメソッドの具体的な例を示しています。

結果解説

カーネルメソッドの例題で生成されるグラフは、データ点、予測値、真の値を表示したものです。

データ点:

青い点で表示されています。
これは、元のデータ点を表しています。

予測値:

赤い実線で表示されています。
これは、カーネルメソッドによる予測値を表しています。

真の値:

青い実線で表示されています。
これは、真の値（この場合は、カーネルメソッドが予測しようとしている関数の値）を表しています。

このグラフは、カーネルメソッドがどのようにデータを高次元空間にマッピングし、その高次元空間上での計算を低次元空間上で行っているかを視覚的に示しています。

また、カーネルメソッドがどの程度真の値を正確に予測できているかも示しています。

制約充足問題

簡単な制約充足問題を解く例として、次の条件を満たす整数の組み合わせを見つける問題を考えましょう。

条件:

1. 3つの整数$ ( x, y, z ) $を見つける。
2. $ ( x, y, z ) $は 1 から 10 までの整数。
3. $ ( x + y + z = 20 ) $の条件を満たす。

これをpython-constraintを使って解きます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

from constraint import Problem

def find_numbers():
    problem = Problem()

    # 変数の定義
    problem.addVariable('x', range(1, 11))
    problem.addVariable('y', range(1, 11))
    problem.addVariable('z', range(1, 11))

    # 制約の設定
    def constraint_function(x, y, z):
        if x + y + z == 20:
            return True

    problem.addConstraint(constraint_function, ['x', 'y', 'z'])

    solutions = problem.getSolutions()
    return solutions

solutions = find_numbers()
if solutions:
    print("整数の組み合わせ:")
    for solution in solutions:
        print(f"x: {solution['x']}, y: {solution['y']}, z: {solution['z']}")
else:
    print("条件を満たす組み合わせは見つかりませんでした。")

このコードは、1から10までの整数$ (x, y, z) $の組み合わせを見つけて、条件$ (x + y + z = 20) $を満たす解を見つけます。

解が見つかれば、整数の組み合わせが出力されます。

[実行結果]

整数の組み合わせ:
x: 10, y: 9, z: 1
x: 10, y: 8, z: 2
x: 10, y: 7, z: 3
x: 10, y: 6, z: 4
x: 10, y: 5, z: 5
x: 10, y: 4, z: 6
x: 10, y: 3, z: 7
x: 10, y: 2, z: 8
x: 10, y: 1, z: 9
x: 9, y: 10, z: 1
x: 9, y: 9, z: 2
x: 9, y: 8, z: 3
x: 9, y: 7, z: 4
x: 9, y: 6, z: 5
x: 9, y: 5, z: 6
x: 9, y: 4, z: 7
x: 9, y: 3, z: 8
x: 9, y: 2, z: 9
x: 9, y: 1, z: 10
x: 8, y: 10, z: 2
x: 8, y: 9, z: 3
x: 8, y: 8, z: 4
x: 8, y: 7, z: 5
x: 8, y: 6, z: 6
x: 8, y: 5, z: 7
x: 8, y: 4, z: 8
x: 8, y: 3, z: 9
x: 8, y: 2, z: 10
x: 7, y: 10, z: 3
x: 7, y: 9, z: 4
x: 7, y: 8, z: 5
x: 7, y: 7, z: 6
x: 7, y: 6, z: 7
x: 7, y: 5, z: 8
x: 7, y: 4, z: 9
x: 7, y: 3, z: 10
x: 6, y: 10, z: 4
x: 6, y: 9, z: 5
x: 6, y: 8, z: 6
x: 6, y: 7, z: 7
x: 6, y: 6, z: 8
x: 6, y: 5, z: 9
x: 6, y: 4, z: 10
x: 5, y: 10, z: 5
x: 5, y: 9, z: 6
x: 5, y: 8, z: 7
x: 5, y: 7, z: 8
x: 5, y: 6, z: 9
x: 5, y: 5, z: 10
x: 4, y: 10, z: 6
x: 4, y: 9, z: 7
x: 4, y: 8, z: 8
x: 4, y: 7, z: 9
x: 4, y: 6, z: 10
x: 3, y: 10, z: 7
x: 3, y: 9, z: 8
x: 3, y: 8, z: 9
x: 3, y: 7, z: 10
x: 2, y: 10, z: 8
x: 2, y: 9, z: 9
x: 2, y: 8, z: 10
x: 1, y: 10, z: 9
x: 1, y: 9, z: 10

ソースコード解説

このコードは、Pythonのpython-constraintライブラリを使用して、特定の制約を満たす整数の組み合わせを見つける例です。

1. ライブラリのインポートと関数定義:

• from constraint import Problem：constraintライブラリからProblemクラスをインポートします。
• find_numbers()関数：制約充足問題を解くための関数を定義します。

2. 制約充足問題の定義:

• problem = Problem()：制約充足問題を定義するためのProblemオブジェクトを作成します。
• problem.addVariable('x', range(1, 11))：変数 x を1から10の整数の範囲で定義します。
• problem.addVariable('y', range(1, 11))：変数 y を1から10の整数の範囲で定義します。
• problem.addVariable('z', range(1, 11))：変数 z を1から10の整数の範囲で定義します。

3. 制約の設定:

• constraint_function()関数：x + y + z = 20 となるような条件を満たす制約関数を定義します。
• problem.addConstraint(constraint_function, ['x', 'y', 'z'])：constraint_functionを使用して、x + y + z = 20 の制約を追加します。

4. 問題の解決:

• problem.getSolutions()：Problemオブジェクトを解き、条件を満たす整数の組み合わせを探します。

5. 解の表示:

• solutionsリストに解が含まれている場合は、各解の x, y, z の値を表示します。
• 解が見つからなかった場合は、条件を満たす組み合わせが見つからなかったことを通知します。