自由落下

物体の自由落下を考えたシミュレーションを考えます。

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

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

# 物体の自由落下を表す微分方程式
def free_fall(y, t):
    g = 9.8  # 重力加速度 (m/s^2)
    dydt = [y[1], -g]  # [速度, 加速度]
    return dydt

# 初期条件
y0 = [0, 0]  # 初期位置と初速度

# 時間の範囲
t = np.linspace(0, 5, 100)

# 微分方程式を解く
sol = odeint(free_fall, y0, t)

# 結果をグラフ化
plt.plot(t, sol[:, 0], label='Position (m)')
plt.plot(t, sol[:, 1], label='Velocity (m/s)')
plt.xlabel('Time (s)')
plt.legend()
plt.title('Free Fall Simulation')
plt.show()

この例では、物体の自由落下を表す微分方程式をSciPyのodeint関数を使用して数値的に解き、その結果をグラフ化しています。

[実行結果]

ソースコード解説

このソースコードは、SciPyを使用して物体の自由落下をシミュレーションし、結果をグラフで視覚化するプログラムです。

以下に、ソースコードの詳細な説明を示します。

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

1
2
3

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

• numpy: 数値計算を行うためのライブラリ。
  数学的な演算や配列操作などに使用されます。
• matplotlib.pyplot: グラフの描画に使用されるライブラリ。
• scipy.integrate.odeint: 常微分方程式を数値的に解くためのSciPyライブラリの関数。

自由落下の微分方程式:

1
2
3
4

def free_fall(y, t):
    g = 9.8  # 重力加速度 (m/s^2)
    dydt = [y[1], -g]  # [速度, 加速度]
    return dydt

• free_fall関数: 物体の自由落下を表す微分方程式を定義します。
  • y: 状態ベクトル（位置と速度の組み合わせ）
  • t: 時間
  • g: 重力加速度
  • dydt: 速度と加速度の変化を表すベクトル

初期条件:

1

y0 = [0, 0]  # 初期位置と初速度

• y0: 物体の初期位置と初速度を示す状態ベクトル。

時間の範囲:

1

t = np.linspace(0, 5, 100)

• t: シミュレーションの時間範囲を示す配列。
  $0$から$5$までの時間を$100$等分した配列。

微分方程式の解法:

1

sol = odeint(free_fall, y0, t)

• odeint: SciPyの関数で、常微分方程式を数値的に解く。
  free_fall関数に対して初期条件y0と時間範囲tを与えて解を計算し、solに格納。

結果のグラフ化:

1
2
3
4
5
6

plt.plot(t, sol[:, 0], label='Position (m)')
plt.plot(t, sol[:, 1], label='Velocity (m/s)')
plt.xlabel('Time (s)')
plt.legend()
plt.title('Free Fall Simulation')
plt.show()

• グラフ描画: 解かれた結果をグラフとして描画。
  • sol[:, 0]: 位置の時間変化を表す曲線。
  • sol[:, 1]: 速度の時間変化を表す曲線。
  • xlabel: x軸のラベル（時間）。
  • legend: 曲線の説明を表示。
  • title: グラフのタイトル。
  • show: グラフの表示。

このプログラムは、物体の自由落下を数値的にシミュレーションし、その結果を位置と速度のグラフで視覚化します。

結果解説

[実行結果]

上記のグラフは、物体の自由落下を表すシミュレーションの結果を示しています。

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

1. Position vs. Time (位置と時間の関係)

• x軸（横軸）: 時間（秒）
• y軸（縦軸）: 物体の位置（メートル）
• 線: 時間に対する物体の位置の変化を表現しています。

グラフは時間の経過に伴って物体がどれだけ下方向に移動するかを示しています。
最初はゼロから始まり、時間が進むにつれて増加します。
自由落下の影響で、物体は加速して地面に向かっています。

2. Velocity vs. Time (速度と時間の関係)

• x軸: 時間（秒）
• y軸: 物体の速度（メートル/秒）
• 線: 時間に対する物体の速度の変化を表現しています。

