複雑な数式

次のような多項式関数を考えてみましょう。

$$
f(x) = x^4 - 6x^3 + 11x^2 - 6x
$$

この多項式関数を Python で解析し、グラフ化することができます。

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

import numpy as np
import matplotlib.pyplot as plt

# 多項式関数の定義
def f(x):
    return x**4 - 6*x**3 + 11*x**2 - 6*x

# xの値を生成
x_values = np.linspace(-1, 4, 400)  # グラフ化する範囲を指定

# 対応するyの値を計算
y_values = f(x_values)

# グラフのプロット
plt.figure(figsize=(8, 6))
plt.plot(x_values, y_values, label='f(x) = x^4 - 6x^3 + 11x^2 - 6x', color='blue')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Graph of a Complex Polynomial Function')
plt.legend()
plt.grid(True)
plt.show()

このコードは、$ ( f(x) = x^4 - 6x^3 + 11x^2 - 6x ) $の関数を生成し、範囲$ ([-1, 4]) $の$x値$に対応する$y値$を計算してグラフ化します。

[実行結果]

このようにして、複雑な数式の挙動を視覚化することができます。

ソースコード解説

このソースコードは、Pythonを使用して複雑な多項式関数のグラフを描画するものです。

以下に、コードの章立てと詳細な説明を示します。

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

1
2

import numpy as np
import matplotlib.pyplot as plt

• numpyは数値計算を行うためのライブラリで、特に配列や行列の操作が得意です。
• matplotlib.pyplotはグラフの描画をサポートするライブラリです。

2. 多項式関数の定義

1
2

def f(x):
    return x**4 - 6*x**3 + 11*x**2 - 6*x

• f(x)は4次の多項式関数を表しています。

3. xの値の生成

1

x_values = np.linspace(-1, 4, 400)

• -1から4までの範囲を等間隔に区切り、400点のxの値を生成しています。

4. 対応するyの値の計算

1

y_values = f(x_values)

• 生成したxの値に対応するyの値を関数f(x)を用いて計算しています。

5. グラフのプロット

1
2
3
4
5
6
7
8

plt.figure(figsize=(8, 6))
plt.plot(x_values, y_values, label='f(x) = x^4 - 6x^3 + 11x^2 - 6x', color='blue')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Graph of a Complex Polynomial Function')
plt.legend()
plt.grid(True)
plt.show()

• plt.figure(figsize=(8, 6))で図のサイズを指定しています。
• plt.plot(x_values, y_values, label='f(x) = x^4 - 6x^3 + 11x^2 - 6x', color='blue')でグラフをプロットしています。ラベルと色も指定しています。
• plt.xlabel('x')とplt.ylabel('f(x)')でx軸とy軸のラベルを指定しています。
• plt.title('Graph of a Complex Polynomial Function')でグラフのタイトルを指定しています。
• plt.legend()で凡例を表示します。
• plt.grid(True)でグリッドを表示します。
• plt.show()でグラフを表示します。

このコードは、指定された多項式関数をグラフに描画し、その形状や特徴を可視化しています。

結果解説

[実行結果]

上記のコードにより描画されるグラフは、4次の多項式関数$ f(x) = x^4 - 6x^3 + 11x^2 - 6x $の形状を示しています。

以下に、グラフの詳細な説明を行います。

1. 関数の形状:

グラフは典型的な4次関数の形状を示しています。
関数の次数が高いため、複雑な曲線が描かれています。

2. xとyの範囲:

$x軸$は$ -1 $から$ 4 $までの範囲で、$y軸$は関数$ f(x) $の値に対応しています。
この範囲で関数がどのように振る舞っているかが確認できます。

3. 関数のゼロ点:

グラフがx軸と交わる点が関数のゼロ点です。
これらの点は、関数の解や根を表しています。

4. 極小点と極大点:

グラフ上の極小点や極大点は、関数の極小値や極大値に対応しています。
これらの点は曲線の谷や山を示しています。

5. 凡例:

グラフには凡例が表示されており、「f(x) = x^4 - 6x^3 + 11x^2 - 6x」というラベルが付いています。
これは描画されている曲線がどの関数に対応しているかを示しています。

6. タイトル:

グラフの上部には「Graph of a Complex Polynomial Function」というタイトルがあり、描画されている関数が複雑な多項式関数であることを示しています。

7. グリッド:

グラフにはグリッドが表示されており、目盛りが曲線上の位置を正確に読み取るのに役立ちます。

このグラフは、数学的な関数の視覚的な表現として使用され、関数の特徴や挙動を理解するのに役立ちます。

DEAP（Distributed Evolutionary Algorithms in Python）

DEAP（Distributed Evolutionary Algorithms in Python）は、進化計算（Evolutionary Computation）や遺伝的アルゴリズム（Genetic Algorithms）を含む進化的アルゴリズムのためのPythonライブラリです。

以下は、deapライブラリの主な特徴と機能です。

1. 遺伝的アルゴリズムのサポート：

deapは主に遺伝的アルゴリズム（Genetic Algorithms）を実装するためのツールボックスを提供しています。
遺伝的アルゴリズムは、生物学的な進化の原則に基づいて問題を解く最適化手法です。

2. 進化計算アルゴリズム：

deapは他の進化計算アルゴリズムもサポートしており、遺伝的プログラミング（Genetic Programming）、進化戦略（Evolution Strategies）、遺伝的プログラミング（Genetic Programming）などを含む様々な進化計算手法を実装できます。

3. 柔軟性と拡張性：

