天気に関する方程式

天気に関する方程式の例として、気温の変化を表す単純なモデルを考えてみましょう。

例えば、一日の時間による気温の変化をモデル化することができます。

$$
T(t) = A \sin(\omega t + \phi) + C
$$

ここで、$ ( T(t) ) $は時間$ ( t ) $における気温を表し、$ ( A ) $は振幅、$ ( \omega ) $は角周波数、$ ( \phi ) $は位相差、$ ( C ) $は平均気温を表します。

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

# パラメータの設定
A = 10  # 振幅 (℃)
omega = 2 * np.pi / 24  # 角周波数 (1時間あたりのラジアン)
phi = 0  # 位相差
C = 20  # 平均気温 (℃)

# 時間の範囲を設定
t_values = np.linspace(0, 24, 100)  # 0時から24時までの範囲を100分割

# 気温の計算
T_values = A * np.sin(omega * t_values + phi) + C

# グラフの描画
plt.plot(t_values, T_values)
plt.title('Daily Temperature Variation')
plt.xlabel('Time (hours)')
plt.ylabel('Temperature (°C)')
plt.grid(True)
plt.show()

このコードでは、一日の時間に応じた気温の変化を示すグラフを描画しています。

[実行結果]

振幅、角周波数、位相差、平均気温などのパラメータを調整することで、異なる天候パターンに対応するグラフを作成できます。

ソースコード解説

このソースコードは、一日の時間に対する気温の変化をモデル化し、それをグラフ化するものです。

以下は、コードの詳細な説明です。

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

1
2

import numpy as np
import matplotlib.pyplot as plt

• numpy パッケージは数値計算を行うためのパッケージであり、配列や数学的な関数を提供します。
• matplotlib.pyplot モジュールは、データの可視化を行うためのライブラリであり、グラフを描画する際に使用されます。

2. パラメータの設定:

1
2
3
4

A = 10  # 振幅 (℃)
omega = 2 * np.pi / 24  # 角周波数 (1時間あたりのラジアン)
phi = 0  # 位相差
C = 20  # 平均気温 (℃)

• A は振幅であり、気温の変化の大きさを表します。
• omega は角周波数であり、1時間あたりのラジアンで表されます。
  ここでは、一日の周期を考慮しています。
• phi は位相差であり、波形の位相を調整します。
  ここでは、$0$としています。
• C は平均気温であり、波形の水平方向のシフトを表します。

3. 時間の範囲の設定:

1

t_values = np.linspace(0, 24, 100)  # 0時から24時までの範囲を100分割

• np.linspace() 関数を使用して、$0$から$24$までの範囲を$100$分割した時間の配列を作成します。

4. 気温の計算:

1

T_values = A * np.sin(omega * t_values + phi) + C

• 時間に対する気温の変化を計算します。
  ここでは、振幅と角周波数を用いた正弦関数を使っています。
  また、位相差と平均気温を加えています。

5. グラフの描画:

1
2
3
4
5
6

plt.plot(t_values, T_values)
plt.title('Daily Temperature Variation')
plt.xlabel('Time (hours)')
plt.ylabel('Temperature (°C)')
plt.grid(True)
plt.show()

• plt.plot() 関数を使用して、時間に対する気温の変化を描画します。
• plt.title()、plt.xlabel()、plt.ylabel() 関数を使用して、タイトルや軸ラベルを設定します。
• plt.grid(True) は、グリッドを表示するための命令です。
• 最後の plt.show() は、グラフを表示します。

結果解説

[実行結果]

このグラフは、一日の時間に対する気温の変化を表しています。

以下は、グラフの各要素の詳細です：

X軸 (Time)：

一日の時間を表します。
単位は時間です。
この軸は、24時間を通じて時間の経過に対する変化を示します。

Y軸 (Temperature)：

気温を表します。
単位は摂氏 (℃) です。この軸は、気温の変化を示します。

グラフの線：

$X$軸の値 (時間) に対する、$Y$軸の値 (気温) を表す曲線です。
この曲線は、一日の時間に応じて気温がどのように変化するかを示しています。
振幅、角周波数、位相差、平均気温などのパラメータによって、この曲線の形状が変化します。

グリッド：

背景に表示される格子状の線で、データの視覚的な評価を支援します。
水平方向と垂直方向の線が交わる点で、それぞれの時間と気温の値が示されます。

このグラフからは、一日の時間帯によって気温が変化する様子がわかります。

クーロンの法則

クーロンの法則は、電荷間の相互作用を記述する法則であり、以下の式で表されます：

$$
F = \frac{k \cdot q_1 \cdot q_2}{r^2}
$$

