生産計画問題

PuLP は Python の線形計画問題を解くためのライブラリです。

ここで示す一例として、生産計画問題を解いてみましょう。

生産計画の最適化はビジネスの現場でよく見る問題の一つで、リソースが限られている中で利益を最大化する生産量を決定します。

問題設定:

• 2つの製品（AとB）を生産します。
• 製品Aと製品Bの利益はそれぞれ20ドル、30ドルです。
• 製品Aを作るのには1時間の作業が、製品Bを作るのには2時間の作業が必要です。
• 製造機械が一週間に稼働できるのは100時間です。
• 市場の需要は製品Aが40個、製品Bが30個までです。
• 製品ごとの生産上限は製品Aは30個、製品Bは20個です。

この情報に基づいて、週の利益を最大化する生産計画を立ててみましょう。

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

import pulp
import matplotlib.pyplot as plt

# 問題のインスタンスを作成します。最大化問題として定義します。
model = pulp.LpProblem("最大利益問題", pulp.LpMaximize)

# 変数を定義します。生産量は0個以上、上限まで可能です。
x1 = pulp.LpVariable('x1', lowBound=0, upBound=30, cat='Continuous')
x2 = pulp.LpVariable('x2', lowBound=0, upBound=20, cat='Continuous')

# 目的関数を定義します。利益を最大化したいので、製品ごとの利益に生産量を掛けたものの和です。
model += 20 * x1 + 30 * x2

# 制約を定義します。
model += x1 + 2 * x2 <= 100  # 機械の作業時間の制限
model += x1 <= 40  # 製品Aの市場需要
model += x2 <= 30  # 製品Bの市場需要

# 問題を解きます。
model.solve()

このコードは PuLP を使用して問題をモデル化し、解を求めます。

そして解（生産量の最適値）を出力します。

1
2
3
4

# 結果の出力
print(f"製品Aの生産量: {x1.value()}")
print(f"製品Bの生産量: {x2.value()}")
print(f"総利益: {pulp.value(model.objective)}")

[実行結果]

製品Aの生産量: 30.0
製品Bの生産量: 20.0
総利益: 1200.0

結果は、与えられた条件と制約のもとで、利益を最大化する生産量の組み合わせを表しています。

ソースコード解説

このソースコードは、線形計画法（Linear Programming）を用いて最大利益を達成する製造計画を最適化する問題を解いています。

以下はコードの章立てと説明です：

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

1
2

import pulp
import matplotlib.pyplot as plt

必要なライブラリ（PuLPとMatplotlib）をインポートします。

2. 問題のインスタンス作成

1

model = pulp.LpProblem("最大利益問題", pulp.LpMaximize)

PuLPを使用して最大利益を求める最適化問題のインスタンスを作成します。

3. 変数の定義

1
2

x1 = pulp.LpVariable('x1', lowBound=0, upBound=30, cat='Continuous')
x2 = pulp.LpVariable('x2', lowBound=0, upBound=20, cat='Continuous')

変数 x1 と x2 を定義し、それぞれが0以上30以下、0以上20以下の実数であることを指定します。

4. 目的関数の定義

1

model += 20 * x1 + 30 * x2

最大化したい目的関数を定義します。

この場合、製品Aと製品Bの単位あたりの利益を最大化します。

5. 制約条件の定義

1
2
3

model += x1 + 2 * x2 <= 100  # 機械の作業時間の制限
model += x1 <= 40  # 製品Aの市場需要
model += x2 <= 30  # 製品Bの市場需要

最適化の際に考慮する制約条件を定義します。

これらは機械の作業時間、製品Aおよび製品Bの市場需要です。

6. 最適化の実行

1

model.solve()

PuLPを使用して最適化問題を解きます。

7. 結果の出力

1
2
3

print(f"製品Aの生産量: {x1.value()}")
print(f"製品Bの生産量: {x2.value()}")
print(f"総利益: {pulp.value(model.objective)}")

最適な製造計画が見つかった後、生産量と総利益を表示します。

ロジスティック回帰モデル

scikit-learnは、Pythonで機械学習モデルを作成するための人気のあるライブラリです。

以下に、scikit-learnを使用してロジスティック回帰モデルを作成し、結果をグラフ化する基本的な手順を示します。

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

1
2
3
4
5
6

from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score, confusion_matrix
import matplotlib.pyplot as plt
import seaborn as sns

次に、データセットをロードし、トレーニングセットとテストセットに分割します：

1
2

dataset = load_breast_cancer()
X_train, X_test, y_train, y_test = train_test_split(dataset.data, dataset.target, test_size=0.2, random_state=42)

次に、ロジスティック回帰モデルを作成し、トレーニングデータに対してフィットします：

1
2

model = LogisticRegression()
model.fit(X_train, y_train)

モデルを使用して予測を行い、予測の精度を計算します：

1
2
3