deapは柔軟かつ拡張性があり、ユーザーが独自の進化的アルゴリズムを構築するための豊富な機能を提供します。
遺伝的操作や適応度評価関数の定義、進化戦略の設定など、各要素を簡単にカスタマイズできます。

4. 統計情報の収集：

deapは進化計算の過程で発生する*統計情報を収集し、進化の様子**を視覚化するためのツールも提供します。
これにより、進化の効率や収束の速さなどを分析することが可能です。

5. 分散進化アルゴリズム：

deapは分散進化計算をサポートし、分散環境での進化計算を実現できます。
これにより、複雑な最適化問題を解決するために複数の計算ノードを利用することが可能です。

サンプルコード

以下は、deapライブラリの基本的な使用例です。

この例は、遺伝的アルゴリズムを使用して1次元の関数の最小値を求めるものです。

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

!pip install deap

from deap import base, creator, tools, algorithms
import random
import numpy as np

# 適応度クラスの作成
creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
creator.create("Individual", np.ndarray, fitness=creator.FitnessMin)

# 遺伝子個体の初期化
toolbox = base.Toolbox()
toolbox.register("attr_float", random.uniform, -5, 5)
toolbox.register("individual", tools.initRepeat, creator.Individual, toolbox.attr_float, 1)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)

# 適応度評価関数
def evaluate(individual):
    return individual[0]**2,

# 遺伝的操作
toolbox.register("mate", tools.cxBlend, alpha=0.5)
toolbox.register("mutate", tools.mutGaussian, mu=0, sigma=1, indpb=0.2)
toolbox.register("select", tools.selTournament, tournsize=3)
toolbox.register("evaluate", evaluate)

# 進化計算の実行
pop = toolbox.population(n=50)
algorithms.eaMuPlusLambda(pop, toolbox, mu=10, lambda_=50, cxpb=0.7, mutpb=0.2, ngen=100, stats=None, halloffame=None, verbose=True)

# 結果の表示
best_ind = tools.selBest(pop, 1)[0]
print("最適な個体:", best_ind)
print("最小値:", best_ind.fitness.values[0])

このコードでは、2乗した値が最小となるような1次元の関数を最適化しています。

[実行結果]

gen    nevals
0      50    
1      45    
2      40    
3      43    
4      45    
5      41    
6      42    
7      45    
8      47    
9      43    
10     43    
11     46    
12     42    
13     49    
14     47    
15     47    
16     43    
17     46    
18     46    
19     45    
20     43    
21     46    
22     47    
23     47    
24     45    
25     44    
26     45    
27     42    
28     41    
29     48    
30     45    
31     42    
32     47    
33     45    
34     45    
35     46    
36     45    
37     47    
38     47    
39     46    
40     44    
41     42    
42     46    
43     39    
44     43    
45     49    
46     45    
47     45    
48     40    
49     46    
50     43    
51     47    
52     45    
53     46    
54     46    
55     46    
56     46    
57     45    
58     44    
59     43    
60     48    
61     48    
62     46    
63     42    
64     46    
65     49    
66     40    
67     44    
68     47    
69     48    
70     48    
71     44    
72     45    
73     45    
74     44    
75     46    
76     44    
77     48    
78     48    
79     43    
80     42    
81     43    
82     42    
83     48    
84     41    
85     42    
86     47    
87     45    
88     45    
89     43    
90     43    
91     43    
92     45    
93     47    
94     47    
95     43    
96     44    
97     48    
98     45    
99     45    
100    44    
最適な個体: [-2.50777571e-07]
最小値: 6.288939021627859e-14

ソースコード解説

このPythonスクリプトは、deapライブラリを使用して進化計算（遺伝的アルゴリズム）を実行し、1次元の関数を最適化する例です。

以下にコードの章立てと詳細な説明を示します。

1. モジュールのインポート:

1
2
3

from deap import base, creator, tools, algorithms
import random
import numpy as np

必要なモジュールやライブラリをインポートしています。
deapライブラリに加え、randomやnumpyも使用します。

2. 適応度クラスの作成:

1
2

creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
creator.create("Individual", np.ndarray, fitness=creator.FitnessMin)

deapで遺伝的アルゴリズムを使用するために必要な適応度クラスと個体クラスを作成しています。

3. 遺伝子個体の初期化:

1
2
3
4

toolbox = base.Toolbox()
toolbox.register("attr_float", random.uniform, -5, 5)
toolbox.register("individual", tools.initRepeat, creator.Individual, toolbox.attr_float, 1)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)

遺伝子個体の初期化を行っています。
個体は1次元の浮動小数点数からなり、その範囲は$-5$から$5$までです。

4. 適応度評価関数:

1
2

def evaluate(individual):
   return individual[0]**2,

適応度関数を定義しています。
この例では、最小値を求めるために個体の1乗を適応度としています。

5. 遺伝的操作:

1
2
3
4

toolbox.register("mate", tools.cxBlend, alpha=0.5)
toolbox.register("mutate", tools.mutGaussian, mu=0, sigma=1, indpb=0.2)
toolbox.register("select", tools.selTournament, tournsize=3)
toolbox.register("evaluate", evaluate)

遺伝的操作（交叉、突然変異、選択）を定義しています。
具体的な操作は、ブレンド交叉、ガウシアン突然変異、トーナメント選択です。

6. 進化計算の実行:

1
2

pop = toolbox.population(n=50)
algorithms.eaMuPlusLambda(pop, toolbox, mu=10, lambda_=50, cxpb=0.7, mutpb=0.2, ngen=100, stats=None, halloffame=None, verbose=True)