ここで、$ ( F ) $は電荷間の力、$ ( k ) $はクーロン定数 $ (( 8.9875 \times 10^9 , \text{N m}^2/\text{C}^2 )）$、$ ( q_1 ) $と$ ( q_2 ) $はそれぞれの電荷、$ ( r ) $は電荷間の距離です。

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

import numpy as np
import matplotlib.pyplot as plt

# クーロンの法則
def coulombs_law(q1, q2, r):
    k = 8.9875e9  # クーロン定数 (N m^2 / C^2)
    return k * q1 * q2 / r**2

# 電荷の設定
q1 = 1e-6  # 電荷1 (C)
q2 = -1e-6  # 電荷2 (C)

# 距離の範囲
r_values = np.linspace(1, 10, 100)  # 1 から 10 の範囲で距離を設定

# 電荷間の力を計算
force_values = coulombs_law(q1, q2, r_values)

# グラフ化
plt.plot(r_values, force_values)
plt.title("Coulomb's Law")
plt.xlabel('Distance (m)')
plt.ylabel('Force (N)')
plt.grid(True)
plt.show()

このコードは、クーロンの法則を用いて電荷間の力を計算し、結果をグラフ化します。

[実行結果]

ソースコード解説

以下は、ソースコードの詳細な説明です。

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

   1 2	import numpy as np import matplotlib.pyplot as plt

• numpy は数値計算を行うためのパッケージです。
  主に配列や行列などの数学的な演算をサポートします。
• matplotlib.pyplot はデータを可視化するためのライブラリです。
  グラフや図を描画する際に使用されます。

2. クーロンの法則を定義する関数:

   1 2 3

   def coulombs_law(q1, q2, r): k = 8.9875e9 # クーロン定数 (N m^2 / C^2) return k * q1 * q2 / r**2

• coulombs_law 関数は、クーロンの法則に基づいて電荷間の力を計算します。
• q1 と q2 はそれぞれの電荷の大きさを表し、r は電荷間の距離です。
• k はクーロン定数であり、$ (8.9875 \times 10^9 , \text{N m}^2/\text{C}^2) $で定義されています。

3. 電荷の設定:

   1 2	q1 = 1e-6 # 電荷1 (C) q2 = -1e-6 # 電荷2 (C)

• q1 と q2 は、それぞれの電荷の大きさを表す変数です。
  ここでは、$ 1e-6 $ C (クーロン) と$ -1e-6 $ Cの電荷を設定しています。

4. 距離の範囲の設定:

   1	r_values = np.linspace(1, 10, 100) # 1 から 10 の範囲で距離を設定

• np.linspace 関数を使用して、$ 1 $から$ 10 $の範囲で距離を等間隔で$100$点取得しています。

5. 電荷間の力の計算:

   1	force_values = coulombs_law(q1, q2, r_values)

• 先ほど定義した coulombs_law 関数を使用して、電荷間の距離に応じた力を計算しています。

6. グラフの作成と表示:

   1 2 3 4 5 6

   plt.plot(r_values, force_values) plt.title("Coulomb's Law") plt.xlabel('Distance (m)') plt.ylabel('Force (N)') plt.grid(True) plt.show()

• plt.plot 関数を使用して、距離に対する力のグラフを描画します。
• plt.title、plt.xlabel、plt.ylabel 関数を使用して、タイトルや軸ラベルを設定します。
• plt.grid(True) は、グリッドを表示するための命令です。
• 最後の plt.show() は、グラフを表示します。

結果解説

[実行結果]

このグラフは、クーロンの法則に基づいて計算された電荷間の力を距離に対してプロットしたものです。

以下はグラフの各要素の詳細です：

$X$軸 (Distance)：

電荷間の距離を表します。
単位はメートル (m) です。
この軸は、電荷同士の間隔がどのように変化するかを示しています。

$Y$軸 (Force)：

電荷間の力を表します。
単位はニュートン (N) です。
この軸は、電荷同士の間に生じる相互作用力がどのように変化するかを示しています。

グラフの線：

$X$軸の値 (距離) に対する、$Y$軸の値 (力) を表す曲線です。
クーロンの法則に従って、電荷間の距離が近づくにつれて力が増加し、距離が離れるにつれて力が減少することが期待されます。

グリッド：

背景に表示される格子状の線で、データの視覚的な評価を支援します。
水平方向と垂直方向の線が交わる点で、それぞれの距離と力の値が示されます。

このグラフからは、電荷間の力が距離の関数としてどのように振る舞うかがわかります。

距離が増加するにつれて、力が減少することが確認できます。

この振る舞いは、クーロンの法則の特徴的な振る舞いを示しています。

球面方程式

３次元方程式の例としては、例えば球面方程式が挙げられます。

球面方程式は以下のように表されます:

$$
x^2 + y^2 + z^2 = r^2
$$

ここで、$ ( r ) $は球の半径です。

以下に、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

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

# Sphere equation parameters
r = 1  # Radius

# Generate points on the sphere
theta = np.linspace(0, 2 * np.pi, 100)
phi = np.linspace(0, np.pi, 100)
x = r * np.outer(np.cos(theta), np.sin(phi))
y = r * np.outer(np.sin(theta), np.sin(phi))
z = r * np.outer(np.ones(100), np.cos(phi))

# Plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x, y, z, color='b', alpha=0.5)

# Set equal aspect ratio
ax.set_box_aspect([1,1,1])

# Set labels
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')

plt.title('Sphere')

plt.show()

このコードは、半径が$1$の球を描画しています。

[実行結果]

必要に応じて、r の値を変更して異なる半径の球を描画することができます。

ソースコード解説

このソースコードは、Pythonで半径が$1$の球面を描画するためのプログラムです。

以下はそれぞれの部分の詳細な説明です。

ライブラリのインポートとモジュールの選択

1
2
3

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

• numpyは、数値計算を行うためのPythonのライブラリです。
  主に多次元配列や行列演算、数学関数などを提供します。
• matplotlib.pyplotは、グラフ描画のためのライブラリです。
  主に2次元のグラフを描画するための機能を提供します。
• mpl_toolkits.mplot3dは、3次元のグラフ描画のためのツールキットです。
  これには、3Dプロットを作成するための機能が含まれています。

球面方程式のパラメータの設定

1

r = 1  # Radius

• rは球の半径を表します。
  この場合、半径は$1$と設定されています。

球面上の点の生成

1
2
3
4
5

theta = np.linspace(0, 2 * np.pi, 100)
phi = np.linspace(0, np.pi, 100)
x = r * np.outer(np.cos(theta), np.sin(phi))
y = r * np.outer(np.sin(theta), np.sin(phi))
z = r * np.outer(np.ones(100), np.cos(phi))

• thetaとphiはそれぞれ、水平方向と垂直方向の角度を表します。
  np.linspaceは、指定された範囲内で等間隔の値を生成します。
• np.outerは、2つのベクトルの外積を計算します。
  ここでは、球面上の点の座標を計算するために使用されます。
• x、y、zは、それぞれ球面上の点の$x$座標、$y$座標、$z$座標を表します。

グラフの描画

1
2
3

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x, y, z, color='b', alpha=0.5)

• plt.figure()は、新しい図を作成します。
• fig.add_subplot(111, projection='3d')は、3Dプロット用のサブプロットを作成します。
• ax.plot_surface(x, y, z, color='b', alpha=0.5)は、球面を描画します。
  x、y、zは、それぞれ球面上の点の座標を表し、colorは色を指定します。
  alphaは透明度を設定します。

グラフの設定

1
2
3
4
5

ax.set_box_aspect([1,1,1])
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.title('Sphere')

• ax.set_box_aspect([1,1,1])は、3Dプロットの表示を正確なアスペクト比で設定します。
• ax.set_xlabel('X')、ax.set_ylabel('Y')、ax.set_zlabel('Z')は、それぞれ$X$軸、$Y$軸、$Z$軸のラベルを設定します。
• plt.title('Sphere')は、グラフのタイトルを設定します。

グラフの表示

1

plt.show()

• plt.show()は、プロットを表示します。

これらの手順が組み合わさって、半径が$1$の球面が3Dグラフ上に描画されます。

結果解説

[実行結果]

以下はグラフに表示される内容の詳細です:

1. 3次元座標系:

グラフ上には、$X$軸、$Y$軸、$Z$軸があります。
これらはそれぞれ球の中心を通る直線です。

2. 球面:

グラフ上には、半径が$1$の球面が描かれます。
球面は、球の中心を中心として座標系の全方向に等しく広がっています。
球面は青色で描かれ、半透明になっています。

3. 等距離線:

球面上の等距離線は、球の中心から球面上の各点までの距離が等しい点の集合を示しています。
これは球面が円状の断面を持つことを視覚的に示しています。

4. プロットの表示:

3Dグラフでは、各軸上の値に対応する点がプロットされます。
この例では、球面上の多くの点がプロットされていますが、それらの点が球面上にあることが示されています。

5. タイトルと軸ラベル:

グラフには、「Sphere」というタイトルが表示され、各軸には「X」、「Y」、「Z」というラベルが付けられています。
これらのラベルは、それぞれ$X$軸、$Y$軸、$Z$軸を表します。

これらの要素が組み合わさって、3次元空間内で半径が$1$の球面が視覚化されます。

プラズマのエネルギー収支方程式

核融合に関する基本的な方程式として、例えばプラズマのエネルギー収支方程式があります。

これはプラズマ内でのエネルギーの生成と損失を表します。

ここでは、プラズマのエネルギー収支方程式を簡単なモデルで表現し、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 solve_ivp
import matplotlib.pyplot as plt

# プラズマのエネルギー収支方程式の定義
def plasma_energy_balance(t, E):
    # パラメータ
    alpha = 1.0  # エネルギー生成率
    beta = 0.5   # エネルギー損失率

    # エネルギーの時間変化
    dEdt = alpha - beta * E

    return dEdt

# 初期エネルギー
E0 = 0.0

# 時間範囲
t_span = (0, 10)  # 0から10までの時間

# 時間点
t_eval = np.linspace(0, 10, 100)

# 常微分方程式の解く
sol = solve_ivp(plasma_energy_balance, t_span, [E0], t_eval=t_eval)

# 結果のグラフ化
plt.figure(figsize=(10, 6))
plt.plot(sol.t, sol.y[0], label='Plasma Energy')
plt.xlabel('Time')
plt.ylabel('Energy')
plt.title('Plasma Energy Balance')
plt.legend()
plt.grid(True)
plt.show()

このコードでは、プラズマのエネルギー収支方程式を定義し、SciPyのsolve_ivp関数を使用してこの微分方程式を解いています。

そして、プラズマのエネルギーが時間とともにどのように変化するかをグラフ化しています。

[実行結果]

ソースコード解説

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

ライブラリのインポート

1
2
3

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

• numpy：数値計算を行うための基本的なライブラリ。
• scipy.integrate.solve_ivp：常微分方程式を解くための関数を提供するSciPyのモジュール。
• matplotlib.pyplot：データの可視化を行うためのグラフ描画ライブラリ。

プラズマのエネルギー収支方程式の定義

1
2
3
4
5
6
7
8
9

def plasma_energy_balance(t, E):
    # パラメータ
    alpha = 1.0  # エネルギー生成率
    beta = 0.5   # エネルギー損失率

    # エネルギーの時間変化
    dEdt = alpha - beta * E

    return dEdt

• plasma_energy_balance 関数：プラズマのエネルギー収支方程式を定義します。
  この関数は、時間 t と現在のエネルギー E を受け取り、時間変化率 dEdt を計算して返します。

初期エネルギーの設定

1

E0 = 0.0

• E0：初期のプラズマのエネルギーを設定します。

時間範囲の設定

1

t_span = (0, 10)  # 0から10までの時間

• t_span：解析する時間範囲を設定します。
  ここでは、$0$から$10$までの時間を指定しています。

時間点の設定

1

t_eval = np.linspace(0, 10, 100)

• t_eval：グラフ化する際の時間点を設定します。
  ここでは、$0$から$10$までの時間を$100$等分した時間点を指定しています。

常微分方程式の解く

1

sol = solve_ivp(plasma_energy_balance, t_span, [E0], t_eval=t_eval)

• solve_ivp 関数：プラズマのエネルギー収支方程式を解きます。
  初期エネルギー E0 から始まり、t_span の時間範囲で解析し、t_eval で指定された時間点で解を計算します。

結果のグラフ化

1
2
3
4
5
6
7
8

plt.figure(figsize=(10, 6))
plt.plot(sol.t, sol.y[0], label='Plasma Energy')
plt.xlabel('Time')
plt.ylabel('Energy')
plt.title('Plasma Energy Balance')
plt.legend()
plt.grid(True)
plt.show()

• プロットの設定：グラフのサイズ、ラベル、タイトル、凡例などを設定します。
• plt.plot：解析結果をグラフにプロットします。
  横軸に時間、縦軸にプラズマのエネルギーを設定し、ラベルを付けています。
• plt.show：グラフを表示します。

これにより、プラズマのエネルギー収支方程式を解いて時間変化をグラフ化する完全なプロセスが実行されます。

結果解説

[実行結果]

このグラフは、時間に対するプラズマのエネルギーの変化を示しています。

• 横軸は時間を表し、縦軸はプラズマのエネルギーを表します。
• グラフ上の曲線は、時間とともにプラズマのエネルギーがどのように変化するかを表します。
• 最初の時点では、初期エネルギーが$0$であるため、グラフの始点は原点になります。
• 曲線の傾きは、プラズマのエネルギーの時間変化率を示しており、この例では定数の生成率と損失率に基づいています。
• 時間が経過するにつれて、エネルギー生成と損失の影響により、プラズマのエネルギーが変化します。
• グラフ全体の形状は、生成率と損失率によって異なります。
  生成率が損失率よりも大きい場合、エネルギーは増加し続ける傾向がありますが、逆の場合はエネルギーが減少します。
• グラフが時間の経過とともにどのように振る舞うかは、定義したエネルギー収支方程式に基づいています。

このようなグラフを通じて、時間に対するプラズマのエネルギーの挙動を直感的に理解することができます。

流体力学 ナビエ・ストークス方程式

流体力学に関連する方程式として、ナビエ・ストークス方程式が代表的です。

これは流体の運動を記述する偏微分方程式であり、以下のように表されます：

1. 質量保存の式（連続の式）:

$$
\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0
$$

2. 運動方程式（運動量保存の式）:

$$
\rho \left( \frac{\partial \mathbf{v}}{\partial t} + \mathbf{v} \cdot \nabla \mathbf{v} \right) = - \nabla p + \mu \nabla^2 \mathbf{v} + \mathbf{f}
$$

ここで、$ (\rho) $は流体の密度、$ (\mathbf{v}) $は流速ベクトル、$ (p) $は圧力、$ (\mu) $は動粘性係数、$ (\mathbf{f}) $は体積力ベクトルです。

この方程式をPythonで解き、結果をグラフ化するには、適切な境界条件と数値解法を選択する必要があります。