グラフは時間が経つにつれて速度が増加し、これもまた自由落下の影響を反映しています。
最初はゼロから始まり、時間が進むにつれて減速せずに増加していきます。
物体が地面に向かって自由落下するためです。

このシミュレーションでは、物体が重力に従って自由落下する様子を可視化しています。
物体の位置は時間とともに変化し、速度も時間に応じて変動します。

この例は物理学的な現象を数値的にモデリングし、SciPyを使用して微分方程式を解くことで、結果をグラフとして視覚化しています。

放射線源の観測データのクラスタリング

宇宙に関する問題を解決するために、scikit-learnを使用してクラスタリングを行います。

具体的なデータセットとして、宇宙からの放射線源の観測データを考えます。

このデータセットを用いて、放射線源のクラスタを特定し、その結果をグラフ化します。

まず、適切なデータセットを用意しましょう。

宇宙からの放射線源のデータセットは、実際のものではなく、架空のデータセットとして扱います。

以下のサンプルコードでは、make_blobs関数を使用して3つのクラスタを持つデータセットを生成します。

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
from sklearn.datasets import make_blobs
from sklearn.cluster import KMeans

# 仮想の宇宙からの放射線源データセット生成
X, y = make_blobs(n_samples=300, centers=3, random_state=42, cluster_std=1.0)

# KMeansクラスタリングを実行
kmeans = KMeans(n_clusters=3)
kmeans.fit(X)

# クラスタリング結果を取得
labels = kmeans.labels_
centers = kmeans.cluster_centers_

# データとクラスターセンターをプロット
plt.scatter(X[:, 0], X[:, 1], c=labels, cmap='viridis', edgecolor='k', s=50)
plt.scatter(centers[:, 0], centers[:, 1], c='red', marker='X', s=200, label='Cluster Centers')
plt.title('宇宙からの放射線源クラスタリング')
plt.xlabel('放射線源のX座標')
plt.ylabel('放射線源のY座標')
plt.legend()
plt.show()

このサンプルコードでは、make_blobs関数を使用して3つのクラスタを持つデータセットを生成し、KMeansクラスタリングアルゴリズムを適用しています。

結果をプロットすると、各クラスタの中心が赤い”X”で表示され、データ点が色分けされたクラスタに属していることがわかります。

[実行結果]

この例は架空のデータを使用していますが、実際の宇宙データセットを用いることで、クラスタリングアルゴリズムを実際の天文データに適用し、異なる放射線源のグループを発見することができます。

ソースコード解説

このソースコードは、仮想の宇宙からの放射線源データセットを生成し、それにKMeansクラスタリングアルゴリズムを適用して結果を可視化するためのものです。

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

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

1
2
3
4

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_blobs
from sklearn.cluster import KMeans

• numpy：数値計算を行うためのライブラリ
• matplotlib.pyplot：グラフの描画や可視化のためのライブラリ
• make_blobs：クラスタリングのための仮想のデータセットを生成するための関数
• KMeans：KMeansクラスタリングアルゴリズムを提供するscikit-learnのクラス

2. 仮想の宇宙からの放射線源データセット生成

1

X, y = make_blobs(n_samples=300, centers=3, random_state=42, cluster_std=1.0)

• make_blobs関数を使用して、3つのクラスタを持つ仮想の宇宙からの放射線源データセットを生成
• n_samples：データセットのサンプル数
• centers：生成するクラスタの数
• random_state：再現性のための乱数シード
• cluster_std：クラスタの標準偏差

3. KMeansクラスタリングを実行

1
2

kmeans = KMeans(n_clusters=3)
kmeans.fit(X)

• KMeansクラスを用いてクラスタリングモデルを初期化し、fitメソッドでデータに適用
• n_clusters：クラスタの数

4. クラスタリング結果を取得

1
2

labels = kmeans.labels_
centers = kmeans.cluster_centers_

• labels：各データポイントが属するクラスタのラベル
• centers：各クラスタの中心座標

5. データとクラスターセンターをプロット

1
2
3
4
5
6
7