初期個体を生成し、進化計算を実行しています。
eaMuPlusLambdaはμ+λ進化アルゴリズムを指します。

7. 結果の表示:

1
2
3

best_ind = tools.selBest(pop, 1)[0]
print("最適な個体:", best_ind)
print("最小値:", best_ind.fitness.values[0])

進化計算の結果として、最適な個体とその適応度を表示しています。

適応度関数の最小値が最適な解です。

結果解説

実行結果は、進化計算の各世代（generation）ごとの統計情報を示しています。

各行には、進化計算の特定の世代における情報が含まれています。

• gen: 世代数
• nevals: 評価された個体数

例えば、最初の行では世代数$0$で$50$個体が評価されました。

最終的には$100$世代まで進化計算が行われています。

ここで注目すべきポイントは、最終的な結果です：

• 最適な個体: [-2.50777571e-07]
• 最小値: 6.288939021627859e-14

最適な個体は1つの次元を持つ浮動小数点数の配列で、その最小値はほぼゼロに非常に近い値です。
進化計算は、遺伝的アルゴリズムを使用してこの最小値に収束しました。

なお、最適な個体がゼロに近い値であることから、この例では1次元関数$ f(x) = x^2 $を最小化しています。
ゼロがこの関数の最小値であり、進化計算がそれに収束したことが示されています。

多項式関数

複雑な数式として、例えば次のような多項式関数を考えてみましょう。

$$
f(x, y) = \sin(x^2 + y^2) \cdot \cos(xy)
$$

このような関数を3Dグラフで表現することができます。

以下は、この関数のグラフをPythonでプロットする例です。

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

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

x = np.linspace(-5, 5, 100)
y = np.linspace(-5, 5, 100)
x, y = np.meshgrid(x, y)
z = np.sin(x**2 + y**2) * np.cos(x * y)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x, y, z, cmap='viridis')

ax.set_xlabel('X-axis')
ax.set_ylabel('Y-axis')
ax.set_zlabel('Z-axis')
ax.set_title('Graph of a Complex Function')

plt.show()

これにより、関数$ f(x, y) = \sin(x^2 + y^2) \cdot \cos(xy) $の3Dグラフが生成されます。

[実行結果]

数式を変更して他の複雑な関数をプロットすることもできます。

ソースコード解説

以下に、コードの各部分を詳しく説明します。

1. numpyライブラリの導入:

1

import numpy as np

数値計算に使用されるNumPyライブラリを導入しています。

2. matplotlib.pyplotおよびmpl_toolkits.mplot3dの導入:

1
2

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

Matplotlibのpyplotモジュールと3DプロットのためのAxes3Dを導入しています。

3. データの生成:

1
2
3
4

x = np.linspace(-5, 5, 100)
y = np.linspace(-5, 5, 100)
x, y = np.meshgrid(x, y)
z = np.sin(x**2 + y**2) * np.cos(x * y)

xとyはそれぞれ$-5$から$5$までの範囲で等間隔に$100$点生成されます。
np.meshgridを使用して2つの1次元配列をグリッドデータに変換し、zは複雑な関数の値が計算されます。
この場合、z = sin(x**2 + y**2) * cos(x * y) です。

4. 3Dプロットの作成:

1
2

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

MatplotlibのFigureオブジェクトと3Dサブプロットを作成します。

5. 3Dサーフェスプロット:

1

ax.plot_surface(x, y, z, cmap='viridis')

ax.plot_surfaceを使用して、x、y、およびzで表される3Dサーフェスプロットを作成します。
cmapパラメータはカラーマップを指定しています。
ここでは’viridis’を使用しています。

6. 軸ラベルとタイトルの設定:

1
2
3
4

ax.set_xlabel('X-axis')
ax.set_ylabel('Y-axis')
ax.set_zlabel('Z-axis')
ax.set_title('Graph of a Complex Function')

X軸、Y軸、およびZ軸のラベルを設定し、グラフにタイトルを追加します。

7. グラフの表示:

1

plt.show()

最後にplt.show()を呼び出して、グラフを表示します。

このコードによって生成されるグラフは、3Dサーフェスプロットで表される関数の形状を示しています。

関数は $x$, $y$ の二乗の和と $x$ と$ y $の積に依存しており、$sin$ および$ cos $関数を組み合わせています。

グラフの形状や色合いは、この複雑な関数の値に基づいています。

クーロンの法則 3Dグラフ

物理学の中でも電磁気学に現れるポテンシャル方程式であるクーロンの法則を題材に、Pythonで3Dグラフ化してみましょう。

クーロンの法則によれば、2つの電荷間の力は各電荷の大きさと距離の二乗に反比例します。

ここでは負電荷の周りの電位分布を3Dで可視化してみましょう。

必要なライブラリをインポートしましょう：

1
2
3

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

次に、電位計算と可視化関数をそれぞれ定義します：

1
2
3
4
5
6
7
8
9
10
11
12
13
14

def potential(x, y, q, pos):
    k = 8.99e9
    return k*q/np.sqrt((x-pos[0])**2 + (y-pos[1])**2)

def plot_potential(q, pos):
    x = np.linspace(-10, 10, 100)
    y = np.linspace(-10, 10, 100)
    x, y = np.meshgrid(x, y)
    z = potential(x, y, q, pos)

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.plot_surface(x, y, z, cmap='jet')
    plt.show()

そして、電荷 $-q$が$-1$、位置が$(0, 0)$の場合の電位をプロットします：

1

plot_potential(-1, [0, 0])

それぞれの関数で、potentialでは電荷 $q$の大きさと位置から電位を計算し、plot_potentialではその電位を3Dグラフでプロットしています。