一般的には有限体積法や有限要素法などの手法が用いられますが、ここでは簡単な場合として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

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
L = 1.0  # ドメインの長さ
Nx = 100  # 空間格子点数
Nt = 1000  # 時間ステップ数
dx = L / Nx  # 空間刻み
dt = 0.001  # 時間刻み
nu = 0.01  # 動粘性係数
rho = 1.0  # 流体の密度
force = 1.0  # 外力

# 初期条件設定
u = np.zeros(Nx)  # 流速
u[int(0.5 / dx):int(0.6 / dx)] = 1.0  # 初期条件: 一部の領域に初速度を与える

# シミュレーション
for t in range(Nt):
    un = np.copy(u)
    for i in range(1, Nx-1):
        u[i] = un[i] - dt / (2 * dx) * (un[i+1]**2 - un[i-1]**2) + \
               nu * dt / dx**2 * (un[i+1] - 2*un[i] + un[i-1]) + \
               dt / rho * force

# 結果のプロット
x = np.linspace(0, L, Nx)
plt.plot(x, u)
plt.xlabel('Position')
plt.ylabel('Velocity')
plt.title('1D Navier-Stokes Simulation')
plt.show()

このコードは、1次元の流れにおけるナビエ・ストークス方程式を解き、結果をプロットします。

初期条件として、ドメインの中央部分にのみ初速度が与えられます。

流速の時間変化をプロットすることで、流れの挙動を可視化することができます。

[実行結果]

ソースコード解説

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

ライブラリのインポート

1
2

import numpy as np
import matplotlib.pyplot as plt

• NumPy: 数値計算を行うためのライブラリ。
• Matplotlib: グラフの描画に使用されるライブラリ。

パラメータ設定

1
2
3
4
5
6
7
8

L = 1.0  # ドメインの長さ
Nx = 100  # 空間格子点数
Nt = 1000  # 時間ステップ数
dx = L / Nx  # 空間刻み
dt = 0.001  # 時間刻み
nu = 0.01  # 動粘性係数
rho = 1.0  # 流体の密度
force = 1.0  # 外力

• L: シミュレーションする空間の長さ。
• Nx: 空間を離散化する格子点の数。
• Nt: 時間を離散化するステップ数。
• dx: 空間刻み。
• dt: 時間刻み。
• nu: 流体の動粘性係数。
• rho: 流体の密度。
• force: 外力の強さ。

初期条件設定

1
2

u = np.zeros(Nx)  # 流速
u[int(0.5 / dx):int(0.6 / dx)] = 1.0  # 初期条件: 一部の領域に初速度を与える

• u: 流速を表す配列。
  最初は全ての要素がゼロに初期化されています。
• u[int(0.5 / dx):int(0.6 / dx)] = 1.0: 流速の一部の領域に初速度を与えています。
  この範囲はドメインの中心付近です。

シミュレーション

1
2
3
4
5
6

for t in range(Nt):
    un = np.copy(u)
    for i in range(1, Nx-1):
        u[i] = un[i] - dt / (2 * dx) * (un[i+1]**2 - un[i-1]**2) + \
               nu * dt / dx**2 * (un[i+1] - 2*un[i] + un[i-1]) + \
               dt / rho * force

• ループを使用して、時間ステップごとに流体の状態を更新します。
• 2重のループを使用して、空間内の各格子点での更新を行います。
• ナビエ・ストークス方程式の離散化バージョンを用いて、流速を更新しています。

結果のプロット

1
2
3
4
5
6

x = np.linspace(0, L, Nx)
plt.plot(x, u)
plt.xlabel('Position')
plt.ylabel('Velocity')
plt.title('1D Navier-Stokes Simulation')
plt.show()

• シミュレーションの結果をグラフ化します。
• x: ドメイン内の位置の配列。
• plt.plot(x, u): 位置に対する流速をプロットします。
• plt.xlabel('Position'): x軸のラベルを設定します。
• plt.ylabel('Velocity'): y軸のラベルを設定します。
• plt.title('1D Navier-Stokes Simulation'): グラフのタイトルを設定します。
• plt.show(): グラフを表示します。

結果解説

[実行結果]

このグラフは、1次元のナビエ・ストークス方程式を解いた結果であり、横軸には空間の位置、縦軸には流速を示しています。

具体的には、以下の内容がグラフに表示されます：

1. 位置（横軸）:

ドメインの長さに対する位置が横軸に表示されます。
ドメインの左端から右端までを表し、単位は任意の長さの単位（例えばメートル）になります。

2. 流速（縦軸）:

流体の速度が縦軸に表示されます。
流速は時間とともに変化し、位置によっても異なります。
流速の単位は長さの単位に時間の逆数を掛けたものになります（例えば m/s）。

3. グラフの形状:

初期条件として一部の領域に初速度が与えられているため、その部分の流速が一定であることが分かります。
その他の部分では、時間とともに流速が変化していることが示されています。
流体の流れがどのように振る舞うかは、初期条件や境界条件、流体の性質（動粘性係数や密度）に依存します。

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

グラフには、図の内容を説明するためのタイトルと軸ラベルが付いています。
タイトルには「1D Navier-Stokes Simulation」とあり、1次元のナビエ・ストークス方程式のシミュレーション結果であることが示されています。

$x軸$には「Position」（位置）が、$y軸$には「Velocity」（流速）が表示されています。

スキャッター方程式（Scatter equation）は、放射輸送方程式の一種であり、放射伝達を記述する非常に重要な方程式です。

スキャッター方程式は、放射粒子（光子など）が物質と相互作用する際の挙動をモデル化します。

放射輸送方程式の一般形式は以下の通りです：

$$
[
\frac{\partial I(\mathbf{r}, \mathbf{\Omega}, \lambda, t)}{\partial s} = -\sigma(\mathbf{r}, \lambda) I(\mathbf{r}, \mathbf{\Omega}, \lambda, t) + \sigma_s(\mathbf{r}, \lambda) \int_{4\pi} P(\mathbf{\Omega}’ \to \mathbf{\Omega}, \lambda) I(\mathbf{r}, \mathbf{\Omega}’, \lambda, t) d\mathbf{\Omega}’ + Q(\mathbf{r}, \mathbf{\Omega}, \lambda, t)
]
$$

ここで、$ (I(\mathbf{r}, \mathbf{\Omega}, \lambda, t)) $は空間位置$ (\mathbf{r}) $、方向$ (\mathbf{\Omega}) $、波長$ (\lambda) $、時間$ (t) $における放射束密度、$ (s) $は放射束が進む距離、$ (\sigma) $は吸収係数、$ (\sigma_s) $は散乱係数、$ (P(\mathbf{\Omega}’ \to \mathbf{\Omega} $,$ \lambda)) $は方向$ (\mathbf{\Omega}’) $から方向$ (\mathbf{\Omega}) $に散乱される確率密度関数、$ (Q) $は放射源項です。