y_pred = model.predict(X_test)
accuracy = accuracy_score(y_test, y_pred)
print("Accuracy: ", accuracy)

最後に、混同行列を作成し、ヒートマップを使用して視覚化します：

1
2
3
4
5
6

cm = confusion_matrix(y_test, y_pred)
plt.figure(figsize=(10,7))
sns.heatmap(cm, annot=True, fmt='d')
plt.xlabel('Predicted')
plt.ylabel('Truth')
plt.show()

このコードは、scikit-learnの基本的な使用方法を示しています。

まず、データセットをロードし、トレーニングセットとテストセットに分割します。

次に、ロジスティック回帰モデルを作成し、トレーニングデータに対してフィットします。

モデルを使用して予測を行い、予測の精度を計算します。

最後に、混同行列を作成し、ヒートマップを使用して視覚化します。

[実行結果]

ソースコード解説

このソースコードは、乳がんデータセットを用いてロジスティック回帰モデルを構築し、その性能を評価するプログラムです。

以下に、詳細な説明を示します。

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

1
2
3
4
5
6

from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score, confusion_matrix
import matplotlib.pyplot as plt
import seaborn as sns

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

• load_breast_cancer: scikit-learnに含まれる乳がんデータセットをロードするための関数。
• train_test_split: データをトレーニングセットとテストセットに分割する関数。
• LogisticRegression: ロジスティック回帰モデルを構築するためのクラス。
• accuracy_score: 分類モデルの精度を計算するための関数。
• confusion_matrix: 混同行列を計算するための関数。
• matplotlib.pyplot: プロットのためのMatplotlibのモジュール。
• seaborn: データ可視化のためのSeabornライブラリ。

2. データセットの読み込みと分割

1
2

dataset = load_breast_cancer()
X_train, X_test, y_train, y_test = train_test_split(dataset.data, dataset.target, test_size=0.2, random_state=42)

乳がんデータセットを読み込み、説明変数 X と目的変数 y に分割します。
トレーニングセットとテストセットは train_test_split 関数を使用して分割されます。

3. ロジスティック回帰モデルの構築と訓練

1
2

model = LogisticRegression()
model.fit(X_train, y_train)

LogisticRegression クラスを用いてロジスティック回帰モデルを構築し、トレーニングデータを使用してモデルを訓練します。

4. テストデータでの予測と精度の計算

1
2
3

y_pred = model.predict(X_test)
accuracy = accuracy_score(y_test, y_pred)
print("Accuracy: ", accuracy)

構築したモデルを用いてテストデータで予測を行い、その予測精度を accuracy_score 関数で計算し、表示しています。

5. 混同行列の可視化

1
2
3
4
5
6

cm = confusion_matrix(y_test, y_pred)
plt.figure(figsize=(10,7))
sns.heatmap(cm, annot=True, fmt='d')
plt.xlabel('Predicted')
plt.ylabel('Truth')
plt.show()

confusion_matrix 関数を使用して混同行列を計算し、SeabornとMatplotlibを使用して混同行列をヒートマップとして可視化しています。

混同行列はモデルの性能を評価するために用いられます。

結果解説

[実行結果]

混同行列（Confusion Matrix）は、分類モデルのパフォーマンスを評価するための重要なツールです。

混同行列は、モデルの予測結果と実際のラベルを比較し、正確な予測（True Positive, True Negative）と誤った予測（False Positive, False Negative）の数を表示します。

以下に、混同行列の各要素の意味を示します：

True Positive (TP):

モデルが正しく正のクラスを予測したインスタンスの数。

True Negative (TN):

モデルが正しく負のクラスを予測したインスタンスの数。

False Positive (FP):

モデルが正しく負のクラスを予測したが、実際には正のクラスであったインスタンスの数。

False Negative (FN):

モデルが正しく正のクラスを予測したが、実際には負のクラスであったインスタンスの数。

これらの要素を使用して、モデルの精度（Accuracy）や精度（Precision）、再現率（Recall）、F1スコアなどのメトリクスを計算することができます。

また、混同行列はヒートマップとして視覚化することも可能です。

ヒートマップは、各要素の値が色で表現され、高い値を赤色で、低い値を青色で表示します。

これにより、モデルのパフォーマンスを直感的に理解することができます。

非線形／多変数方程式

以下のような関数を考えてみましょう。

この関数は非線形で多変数の方程式であり、解析的な解法は難しいとされています。