plt.scatter(X[:, 0], X[:, 1], c=labels, cmap='viridis', edgecolor='k', s=50)
plt.scatter(centers[:, 0], centers[:, 1], c='red', marker='X', s=200, label='Cluster Centers')
plt.title('宇宙からの放射線源クラスタリング')
plt.xlabel('放射線源のX座標')
plt.ylabel('放射線源のY座標')
plt.legend()
plt.show()

• plt.scatterを使用してデータポイントを散布図としてプロット
• c=labelsによりクラスタごとに色分けし、cmap='viridis'でカラーマップを指定
• plt.scatterでクラスターセンターを赤い”X”でプロット
• グラフのタイトル、軸ラベル、凡例を設定
• plt.show()でグラフを表示

このソースコードは、データの生成からKMeansクラスタリングの実行、結果の可視化までを一貫して実施しています。

最終的なグラフでは、異なるクラスタのデータポイントと各クラスタの中心が視覚的に表示され、クラスタリングの効果が確認できます。

結果解説

[実行結果]

上記のグラフは、宇宙からの放射線源のクラスタリング結果を可視化したものです。

以下に各要素の説明を記載します。

1. データポイント:

• 散布図上に散らばる点は、仮想の宇宙からの放射線源を表しています。
  各点はデータセット内の1つの観測データを表し、その座標は宇宙の座標を模しています。

2. クラスタ:

• 色分けされた領域は、KMeansクラスタリングアルゴリズムによって同じクラスタに分類されたデータポイントを示しています。
  この例では3つのクラスタがあり、それぞれ異なる色で表されています。
  各クラスタは同じ色で示され、その色はcmap='viridis'によって指定されています。

3. クラスターセンター:

• 赤い”X”は各クラスタの中心を示しています。
  これらの点は、KMeansアルゴリズムによって見つかったクラスタの中心座標です。

4. グラフタイトルと軸ラベル:

• グラフのタイトルは宇宙からの放射線源クラスタリングとなっており、軸ラベルはそれぞれ放射線源のX座標と放射線源のY座標です。

このグラフを見ると、KMeansクラスタリングがデータセットを3つの異なるクラスタに効果的に分類していることがわかります。

各クラスタの中心（赤い”X”）は、そのクラスタのデータポイントが集まる中心点を表しています。

このような可視化を通じて、異なる放射線源のグループがどのように分布しているかを把握することができます。

ターミナル速度方程式

ロケット工学に関連する方程式の一例として、ロケットの運動方程式であるターミナル速度方程式を挙げて、それをPythonで3Dグラフ化してみましょう。

ターミナル速度方程式：

$$
v_t = \sqrt{\frac{2 \cdot g \cdot I_{sp} \cdot \ln\left(\frac{m_0}{m_f}\right)}{A \cdot \rho_c \cdot C_d}}
$$

ここで、各変数の意味は以下の通りです：

• (v_t)：ターミナル速度（Terminal Velocity）
• (g)：重力加速度（通常は地球の場合は約9.8 m/s(^2)）
• (I_{sp})：比推力（Specific Impulse）
• (m_0)：初期質量
• (m_f)：最終質量
• (A)：断面積
• (\rho_c)：大気密度
• (C_d)：抗力係数

以下は、この方程式を3Dグラフ化するための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

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

# ターミナル速度方程式
def terminal_velocity(m0, mf, Isp, A, rho_c, Cd):
    g = 9.8  # 重力加速度 (m/s^2)
    return np.sqrt((2 * g * Isp * np.log(m0 / mf)) / (A * rho_c * Cd))

# パラメータの設定（例）
m0 = 1000  # 初期質量 (kg)
mf = 500   # 最終質量 (kg)
Isp = 300  # 比推力 (s)
A = 5      # 断面積 (m^2)
rho_c = 1  # 大気密度 (kg/m^3)
Cd = 0.5   # 抗力係数

# 初期質量と最終質量の範囲を設定
m0_values = np.linspace(800, 1200, 50)
mf_values = np.linspace(400, 600, 50)

# 3Dグリッドの生成
m0_grid, mf_grid = np.meshgrid(m0_values, mf_values)

# ターミナル速度の計算
vt_values = terminal_velocity(m0_grid, mf_grid, Isp, A, rho_c, Cd)