これを適切な境界条件とともに解くことで、放射場の挙動を理解することができます。

以下に、スキャッター方程式を解くためのPythonコードを示します。

このコードでは、単純な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

import numpy as np
import matplotlib.pyplot as plt

# 定数とパラメータ
L = 1.0  # 空間の長さ
N = 100  # 空間の分割数
sigma_s = 0.1  # 散乱係数
sigma_a = 0.2  # 吸収係数
Q = 1.0  # 放射源項
dx = L / N  # 空間の分割幅

# 初期条件
I0 = np.zeros(N)  # 放射束密度
I0[0] = 1.0  # 入射する放射束

# スキャッター方程式を解く関数
def solve_scatter_eq(I0, sigma_s, sigma_a, Q, dx):
    I = np.copy(I0)
    for i in range(N):
        I[i] = (Q + sigma_s * sum(I0[:i]) * dx) / (sigma_a + sigma_s)
    return I

# スキャッター方程式を解く
I = solve_scatter_eq(I0, sigma_s, sigma_a, Q, dx)

# 結果のプロット
x_values = np.linspace(0, L, N)
plt.plot(x_values, I)
plt.title('放射束密度の分布')
plt.xlabel('空間の位置')
plt.ylabel('放射束密度')
plt.grid(True)
plt.show()

このコードでは、1次元の空間内でスキャッター方程式を解き、放射束密度の分布をプロットしています。

[実行結果]

ソースコード解説

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

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

1
2

import numpy as np
import matplotlib.pyplot as plt

• numpy：数値計算を行うためのライブラリ。配列処理などをサポートします。
• matplotlib.pyplot：データの可視化（グラフ化）のためのライブラリ。

2. 定数とパラメータの設定

1
2
3
4
5
6

L = 1.0  # 空間の長さ
N = 100  # 空間の分割数
sigma_s = 0.1  # 散乱係数
sigma_a = 0.2  # 吸収係数
Q = 1.0  # 放射源項
dx = L / N  # 空間の分割幅

• L：空間の長さ。
• N：空間をいくつの分割に区切るかを示す値。
• sigma_s：散乱係数。
• sigma_a：吸収係数。
• Q：放射源項。
• dx：空間の分割幅。

3. 初期条件の設定

1
2

I0 = np.zeros(N)  # 放射束密度
I0[0] = 1.0  # 入射する放射束

• I0：初期の放射束密度。
  全ての要素がゼロで初期化されます。
• I0[0] = 1.0：入射する放射束。
  空間の最初のセルにのみ入射する放射束が1.0と設定されます。

4. スキャッター方程式を解く関数の定義

1
2
3
4
5

def solve_scatter_eq(I0, sigma_s, sigma_a, Q, dx):
    I = np.copy(I0)
    for i in range(N):
        I[i] = (Q + sigma_s * sum(I0[:i]) * dx) / (sigma_a + sigma_s)
    return I

• solve_scatter_eq：与えられた初期条件とパラメータに基づいて、スキャッター方程式を解く関数。
• 方程式はループを使って数値的に解かれます。
• 解析解ではなく、数値解法を使用しています。

5. スキャッター方程式を解く

1

I = solve_scatter_eq(I0, sigma_s, sigma_a, Q, dx)

• solve_scatter_eq 関数を使ってスキャッター方程式を解き、放射束密度の配列 I に結果を格納します。

6. 結果のプロット

1
2
3
4
5
6
7

x_values = np.linspace(0, L, N)
plt.plot(x_values, I)
plt.title('放射束密度の分布')
plt.xlabel('空間の位置')
plt.ylabel('放射束密度')
plt.grid(True)
plt.show()

• plt.plot(x_values, I)：放射束密度の分布を折れ線グラフとしてプロットします。
• plt.title('放射束密度の分布')：グラフのタイトルを設定します。
• plt.xlabel('空間の位置')：$x軸$のラベルを設定します。
• plt.ylabel('放射束密度')：$y軸$のラベルを設定します。
• plt.grid(True)：グリッドを表示します。
• plt.show()：グラフを表示します。

これらの手順によって、スキャッター方程式の解析結果が放射束密度の分布として可視化されます。

結果解説

下記のグラフは、放射束密度の空間内での分布を表しています。

[実行結果]

以下は、グラフに表示される内容の詳細な説明です：

• 横軸（x軸）は空間の位置を表します。
  1次元の空間が均等に分割されており、それぞれの分割（セル）の中心が表示されています。
• 縦軸（y軸）は放射束密度を表します。
  放射束密度は各位置における放射の強度を示し、単位長さあたりの放射エネルギーを表します。
• グラフは、放射源から放射束が入射し、空間内を進んでいく様子を示しています。
  放射源からの入射は、グラフの左端にある位置から始まります。
• 入射した放射束が空間内を進むにつれて、吸収や散乱によって減衰したり増加したりします。
  これにより、放射束密度の分布が変化します。
• 散乱係数 $（(\sigma_s)）$と吸収係数 $（(\sigma_a)）$の値によって、放射の挙動が変化します。
  例えば、散乱が支配的な場合は放射が拡散し、吸収が支配的な場合は放射が減衰します。
• 放射源項 $（(Q)）$は、空間内に追加の放射を与える項です。
  放射源からの入射だけでなく、空間内での放射生成を考慮しています。

このグラフは、スキャッター方程式に基づいて、放射輸送の基本的な挙動を可視化しています。

有名な数学的方程式の１つとして、シュレーディンガー方程式があります。

これは量子力学における基本的な方程式であり、粒子の波動関数の時間変化を記述します。

一次元の自由粒子のシュレーディンガー方程式は以下のように表されます：

$$
i\hbar \frac{\partial}{\partial t} \Psi(x, t) = -\frac{\hbar^2}{2m} \frac{\partial^2}{\partial x^2} \Psi(x, t)
$$

ここで、$ ( \Psi(x, t) ) $は波動関数、$ ( x ) $は空間座標、$ ( t ) $は時間、$ ( \hbar ) $はディラック定数、$ ( m ) $は粒子の質量です。

この方程式を解いて、時間と空間に依存する波動関数をプロットします。

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

import numpy as np
import matplotlib.pyplot as plt

# 定数の設定
h_bar = 1.0  # ディラック定数
m = 1.0      # 粒子の質量

# 空間および時間の設定
x = np.linspace(-5, 5, 400)
t = np.linspace(0, 2*np.pi, 400)

# 波動関数の計算
X, T = np.meshgrid(x, t)
Psi = np.exp(-(X**2)/2) * np.exp(1j * T)