電荷の位置は$0,0$とし、$x$と$y$の範囲を$-10$から$10$までとしています。

また、電位の計算にはクーロンの定数 $k$を用いています。

[実行結果]

ソースコード解説

このソースコードは、電荷が生じる空間における電位を3Dプロットするためのものです。

以下に、ソースコードの章立てと説明を示します。

1. 必要なライブラリのインポート

1
2
3

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

• NumPy: 数学的な演算を効率的に行うためのライブラリ。
• Matplotlib: グラフ描画のためのライブラリ。
• mpl_toolkits.mplot3d: 3Dプロットを生成するためのツールキット。

2. 電位を計算する関数 potential の定義

1
2
3

def potential(x, y, q, pos):
   k = 8.99e9
   return k*q/np.sqrt((x-pos[0])**2 + (y-pos[1])**2)

• potential 関数は、与えられた座標 (x, y) における電位を計算します。
• q: 電荷の大きさ。
• pos: 電荷の位置。

3. 電位を3Dプロットする関数 plot_potential の定義

1
2
3
4
5
6
7
8
9
10

def plot_potential(q, pos):
   x = np.linspace(-10, 10, 100)
   y = np.linspace(-10, 10, 100)
   x, y = np.meshgrid(x, y)
   z = potential(x, y, q, pos)

   fig = plt.figure()
   ax = fig.add_subplot(111, projection='3d')
   ax.plot_surface(x, y, z, cmap='jet')
   plt.show()

• plot_potential 関数は、指定された電荷 (q) と位置 (pos) に基づいて、空間内の電位を3Dサーフェスプロットとして描画します。
• x と y は座標軸の範囲を表し、これらの座標で x, y 平面上の格子点を生成します。
• potential 関数を使用して各格子点で電位を計算し、それを3Dサーフェスプロットとして描画します。

4. 関数の呼び出し

1

plot_potential(-1, [0, 0])

• -1 の電荷が原点 [0, 0] にある場合の電位をプロットします。

このコードは、電場のポテンシャルの視覚化を行います。

結果解説

[実行結果]

このグラフは、電荷が生じる空間における電位を表しています。

具体的には、電荷が原点 (0, 0) にあり、その電荷が周囲の空間に生じる電位を3Dサーフェスプロットで視覚化しています。

以下はグラフの詳細な説明です：

1. 座標軸と目盛り：

• X軸とY軸はそれぞれ -10 から 10 の範囲で、2次元平面を表します。
• Z軸は電位を表し、色で表示されています。

2. サーフェスプロット：

• 3Dサーフェスプロットは、座標平面上の各点における電位を表しています。
• カラーマップ “jet” が使用されており、色の濃淡が電位の変動を示しています。

3. 電荷の位置と影響範囲：

• 電荷は原点$ (0, 0) $にあります。
  この点から離れるほど電位が変動します。
• グラフ上で、原点の周りに盛り上がりやくぼみが見られ、これが電荷による電位の影響を示しています。

4. 濃淡と電位の関係：

• カラーマップにより、電位が濃い色ほど高いことを示しています。
  つまり、電位が高い領域ではより濃い色調が見られます。

このグラフは、電荷が空間に与える電位の変動を見ることができるものであり、電場の概念を直感的に理解するのに役立ちます。

モジュロ関数

複雑な数式の例として、複素数のモジュロ関数を考えてみましょう。

この関数は、複素数 $z$に対して、$z$の絶対値（モジュロ）を計算します。

複素数$z$は、実部（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

import numpy as np
import matplotlib.pyplot as plt

# 複素数zの範囲を設定
N = 500
lim = 4

# 複素数zの実部と虚部のグリッドを生成
x, y = np.meshgrid(np.linspace(-lim, lim, N), np.linspace(-lim, lim, N))
z = x + 1j*y

# 複素数zのモジュロを計算
f = abs(1/(z**2-4))
f[f > 1.3] = 1.3 # 極（無限大）の近くで関数を切る

# 3Dプロットの設定
fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111, projection="3d")
ax.plot_surface(x, y, f, cmap="viridis", shade=True, alpha=1)
ax.set_xlabel("Re(z)", size=14)
ax.set_ylabel("Im(z)", size=14)
ax.set_title("|f(z)| = |1/(z^2-4)|", size=18, pad=30)

plt.show()

このコードでは、numpyのmeshgrid関数を使用して複素数 $z$の実部と虚部のグリッドを生成し、それを使用して複素数 $z$を計算しています。

そして、複素数 $z$のモジュロを計算し、それを3Dプロットに使用しています。

[実行結果]

ソースコード解説

ソースコードを解説します。

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

numpy と matplotlib を使用します。
numpy は数値計算を行うために使用され、matplotlib はグラフのプロットに使用されます。

2. 複素数zの範囲の設定:

• N = 500: グリッドの数を定義します。
• lim = 4: 複素平面の範囲を設定します。
  実部と虚部の両方が$ -4 $から$ 4 $の範囲で設定されます。

3. 複素数zの実部と虚部のグリッドを生成:

• np.meshgrid() を使って、$ -4 $から$ 4 $までの実部と虚部の範囲でグリッドを生成します。

4. 複素数zのモジュロ計算:

• $ ( f(z) = \left| \frac{1}NaN \right| ) $の計算を行います。
  ここでは abs() を使用して、モジュロ（絶対値）を求めます。
• 極（無限大）の近くで値が非常に大きくなるのを防ぐために、一定の値（ここでは$1.3$）を超える部分を$1.3$にクリッピングしています。

5. 3Dプロットの設定:

• plt.figure() で図を作成し、projection="3d" を指定して3Dプロットを作成します。
• ax.plot_surface() で複素関数の3Dプロットを作成します。
  カラーマップは “viridis” を使っています。
• ax.set_xlabel()、ax.set_ylabel()、ax.set_title() でそれぞれx軸、y軸、タイトルのラベルを設定しています。

6. グラフの表示:

• plt.show() でグラフを表示します。

このコードは、複素関数の特定の領域におけるモジュロの振る舞いを可視化するものです。

結果解説

[実行結果]

上記のPythonコードで生成される3Dグラフは、複素数$z$のモジュロを計算した結果を表示しています。
ここで$z$は複素数で、実部（x座標）と虚部（y座標）から成り立ちます。

具体的には、x軸とy軸は複素数$z$の実部と虚部を表し、z軸は複素数$z$のモジュロを表しています。
つまり、3Dグラフの一つの点は、複素数$z$ （実部と虚部のペア）とそのモジュロを表しています。

また、このグラフは$ |1/(z^2-4)| $の形をしています。

これは、複素数$z$の二乗から4を引いたものの逆数の絶対値を計算したものです。

この関数は、$ z^2-4 $が$0$に等しいとき（つまり、$z$が$±√4=±2のとき）$には無限大となります。

そのため、このグラフでは、$z$が$±2$の近くで関数値が$1.3$に切られています。

複素数関数 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

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

# グリッドを作成
real = np.linspace(-5, 5, 100)
imag = np.linspace(-5, 5, 100)
real, imag = np.meshgrid(real, imag)
z = real + 1j * imag

# 方程式
eqn = z**3 - 1  # 方程式の例

# 3Dグラフの作成
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(np.real(z), np.imag(z), np.real(eqn), c=np.real(eqn), cmap='viridis')

# グラフのラベル設定
ax.set_xlabel('Real axis')
ax.set_ylabel('Imaginary axis')
ax.set_zlabel('Real part of the equation')
ax.set_title('Graph of a complex equation')

plt.show()

この例では、$ (z^3 - 1) $という複素数の3次方程式を使用しています。

これは複素数平面上の方程式です。

3Dプロットは、実部と虚部の関係を示し、方程式の解の振る舞いを可視化します。

[実行結果]

ソースコード解説

このコードは、複素平面上で$ (z^3 - 1) $のグラフを3Dで描画しています。

ここでの$ (z) $は、複素数$ (z = x + yi) $で表され、実部が$ (x) $軸、虚部が$ (y) $軸にマッピングされています。

1. グリッドの作成:

np.linspace を用いて、実数と虚数の範囲を指定し、それぞれ$ (x) $軸および$ (y) $軸の値を生成し、meshgrid を使って複素平面上の格子点を作成しています。

2. 方程式の定義:

z**3 - 1 という方程式を定義し、eqn に格納しています。これは$ (z^3 - 1 = 0) $という方程式を表しています。

3. 3Dグラフの描画:

matplotlib を用いて、3Dの散布図を描画しています。
scatter 関数を使い、複素数$ (z) $の実部と虚部に対応する点をプロットし、色は$ (z^3 - 1) $の実部によって変わります。

4. 軸ラベルとタイトルの設定:

set_xlabel、set_ylabel、set_zlabel を用いてそれぞれ$ (x) $軸、$ (y) $軸、$ (z) $軸のラベルを設定し、set_title でグラフのタイトルを設定しています。

5. グラフの表示:

plt.show() を使って、作成した3Dグラフを表示しています。

このグラフは、複素平面上で$ (z^3 - 1) $の方程式を表し、その実部がゼロになる$ (z) $の値を表しています。

それにより得られる図形は、3次元曲線になります。

結果解説

[実行結果]

このグラフは、複素平面上で方程式$ ( z^3 - 1 ) $を表現しています。

グラフは複素数 $ ( z = x + yi ) $の実部、虚部、および方程式 $ ( z^3 - 1 ) $の実部を3Dで表現しています。

具体的には、複素平面内で実部が$ ( x ) $軸、虚部が$ ( y ) $軸、そして$ ( z^3 - 1 ) $の実部が$ ( z ) $軸にマッピングされています。

グラフ上の点は複素数 $ ( z ) $の実部と虚部に基づいて配置され、色は$ ( z^3 - 1 ) $の実部を示しています。

このグラフは$ ( z^3 - 1 ) $の実部が$ ( 0 ) $になるような$ ( z ) $の値を示しており、その結果得られる図形は特定の形を持つ3次元曲線です。

この方程式が複素平面上でどのように振る舞うかを示しています。

弱肉強食 最適化問題

弱肉強食に関する最適化問題を数学的に表現するのは難しいですが、関連する考え方を最適化問題に当てはめることはできます。

例えば、生態系における捕食者と被食者の関係をモデル化し、最適な捕食率や資源の利用率を求める問題が考えられます。

以下は捕食者と被食者の数を表すロジスティック方程式を利用した簡単なモデルです。

これを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

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

# ロジスティック方程式（捕食者-被食者モデル）
def predator_prey_system(y, t, alpha, beta, delta, gamma):
    x, y = y
    dxdt = alpha * x - beta * x * y
    dydt = delta * x * y - gamma * y
    return [dxdt, dydt]

# パラメータ設定
alpha = 0.1  # 捕食者の増加率
beta = 0.02  # 捕食者と被食者の相互作用率
delta = 0.02  # 被食者の増加率
gamma = 0.1  # 捕食者の減少率

# 初期値
initial_population = [40, 9]  # 捕食者と被食者の初期数

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