# 3Dプロット
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(m0_grid, mf_grid, vt_values, cmap='viridis')

# 軸ラベルとタイトルの設定
ax.set_xlabel('Initial Mass (kg)')
ax.set_ylabel('Final Mass (kg)')
ax.set_zlabel('Terminal Velocity (m/s)')
ax.set_title('Terminal Velocity vs. Initial and Final Mass')

plt.show()

この例では、初期質量と最終質量を変化させたときのターミナル速度を3Dグラフで表示しています。

[実行結果]

パラメータや方程式を実際の状況に合わせて変更することが重要です。

ソースコード解説

以下に、コードの主な要素を章立てして説明します。

1. 関数 terminal_velocity の定義:

1
2
3

def terminal_velocity(m0, mf, Isp, A, rho_c, Cd):
   g = 9.8  # 重力加速度 (m/s^2)
   return np.sqrt((2 * g * Isp * np.log(m0 / mf)) / (A * rho_c * Cd))

• ロケットのターミナル速度を計算する関数です。
  関数は初期質量 m0、最終質量 mf、比推力 Isp、断面積 A、大気密度 rho_c、抗力係数 Cdを引数に取ります。

2. パラメータの設定:

1
2
3
4
5
6

m0 = 1000  # 初期質量 (kg)
mf = 500   # 最終質量 (kg)
Isp = 300  # 比推力 (s)
A = 5      # 断面積 (m^2)
rho_c = 1  # 大気密度 (kg/m^3)
Cd = 0.5   # 抗力係数

• ロケットと環境に関するパラメータを設定します。

3. 初期質量と最終質量の範囲の設定:

1
2

m0_values = np.linspace(800, 1200, 50)
mf_values = np.linspace(400, 600, 50)

• グラフに表示する初期質量と最終質量の範囲を設定します。

4. 3Dグリッドの生成:

1

m0_grid, mf_grid = np.meshgrid(m0_values, mf_values)

• 3Dグラフの横軸と縦軸に対応する初期質量と最終質量の組み合わせを作成します。

5. ターミナル速度の計算:

1

vt_values = terminal_velocity(m0_grid, mf_grid, Isp, A, rho_c, Cd)

• 3Dグリッド上の各点でのターミナル速度を計算します。

6. 3Dプロット:

1
2
3

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(m0_grid, mf_grid, vt_values, cmap='viridis')

• matplotlib を使用して、初期質量、最終質量、および対応するターミナル速度の3Dプロットを作成します。

7. グラフの装飾:

1
2
3
4

ax.set_xlabel('Initial Mass (kg)')
ax.set_ylabel('Final Mass (kg)')
ax.set_zlabel('Terminal Velocity (m/s)')
ax.set_title('Terminal Velocity vs. Initial and Final Mass')

• グラフに軸ラベルやタイトルを追加します。

8. グラフの表示:

1

plt.show()

• 最終的に作成された3Dグラフを表示します。

結果解説

[実行結果]

この3Dグラフは、ロケットの運動方程式から計算されるターミナル速度を初期質量と最終質量の両方に依存して表示しています。

以下にグラフの主な要素を説明します。

X軸（Initial Mass）:

グラフの底面にあたり、ロケットの初期質量を表します。
この軸を変化させると、グラフの左右方向が変わります。

Y軸（Final Mass）:

グラフの底面にあたり、ロケットの最終質量を表します。
この軸を変化させると、グラフの前後方向が変わります。

Z軸（Terminal Velocity）:

グラフの縦方向に伸びる軸で、計算されたターミナル速度を表します。
この軸を変化させると、グラフの上下方向が変わります。

グラフの表面は、各点でのターミナル速度を色で表しており、色が濃いほど速度が大きいことを示しています。

このような3Dグラフを通じて、初期質量と最終質量がターミナル速度に与える影響を視覚的に理解できます。

複雑な数式

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

$$
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）の数の時間変化を示しています。
捕食者と被食者の関係を表すロジスティック方程式を用いて、彼らの相互作用をモデル化しました。

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

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

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