# 結果のプロット
plt.figure(figsize=(10, 6))
plt.imshow(np.abs(Psi)**2, extent=(-5, 5, 0, 2*np.pi), aspect='auto', cmap='hot')
plt.colorbar(label='Probability Density')
plt.title('Probability Density of a Free Particle')
plt.xlabel('Position')
plt.ylabel('Time')
plt.show()

このコードは、一次元の自由粒子のシュレーディンガー方程式を解き、時間と空間に依存する波動関数の確率密度をプロットします。

確率密度は波動関数の絶対値の二乗として与えられ、波動関数が存在する確率を表します。

[実行結果]

ソースコード解説

以下にソースコードを詳しく説明します。

1. import numpy as np:

• NumPy ライブラリを np としてインポートします。
  NumPy は、数値計算を行うための強力なライブラリであり、配列操作や数学関数の計算に便利です。

2. import matplotlib.pyplot as plt:

• Matplotlib ライブラリの pyplot モジュールを plt としてインポートします。
  Matplotlib は、グラフを描画するためのライブラリであり、pyplot モジュールはその中でも一般的に使用されます。

3. h_bar = 1.0 および m = 1.0:

• シュレーディンガー方程式に登場する定数$ ( \hbar )$（ディラック定数）と粒子の質量$ ( m ) $の値を設定します。
  ここでは両方とも値を$ 1.0 $に設定しています。

4. x = np.linspace(-5, 5, 400) および t = np.linspace(0, 2*np.pi, 400):

• 空間座標$ ( x ) $および時間$ ( t ) $の範囲をそれぞれ -5 から 5 までの区間と 0 から 2π までの区間に分割し、それぞれ$ 400 $個の等間隔な点で表現します。
  これにより、グラフ上で表示される座標の範囲が決まります。

5. X, T = np.meshgrid(x, t):

• NumPy の meshgrid() 関数を使って、空間座標 x と時間 t を格子状の座標系に変換します。
  これにより、各点での波動関数の値を計算するための格子が作成されます。

6. Psi = np.exp(-(X**2)/2) * np.exp(1j * T):

• シュレーディンガー方程式の解である波動関数$ ( \Psi(x, t) ) $を計算します。
  ここで、np.exp() は指数関数を計算し、1j は虚数単位を表します。
  この波動関数は、時間と空間に依存する複素数で表現されます。

7. plt.figure(figsize=(10, 6)):

• Matplotlib で図を作成します。
  図のサイズを幅$10$インチ、高さ$6$インチに設定しています。

8. plt.imshow(np.abs(Psi)**2, extent=(-5, 5, 0, 2*np.pi), aspect='auto', cmap='hot'):

• imshow() 関数を使って、波動関数の絶対値の二乗である確率密度をイメージとして表示します。
  extent 引数は$ x $軸と$ y $軸の範囲を指定し、aspect 引数はアスペクト比を自動調整するように指定します。
  cmap 引数はカラーマップを指定します。

9. plt.colorbar(label='Probability Density'):

• カラーバーをグラフに追加し、確率密度の色に対応する値を表示します。

10. plt.title('Probability Density of a Free Particle'):

• グラフのタイトルを設定します。

11. plt.xlabel('Position') および plt.ylabel('Time'):

• $ x $軸と$ y $軸のラベルを設定します。

12. plt.show():

• グラフを表示します。

結果解説

下記のグラフは、一次元自由粒子のシュレーディンガー方程式の解に基づいて、時間と空間に依存する波動関数の確率密度を表しています。

[実行結果]

横軸は空間座標$ ( x ) $を表し、縦軸は時間$ ( t ) $を表しています。

グラフの各点における色の濃さは、その空間座標と時間における波動関数の確率密度を示しています。

色が濃い領域ほど確率密度が高く、粒子が存在する可能性が高いことを意味します。

このグラフは、波動関数の時間変化を視覚的に理解するのに役立ちます。

時間が経過するにつれて、波動関数の形状が変化し、特定の空間座標で粒子の存在確率が変動します。
特に、波動関数の振動が時間によってどのように進行し、確率密度がどのように変化するかが示されています。

宇宙の膨張 フリードマン方程式

宇宙の膨張を表すフリードマン方程式を考えます。

フリードマン方程式は、宇宙のスケールファクター $ ( a(t) ) $とエネルギー密度 $ ( \rho(t) ) $に関連します。

方程式は以下の通りです：

$$
H^2 = \left(\frac{\dot{a}}{a}\right)^2 = \frac{8\pi G}{3}\rho - \frac{k}{a^2}
$$

ここで、$ ( H ) $はハッブル定数、$ ( G ) $は重力定数、$ ( k ) $は空間の曲率を表します。

また、ドットは時間に関する微分を表しています。

この方程式を解くには、初期条件やエネルギー密度などが必要です。

一般的には、標準的な宇宙論の初期条件を使用します。

以下は、Pythonを使用してこの方程式を解き、グラフ化する簡単な例です。

この例では、ハッブル定数 $ ( H_0 ) $を$1$、重力定数 $ ( G ) $を$1$、曲率 $ ( k ) $を$0$と仮定します。

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

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

# フリードマン方程式
def friedmann_equation(y, t):
    H0 = 1.0  # ハッブル定数
    G = 1.0   # 重力定数
    k = 0.0   # 曲率
    rho = y[0]  # エネルギー密度

    # フリードマン方程式
    dydt = [y[1], (8 * np.pi * G / 3) * rho - k / y[0]**2]
    return dydt

# 初期条件
initial_condition = [1.0, 0.0]  # 初期エネルギー密度、初期速度

# 時間の範囲
t = np.linspace(0, 10, 1000)

# 方程式を解く
solution = odeint(friedmann_equation, initial_condition, t)

# 結果をプロット
plt.plot(t, solution[:, 0], label='Scale Factor (a(t))')
plt.xlabel('Time')
plt.ylabel('Scale Factor')
plt.title('Friedmann Equation Solution')
plt.legend()
plt.show()

この例では、初期エネルギー密度を$1$、初期速度を$0$としています。

得られた結果は、時間に対する宇宙のスケールファクターの変化を示しています。

これは、宇宙が膨張していく様子を表しています。

[実行結果]

ソースコード解説

このソースコードは、Pythonを使用して宇宙の膨張を表すフリードマン方程式を解き、その解をプロットするためのものです。

以下はソースコードの章立てと詳細な説明です：

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

1
2
3

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

• numpy：数値計算用のライブラリ。
  数学的な演算や配列操作に使用されます。
• matplotlib.pyplot：データの可視化のためのライブラリ。
  グラフの作成や表示に利用されます。