# モデルの解を計算
sol = odeint(predator_prey_system, initial_population, t, args=(alpha, beta, delta, gamma))
predator, prey = sol[:, 0], sol[:, 1]

# グラフ化
plt.figure(figsize=(8, 6))
plt.plot(t, predator, label='Predator')
plt.plot(t, prey, label='Prey')
plt.xlabel('Time')
plt.ylabel('Population')
plt.title('Predator-Prey Dynamics')
plt.legend()
plt.grid(True)
plt.show()

このコードは、捕食者（Predator）と被食者（Prey）の数の時間変化を示します。

[実行結果]

捕食者と被食者の相互作用を表すパラメータを調整することで、弱肉強食の関係を模倣することができます。

ソースコード解説

このコードは、捕食者-被食者（Predator-Prey）モデルを用いて、捕食者と被食者の数の時間変化をシミュレーションしています。

これを行うために、次の手順が含まれています：

1. 関数の定義:

predator_prey_system 関数では、捕食者と被食者の数の時間変化を表す微分方程式（ロジスティック方程式）が定義されています。
この方程式は、捕食者と被食者の増減率、相互作用率に基づいて、それぞれの数の変化を記述します。

2. パラメータ設定:

モデル内の各パラメータ（捕食者の増加率、相互作用率など）が定義され、モデルに適用されます。
これらのパラメータは生態系内での捕食者と被食者の相互作用を表現します。

3. 初期値の設定:

捕食者と被食者の初期値が設定されます。
この値は、シミュレーションの開始時点での捕食者と被食者の数を表します。

4. 時間の定義:

シミュレーションする時間軸が定義されます。
この場合、$0$から$200$までの時間を$1000$ステップで区切っています。

5. 微分方程式の数値解析:

odeint を使用して微分方程式を解きます。
初期値、時間、および定義した微分方程式の関数が渡され、解析された結果が sol に格納されます。

6. グラフ化:

得られた結果をグラフに描画します。
時間に対する捕食者と被食者の数がそれぞれ線グラフとして表示され、どのように増減するかが視覚化されます。

このコードを実行することで、捕食者と被食者の数の時間変化を表すグラフが表示されます。

時間が経過するにつれて、捕食者と被食者の数がどのように変化するかが可視化されます。

結果解説

[実行結果]

このグラフは、捕食者（Predator）と被食者（Prey）の数の時間変化を示しています。
捕食者と被食者の関係を表すロジスティック方程式を用いて、彼らの相互作用をモデル化しました。

• 初期値に基づき、捕食者の数は増加し、被食者の数もそれに応じて増加します。
• しかし、被食者の増加に伴って捕食者の数が増え、それによって被食者の数が減少します。
• その後、捕食者の数が減少すると被食者の数が再び増加します。

このように、捕食者と被食者の数は時間とともに振動する関係を持ちます。

捕食者と被食者の相互作用率や増加率などのパラメータを調整すると、この関係性がどのように変化するかを観察できます。
このようなモデルは、生態系における捕食者と被食者の関係性やそのダイナミクスを理解するために使用されます。

栄養学 最適化問題

栄養学の最適化問題の一例として、食事計画を最適化する問題を考えてみましょう。

例えば、特定の栄養素を最大化または最小化するという目標を持った食事プランを作成することができます。

以下のコードは、特定の栄養素（タンパク質、脂質、炭水化物など）を最大化する食事プランを最適化する問題を扱います。

SciPyのminimize関数を使用します。

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

from scipy.optimize import minimize

# 目的関数：最大化する栄養素
def objective_function(x):
    # 栄養素の吸収率（例としてランダムな値を使用）
    protein_absorption = 0.8
    fat_absorption = 0.7
    carb_absorption = 0.9
    
    # 栄養素の含有量（仮の値を使用）
    protein_content = x[0]
    fat_content = x[1]
    carb_content = x[2]
    
    # 含有する栄養素の量に対して吸収率を乗算して目的関数を設定
    return -(protein_content * protein_absorption + fat_content * fat_absorption + carb_content * carb_absorption)

# 初期推定値
initial_guess = [50, 80, 200]  # タンパク質、脂質、炭水化物の初期量

# 制約条件
constraints = ({'type': 'ineq', 'fun': lambda x: 200 - sum(x)},  # 合計カロリー制約
               {'type': 'ineq', 'fun': lambda x: x[0] - 50},     # タンパク質最低量制約
               {'type': 'ineq', 'fun': lambda x: x[1] - 60},     # 脂質最低量制約
               {'type': 'ineq', 'fun': lambda x: x[2] - 150})    # 炭水化物最低量制約

# 最適化の実行
result = minimize(objective_function, initial_guess, constraints=constraints, method='COBYLA')

# 最適化結果
optimal_meal_plan = result.x
maximized_nutrient = -result.fun

print(f"Optimal Meal Plan (Protein, Fat, Carbs): {optimal_meal_plan}")
print(f"Maximized Nutrient: {maximized_nutrient}")

このコードは、最大化したい栄養素を設定し、各栄養素の含有量を変数として最適化します。

制約条件を設定し、合計カロリーと各栄養素の最低摂取量を制限しています。

結果として、最適な食事プランと最大化された栄養素量が出力されます。

[実行結果]

Optimal Meal Plan (Protein, Fat, Carbs): [ 35.  45. 135.]
Maximized Nutrient: 181.0

ソースコード解説

このコードは、栄養学の観点から最適な食事プランを見つける最適化問題を扱います。

1. 目的関数の定義:

• objective_function関数は、最大化したい栄養素（タンパク質、脂質、炭水化物など）を表しています。
• 各栄養素の含有量とそれぞれの吸収率を考慮し、それらの積の総和を最大化することがこの関数の目的です。

2. 初期推定値:

• initial_guessは、最適化アルゴリズムによって使用される栄養素の初期量を示しています。

3. 制約条件:

• constraintsでは、合計カロリーと各栄養素の最低摂取量に関する制約条件を設定しています。
• この場合、合計カロリーは200未満でなければならず、タンパク質、脂質、炭水化物はそれぞれ50g、60g、150g以上である必要があります。

4. 最適化の実行:

• minimize関数を使用して、定義した目的関数と制約条件を考慮して最適な栄養素を見つけます。
• ここではCOBYLA法を使用して最適化を行っています。

5. 最適化結果の表示:

• 最適な食事プラン（タンパク質、脂質、炭水化物の含有量）と、それによって最大化された栄養素の量が出力されます。

結果解説

この結果は、特定の栄養素（タンパク質、脂質、炭水化物など）を最大化する食事プランを見つけたことを示しています。

[実行結果]

Optimal Meal Plan (Protein, Fat, Carbs): [ 35.  45. 135.]
Maximized Nutrient: 181.0

• 最適な食事プラン:
  最適な食事プランでは、タンパク質が35g、脂質が45g、炭水化物が135g含まれています。
  これは、目標としていた栄養素を最大化するために、それぞれの栄養素が含まれる量を示しています。
  これらの値は最適化アルゴリズムが見つけた最良のバランスを表しています。

• 最大化された栄養素量:
  最大化された栄養素の値は181.0です。
  この値は、目的関数が最大化した栄養素の総量を示しています。
  最適化の目標は、この数値をできるだけ高くすることでした。

宇宙工学 最適化問題

宇宙工学の最適化問題の一例として、宇宙船の軌道設計を考えてみましょう。

例えば、地球周回軌道への衛星の打ち上げを考えます。

具体的には、衛星の打ち上げコストと到達する軌道の高度、または到達するまでの時間などの複数の要素を最適化する問題があります。

ここでは、簡単な例として、地球周回軌道への衛星の打ち上げコストと軌道の高度の関係を考えてみます。

軌道の高度が増すにつれて打ち上げコストが増えると仮定します。

以下のコードは、この問題を最適化するためにScipyライブラリを使用しています。

具体的には、打ち上げコストを最小化するための最適な軌道高度を探します。

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

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

# 目的関数：打ち上げコストの最小化
def launch_cost(height):
    # 仮想的な打ち上げコスト関数（ここでは単純な2次関数）
    return (height - 500) ** 2  # 高度と500kmの差の二乗をコストと仮定

# 初期推定値
initial_guess = 1000  # 初期の軌道高度の推定値

# 最適化の実行
result = minimize(launch_cost, initial_guess, method='BFGS')

# 最適化結果
optimal_height = result.x
min_launch_cost = result.fun

# グラフ化
heights = np.linspace(400, 600, 100)  # 400kmから600kmまでの高度を考慮
costs = (heights - 500) ** 2  # 打ち上げコスト関数

plt.plot(heights, costs, label='Launch Cost')
plt.scatter(optimal_height, min_launch_cost, color='red', label='Optimal Height')
plt.xlabel('Orbital Height (km)')
plt.ylabel('Launch Cost')
plt.title('Optimization of Satellite Launch Height')
plt.legend()
plt.grid(True)
plt.show()

print(f"Optimal Height for Minimum Launch Cost: {optimal_height} km")
print(f"Minimum Launch Cost: {min_launch_cost}")

このコードは、打ち上げコストを最小化するための最適な軌道高度を見つけ、その結果をグラフ化しています。

グラフでは、軌道高度に対する打ち上げコストの関係が示され、最適な軌道高度が赤い点として示されます。

[実行結果]

ソースコード解説

このコードは、宇宙工学における最適化問題の一例を示しています。

以下がソースコードの詳細です。

1. 目的関数の定義:

• launch_cost関数は、仮想的な打ち上げコスト関数を定義しています。
  ここでは、高度と500kmの差の二乗をコストとしています。

2. 初期推定値の設定:

• initial_guessは、最適化アルゴリズムの初期の推定値として1000kmを設定しています。

3. 最適化の実行:

• minimize関数を使って、launch_cost関数を最小化することで最適な軌道高度を探します。
  ここではBFGS法を用いて最適化を行っています。

4. 最適化結果の取得:

• 最適な軌道高度とその時の打ち上げコストを取得しています。

5. グラフ化:

• heightsは400kmから600kmまでの範囲の高度を考慮するための配列です。
• costsはそれぞれの高度における打ち上げコストを示します。
• plt.plot()で高度に対する打ち上げコストの関係を表すグラフをプロットしています。
• plt.scatter()で最小の打ち上げコストを達成する最適な高度を赤い点として表示しています。
• その他、グラフのタイトルや軸ラベル、凡例、グリッドを設定しています。

6. 結果の表示:

• 最適な軌道高度とその時の最小打ち上げコストを出力しています。

結果解説

[実行結果]

この結果は、打ち上げコストを最小化するための最適な軌道高度を見つけたことを示しています。
結果として得られた最適な軌道高度は約500kmです。
しかし、数値計算の誤差のため、厳密には499.99999999kmという値が得られています。

さらに、この最適な高度での打ち上げコストは非常に小さく、ほぼゼロに近い値となっています（6.48426847315948e-17）。
これは浮動小数点の誤差などによるもので、ほぼゼロに等しい非常に小さな値です。