$$
f(x, y) = \sin(\sqrt{x^2 + y^2}) - \frac{x}{2} - \frac{y}{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
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40

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

# 非線形方程式の定義
def equation(vars):
    x, y = vars
    return [np.sin(np.sqrt(x**2 + y**2)) - x/2 - y/3, x + y]

# 初期値の設定
initial_guess = [1, 1]

# 非線形方程式の数値解を求める
solution = fsolve(equation, initial_guess)

print("Numerical Solution:", solution)

# 3Dグラフの描画
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')

# 関数の値を計算してプロット
x_vals = np.linspace(-5, 5, 100)
y_vals = np.linspace(-5, 5, 100)
X, Y = np.meshgrid(x_vals, y_vals)
Z = np.sin(np.sqrt(X**2 + Y**2)) - X/2 - Y/3

ax.plot_surface(X, Y, Z, cmap='viridis', alpha=0.8)
ax.scatter(solution[0], solution[1], np.sin(np.sqrt(solution[0]**2 + solution[1]**2)) - solution[0]/2 - solution[1]/3, color='red', marker='o', s=100, label='Solution')

# グラフの設定
ax.set_title('3D Plot of a Nonlinear Equation')
ax.set_xlabel('X-axis')
ax.set_ylabel('Y-axis')
ax.set_zlabel('Z-axis')
ax.legend()

# グラフの表示
plt.show()

この例では、Scipyのfsolveを使用して非線形方程式の数値解を求め、得られた解を3Dグラフ上にプロットしています。

[実行結果]

ソースコード解説

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

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

1
2
3
4

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

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

• numpyは数値計算を行うために使用されます。
• matplotlib.pyplotはグラフ描画のために使用されます。
• scipy.optimize.fsolveは非線形方程式の数値解を求めるために使用されます。
• mpl_toolkits.mplot3dは3Dグラフ描画のために使用されます。

2. 非線形方程式の定義

1
2
3

def equation(vars):
    x, y = vars
    return [np.sin(np.sqrt(x**2 + y**2)) - x/2 - y/3, x + y]

equation関数では、2つの変数$ (x) $と$ (y) $に対する非線形方程式が定義されています。
この例では、$ [ \sin(\sqrt{x^2 + y^2}) - \frac{x}{2} - \frac{y}{3} ] $および$ (x + y) $という2つの方程式が組み合わさっています。

3. 初期値の設定

1

initial_guess = [1, 1]

非線形方程式の数値解を求める際に使用する初期の予測値を設定しています。

4. 非線形方程式の数値解を求める

1
2

solution = fsolve(equation, initial_guess)
print("Numerical Solution:", solution)

fsolve関数を用いて、equation関数に定義された非線形方程式の数値解を求めます。
解は solution に格納され、その値が表示されます。

5. 3Dグラフの描画

1
2

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

Matplotlibを使用して3Dグラフを描画するための準備を行います。

6. 関数の値を計算してプロット

1
2
3
4
5

x_vals = np.linspace(-5, 5, 100)
y_vals = np.linspace(-5, 5, 100)
X, Y = np.meshgrid(x_vals, y_vals)
Z = np.sin(np.sqrt(X**2 + Y**2)) - X/2 - Y/3
ax.plot_surface(X, Y, Z, cmap='viridis', alpha=0.8)

指定された範囲内の$ (x) $と$ (y) $の値を生成し、それに対応する関数の値を計算して3Dプロットします。

7. 解の点をプロット

1

ax.scatter(solution[0], solution[1], np.sin(np.sqrt(solution[0]**2 + solution[1]**2)) - solution[0]/2 - solution[1]/3, color='red', marker='o', s=100, label='Solution')

先ほど求めた数値解を3Dプロットに赤い点として追加します。

8. グラフの設定

1
2
3
4
5

ax.set_title('3D Plot of a Nonlinear Equation')
ax.set_xlabel('X-axis')
ax.set_ylabel('Y-axis')
ax.set_zlabel('Z-axis')
ax.legend()

グラフのタイトルや軸のラベルを設定し、凡例を表示します。

9. グラフの表示

1

plt.show()

最後に、Matplotlibを使用して描画したグラフを表示します。

これにより、非線形方程式の数値解とそのグラフが視覚的に確認できます。

結果解説

[実行結果]

表示されたNumerical Solution ([5.17944401, -5.17944401]) は、非線形方程式

$$
f(x, y) = \sin(\sqrt{x^2 + y^2}) - \frac{x}{2} - \frac{y}{3}
$$

の数値解です。

この解は、関数$ (f(x, y)) $がゼロになる点を表しています。
具体的には、$ (x = 5.17944401) $および$ (y = -5.17944401) $で、この点で方程式が満たされています。

グラフの3Dプロットでは、関数$ (f(x, y)) $の曲面を表示しています。
この曲面上で、赤い点が数値解を表しています。
赤い点が曲面と交差する点が方程式の解であり、その座標が表示された数値解と一致しています。

数値解を求めるプロセスは、初期の予測値（initial_guess）から始まり、Scipyのfsolveが方程式に対して最適な値を見つけるために反復的な手法を使用します。
このプロセスは非線形方程式に対して一般的に使用され、数値的に解を見つけることができます。

この例では、与えられた非線形方程式を解くためにScipyを利用し、得られた数値解を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

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. 濃淡と電位の関係：

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

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