• scipy.integrate.odeint：常微分方程式を解くための関数が含まれています。

2. フリードマン方程式の定義

1
2
3
4
5
6
7
8
9

def friedmann_equation(y, t):
    H0 = 1.0  # ハッブル定数
    G = 1.0   # 重力定数
    k = 0.0   # 曲率
    rho = y[0]  # エネルギー密度

    # フリードマン方程式
    dydt = [y[1], (8 * np.pi * G / 3) * rho - k / y[0]**2]
    return dydt

• friedmann_equation 関数：

フリードマン方程式を定義します。
この関数は、時間 t と状態ベクトル y を受け取り、その微分を返します。

• H0：ハッブル定数。
  宇宙の膨張を示す定数。
• G：重力定数。
  重力の強さを示す定数。
• k：曲率。
  宇宙の空間が平坦か、閉じているか、開いているかを示すパラメータ。
• rho：エネルギー密度。
  宇宙のエネルギー分布を表す。
• dydt：微分方程式の右辺。
  宇宙の膨張を示すフリードマン方程式の微分方程式。

3. 初期条件の設定

1

initial_condition = [1.0, 0.0]  # 初期エネルギー密度、初期速度

• initial_condition：フリードマン方程式の初期条件。
  初期エネルギー密度と初期速度が設定されています。

4. 時間の範囲の設定

1

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

• t：時間の範囲。$0$から$10$までの範囲を等間隔に区切った$1000$のデータポイントを生成します。

5. 方程式の解の計算

1

solution = odeint(friedmann_equation, initial_condition, t)

• odeint 関数：常微分方程式を解くための関数。
  friedmann_equation 関数、初期条件 initial_condition、時間 t を引数として渡し、方程式の解を計算します。

6. 結果のプロット

1
2
3
4
5
6

plt.plot(t, solution[:, 0], label='Scale Factor (a(t))')
plt.xlabel('Time')
plt.ylabel('Scale Factor')
plt.title('Friedmann Equation Solution')
plt.legend()
plt.show()

• plt.plot：解をプロットします。
  横軸に時間 t、縦軸にスケールファクター a(t) をプロットします。
• plt.xlabel：X軸のラベルを設定。
• plt.ylabel：Y軸のラベルを設定。
• plt.title：グラフのタイトルを設定。
• plt.legend：凡例を表示。
• plt.show：グラフを表示します。

このコードは、フリードマン方程式を解いて宇宙の膨張を表すグラフを生成します。

初期条件やパラメータの変更によって、グラフの形状が変わります。

結果解説

[実行結果]

上記のPythonコードで生成されるグラフは、フリードマン方程式を解いた結果、時間に対する宇宙のスケールファクター$ (a(t)) $の変化を表しています。

以下は、グラフに表示される内容の詳細な説明です：

1. 横軸（X軸） - 時間：

グラフの横軸は時間を表しています。
この例では、時間の範囲を$ (0) $から$ (10) $までとし、$ (1000) $の等間隔なデータポイントで構成されています。

2. 縦軸（Y軸） - スケールファクター $ (a(t)) $：

グラフの縦軸は宇宙のスケールファクター $ (a(t)) $を表しています。
スケールファクターは宇宙の膨張を示し、大きくなるほど宇宙が膨張していることを示します。

3. グラフの形状：

グラフは時間に対する宇宙の膨張を示しており、スケールファクターが増加していく様子を観察できます。
これは、ビッグバン理論に基づく宇宙の進化を示しています。
初期条件によってグラフの形状が変わります。

4. 初期条件：

この例では初期エネルギー密度を$ (1) $、初期速度を$ (0) $と仮定しています。
これらの初期条件に基づいて、時間が経過するにつれて宇宙がどのように膨張しているかが反映されています。

このグラフは、ビッグバンモデルにおける宇宙の膨張をシミュレートしたものであり、初期条件やパラメータの変更によって異なる結果が得られます。

常微分方程式

微分方程式の例として、単純な1階の常微分方程式を取り上げましょう。

以下の例では、微分方程式 $ ( \frac{dy}{dt} = -ky ) $を考えます。

この方程式は指数減衰を表します。

これを解いて、その結果をグラフ化します。

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

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

# 微分方程式の定義
def model(y, t, k):
    dydt = -k * y
    return dydt

# 初期条件
y0 = 5.0

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

# パラメータ
k = 0.1

# 微分方程式を解く
solution = odeint(model, y0, t, args=(k,))

# 結果のグラフ化
plt.figure(figsize=(8, 6))
plt.plot(t, solution, label='y(t)')
plt.title('指数減衰：$\\frac{dy}{dt} = -ky$')
plt.xlabel('時間 (t)')
plt.ylabel('y(t)')
plt.legend()
plt.grid(True)
plt.show()

このコードでは、SciPyのodeint関数を使用して微分方程式を解き、結果をグラフ化しています。

[実行結果]

指数減衰が確認できるでしょう。

パラメータ$ (k) $を変更することで解の挙動が変わります。

ソースコード解説

以下はソースコードの詳細な説明です。

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

1
2
3

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

• numpy: 数値計算用ライブラリで、数値配列や数学関数などを提供します。
• matplotlib.pyplot: グラフの描画に使用されるライブラリ。
• scipy.integrate.odeint: 常微分方程式を数値的に解くための関数が含まれているSciPyライブラリの一部。

2. 微分方程式の定義:

1
2
3

def model(y, t, k):
   dydt = -k * y
   return dydt

• model関数は、微分方程式$ ( \frac{dy}{dt} = -ky ) $の右辺を定義しています。
  この関数は odeint 関数によって呼び出されます。

3. 初期条件の設定:

1

y0 = 5.0

• 微分方程式の初期条件$ (y(0)) $を設定します。

4. 時間の範囲の設定:

1

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

• 時間の範囲を$0$から$10$までの区間を$100$点で均等に区切って定義します。

5. パラメータの設定:

1

k = 0.1

• 微分方程式のパラメータ$ (k) $を設定します。

6. 微分方程式を解く:

1

solution = odeint(model, y0, t, args=(k,))

• odeint 関数を使用して微分方程式を解きます。
  初期条件 y0、時間の範囲 t、およびパラメータ k を指定しています。

7. 結果のグラフ化:

1
2
3
4
5
6
7
8

plt.figure(figsize=(8, 6))
plt.plot(t, solution, label='y(t)')
plt.title('指数減衰：$\\frac{dy}{dt} = -ky$')
plt.xlabel('時間 (t)')
plt.ylabel('y(t)')
plt.legend()
plt.grid(True)
plt.show()