グラフでは、軌道高度が高くなるにつれて打ち上げコストが増加することが示されています。
そして、その関係を表す2次関数的なカーブが描かれており、最小値が500km付近であることが視覚的に示されています。
赤い点は最適な高度での最小コストを示しています。

リカードの比較優位

貿易理論には、複数の方程式やモデルが存在しますが、最も基本的なものの一つにリカードの比較優位に関連する方程式があります。

リカードの比較優位理論は、異なる国が特定の商品を生産する際の優位性を分析します。

これには特に次の数式が関連しています。

ここで、$ (P_x) $と$ (P_y) $はそれぞれ商品$ (X) $と商品$ (Y) $の価格を示し、$ (A_x) $と$ (A_y) $は生産コストまたは生産に要する労働時間などを示します。

この不等式が成り立つ場合、つまり商品$ (X) $の国内価格が国際価格よりも低いとき、その国は商品$ (X) $の生産に比較的適しており、他国との間で特化と貿易を行うことで両国とも利益を最大化できる可能性が高くなります。

この方程式はリカードの比較優位を示すものであり、貿易理論における基本的なモデルの一部です。

例題

リカードの比較優位に関連した最適化問題の一例として、2つの国（国Aと国B）が2つの商品（商品Xと商品Y）を生産する場合を考えてみましょう。

以下はその問題設定です。

国Aで商品Xの生産に1時間かかり、商品Yの生産に2時間かかるとします。
国Bでは商品Xの生産に3時間かかり、商品Yの生産に1時間かかるとします。
また、国際価格における商品Xと商品Yの価格比率が2:1であるとします。

この場合、どの国がどの商品を生産すべきか、そして両国間での貿易が有益かを最適化する問題を考えることができます。

以下はPythonでこの問題を解くコードです。
ここではScipyライブラリを使って最適化を行い、結果をグラフ化します。

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

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

# 各国の生産にかかる時間（単位：時間）
# [国Aの商品X生産時間, 国Aの商品Y生産時間]
# [国Bの商品X生産時間, 国Bの商品Y生産時間]
production_times = np.array([[1, 2], [3, 1]])

# 国際価格比率（商品X価格 / 商品Y価格）
international_price_ratio = 2

# 目的関数：貿易利益の最大化（国Aと国Bの生産量）
def objective_function(production_quantities):
    trade_profit = -np.sum(production_quantities * production_times)
    return trade_profit

# 制約条件：国際価格比率と生産量の関係
def constraint_function(production_quantities):
    return international_price_ratio - production_quantities[0] / production_quantities[1]

# 初期推定値
initial_guess = [1, 1]

# 最適化の実行
constraint = {'type': 'eq', 'fun': constraint_function}
result = minimize(objective_function, initial_guess, constraints=constraint)

# 最適化結果
optimal_production_quantities = result.x
trade_profit = -result.fun

# グラフ化
labels = ['国Aの商品X', '国Aの商品Y']
fig, ax = plt.subplots()
ax.bar(labels, optimal_production_quantities)
ax.set_ylabel('生産量')
ax.set_title('最適生産量')
plt.show()

print(f"最適な生産量: {optimal_production_quantities}")
print(f"貿易利益: {trade_profit}")

このコードは、Scipyのminimizeを使って貿易利益を最大化する生産量を見つけ、結果を棒グラフで表示します。

最適な生産量と貿易利益が表示されます。

ただし、このコードは簡略化されたモデルですので、実際の貿易問題には多くの要因が含まれる場合があります。

[実行結果]

ソースコード解説

このコードは、2つの国（国Aと国B）が2つの商品（商品Xと商品Y）を生産する場合の貿易最適化問題を解くためのPythonスクリプトです。

以下がコードの詳細です。

1. 各国の生産時間と国際価格比率の設定:

• production_timesは2つの国での商品の生産時間を表す行列です。
  例えば、国Aでは商品Xの生産に1時間かかり、商品Yの生産に2時間かかります。
• international_price_ratioは商品Xと商品Yの価格比率を示します。

2. 目的関数と制約条件の定義:

• objective_functionは貿易利益を最大化するための目的関数です。
  生産量と生産時間を掛け合わせ、その合計をマイナスにしています。
• constraint_functionは国際価格比率と生産量の関係を示す制約条件です。

3. 最適化の実行:

• minimize関数を使って、目的関数を最大化するための最適化を行います。
  初期推定値はinitial_guessで設定されており、制約条件はconstraintに定義されています。

4. 最適化結果の取得:

• result.xは最適な生産量を示しています。
• result.funは最適な目的関数の値で、貿易利益を示しています。

5. グラフ化:

• 最適な生産量を棒グラフとして表示しています。

最終的に、このコードは貿易利益を最大化するための最適な生産量を計算し、それをグラフ化して表示しています。

また、最適な生産量と貿易利益をコンソールに出力しています。

結果解説

[実行結果]

最適な生産量は、国Aが商品Xをおよそ1.37兆単位、商品Yをおよそ0.69兆単位生産することを示しています。
一方、国Bは商品Xをほぼ0単位、商品Yをおよそ0.69兆単位生産することを示しています。

この生産量の結果、貿易利益は755兆USDに相当します。
この利益は、国際価格比率と生産コストの違いを利用して、国Aが商品Xを生産し、国Bが商品Yを生産することで生じる貿易利益を表しています。

しかし、結果が非常に大きな数字になっており、生産量のバランスや制約条件の設定に誤差がある可能性があります。
このような場合、モデルや数値計算の調整が必要になるかもしれません。