• matplotlib を使用して、微分方程式の解 solution を時間に対してプロットし、グラフを表示します。
• グラフのタイトル、軸ラベル、凡例、グリッドが設定されています。

このプログラムは、指数減衰を示す微分方程式を数値的に解き、その解をグラフで視覚化する基本的な例です。

結果解説

[実行結果]

以下はグラフに表示される内容の詳細な説明です。

横軸 (X軸):時間$ (t) $

$0$から$10$の範囲を$100$点で等分しています。

縦軸 (Y軸):変数$ (y(t))$

微分方程式の解であり、指数減衰の挙動を示します。
初期条件$ (y(0)) $は$ 5.0 $としています。

グラフの線:

微分方程式を解いた結果の$ (y(t)) $の変化を表しています。
指数減衰の特徴的な形状で、初めは急速に減少し、後に緩やかな減少に移ります。

グラフのタイトル:指数減衰：$ \frac{dy}{dt} = -ky $

微分方程式の性質と解いたものが指数減衰を示すことを示しています。

X軸およびY軸のラベル:

時間（時間$ (t)$）および変数$ (y(t)) $の説明があります。

凡例 (Legend):$y(t)$

グラフ内の線が何を示しているかを示すラベルです。

グリッド:

グラフ内に格子状の線が描かれ、視認性を向上させています。

このグラフは、微分方程式の解を通して、指数減衰が時間とともにどのように進行するかを示しています。

指数減衰は、物理学や生態学などのさまざまな分野で重要な概念であり、この例ではPythonを用いて簡単に表現されています。

最適化問題

最適化問題の例として、特定の制約条件下での関数の最小化を考えてみましょう。

以下は、2変数の関数の最小化問題です。

この例では、scipy.optimize ライブラリを使用して最適化を行います。

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

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

# 最小化する関数
def objective_function(x):
    return x[0]**2 + x[1]**2  # 例として x^2 + y^2 を最小化

# 制約条件
constraints = ({'type': 'eq', 'fun': lambda x: x[0] + x[1] - 5})  # 制約条件: x + y = 5

# 初期値
initial_guess = [0, 0]

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

# 結果表示
print("最適解:", result.x)
print("最小値:", result.fun)

# グラフの作成
x_vals = np.linspace(-2, 7, 100)
y_vals = np.linspace(-2, 7, 100)
X, Y = np.meshgrid(x_vals, y_vals)
Z = X**2 + Y**2

# 制約条件を描画
plt.plot(x_vals, 5 - x_vals, label='制約条件: x + y = 5', color='red')

# 最小値を描画
plt.scatter(result.x[0], result.x[1], color='green', marker='x', label='最小値')

# 関数の等高線を描画
contour = plt.contour(X, Y, Z, levels=15, cmap='viridis')
plt.clabel(contour, inline=True, fontsize=8)

# グラフの装飾
plt.title('2変数関数の最小化と制約条件')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True)

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

このコードでは、関数$ x^2 + y^2 $を最小化する問題を考え、制約条件として$ x + y = 5 $を設定しています。

最適化結果や最小化される関数のグラフが表示されます。

[実行結果]

異なる関数や制約条件を試すことで、様々な最適化問題に対応できます。

ソースコード解説

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

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. 最小化する関数の定義:

1
2

def objective_function(x):
    return x[0]**2 + x[1]**2  # 例として x^2 + y^2 を最小化

• x は2要素の配列で、x[0] と x[1] が変数となります。
• $ x[0]^2 + x[1]^2 $を最小化する関数を定義しています。

3. 制約条件の設定:

1

constraints = ({'type': 'eq', 'fun': lambda x: x[0] + x[1] - 5})  # 制約条件: x + y = 5

• 制約条件として、$ x + y = 5 $を設定しています。

4. 初期値の設定:

1

initial_guess = [0, 0]

• 最適化の初期値を [0, 0] に設定します。

5. 最適化の実行:

1

result = minimize(objective_function, initial_guess, constraints=constraints)

• minimize 関数を使用して最適化を実行します。
• objective_function を最小化し、制約条件 constraints の下で最適解を見つけます。

6. 結果の表示:

1
2

print("最適解:", result.x)
print("最小値:", result.fun)

• 最適解と最小値を表示します。

7. グラフの作成:

1
2
3
4

x_vals = np.linspace(-2, 7, 100)
y_vals = np.linspace(-2, 7, 100)
X, Y = np.meshgrid(x_vals, y_vals)
Z = X**2 + Y**2

• グラフ描画用の x, y 座標を生成します。
• Z には$ X^2 + Y^2 $の関数値が格納されます。

8. 制約条件の描画:

1

plt.plot(x_vals, 5 - x_vals, label='制約条件: x + y = 5', color='red')

• 制約条件を赤い線で描画します。

9. 最小値の描画:

1

plt.scatter(result.x[0], result.x[1], color='green', marker='x', label='最小値')

• 最小値の点を緑の “x” マークで描画します。

10. 等高線プロットの描画:

1
2

contour = plt.contour(X, Y, Z, levels=15, cmap='viridis')
plt.clabel(contour, inline=True, fontsize=8)

• 関数の等高線プロットを描画します。

11. グラフの装飾:

1
2
3
4
5

plt.title('2変数関数の最小化と制約条件')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True)

• グラフのタイトルや軸ラベル、凡例、グリッドを設定します。

12. グラフの表示:

1

plt.show()

• これまでの設定で作成したグラフを表示します。

このコードは、2変数の関数を最小化し、その結果をグラフに可視化しています。

特に、制約条件下での最小値を探索している点がポイントです。

結果解説

[実行結果]

この最適化問題は、関数$ ( f(x, y) = x^2 + y^2 ) $を最小化するというもので、制約条件として$ ( x + y = 5 ) $が課されています。

最適化アルゴリズムは、この関数を最小化する$ ( x ) と ( y ) $の値を見つけることを試みます。

結果は以下の通りです：

• 最適解: [2.5, 2.5]
• 最小値: 12.5

最適解 $ [2.5, 2.5] $は、制約条件 $ ( x + y = 5 ) $を満たす中で$ ( f(x, y) ) $を最小化する点です。

最小値 $ 12.5 $は、この最小化された関数の値です。

グラフには、等高線プロットが表示されています。

等高線は、関数の高さを示し、最小値がグリーンの “x” マークで示されています。

また、制約条件 $ ( x + y = 5 ) $を赤い線で示しています。

最適解がこの赤い線上にあり、かつ関数を最小化する点であることが確認できます。

この例は非常にシンプルですが、最適化問題の基本的な理解を提供しています。

実際の応用では、より複雑な関数や複数の制約条件がありますが、同様のアプローチが使われます。