複利計算

投資に関する方程式の例として、単純な複利計算を考えてみましょう。

複利計算の方程式は以下のように表されます:

$$
A = P \times (1 + r)^n
$$

• $ (A) $は最終的な金額（将来価値）
• $ (P) $は元本（現在価値）
• $ (r) $は年利率
• $ (n) $は投資期間（年数）

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

import numpy as np
import matplotlib.pyplot as plt

def compound_interest(P, r, n):
    """
    複利計算を行う関数
    
    Parameters:
        P (float): 元本（現在価値）
        r (float): 年利率
        n (int): 投資期間（年数）
    
    Returns:
        np.ndarray: 各年の金額の配列
    """
    years = np.arange(1, n+1)
    A = P * (1 + r)**years
    return A

# パラメータの設定
P = 1000  # 元本
r = 0.05  # 年利率 (5%)
n = 10    # 投資期間（10年）

# 複利計算を実行
investment_values = compound_interest(P, r, n)

# グラフの描画
plt.plot(np.arange(1, n+1), investment_values, marker='o')
plt.title('Compound Interest over {} Years'.format(n))
plt.xlabel('Years')
plt.ylabel('Investment Value')
plt.grid(True)
plt.show()

このコードは、元本$ (P) $、年利率$ (r) $、投資期間$ (n) $を与えられたときに、各年の投資金額$ (A) $を計算し、それをグラフ化しています。

[実行結果]

ソースコード解説

このソースコードは、Pythonを使用して複利計算を行い、その結果をグラフ化するプログラムです。

以下では、各部分の詳細を説明します。

ライブラリのインポート

1
2

import numpy as np
import matplotlib.pyplot as plt

• numpy は数値計算を行うためのライブラリです。
  ここでは主に配列操作や数学的な計算に使用されます。
• matplotlib.pyplot はグラフの描画に使用されるライブラリです。
  plt として一般的にエイリアスされます。

複利計算の関数定義

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

def compound_interest(P, r, n):
    """
    複利計算を行う関数
    
    Parameters:
        P (float): 元本（現在価値）
        r (float): 年利率
        n (int): 投資期間（年数）
    
    Returns:
        np.ndarray: 各年の金額の配列
    """
    years = np.arange(1, n+1)
    A = P * (1 + r)**years
    return A

• compound_interest 関数は、与えられた元本$ (P) $、年利率$ (r) $、投資期間$ (n) $に基づいて複利計算を行います。
  各年の金額の配列を返します。

パラメータの設定

1
2
3

P = 1000  # 元本
r = 0.05  # 年利率 (5%)
n = 10    # 投資期間（10年）

• P は元本（投資の始めの金額）を表します。
• r は年利率を表します。ここでは$ 5% $を表す$ 0.05 $が設定されています。
• n は投資期間（年数）を表します。ここでは$ 10 $年が設定されています。

複利計算の実行

1

investment_values = compound_interest(P, r, n)

• compound_interest 関数を使用して、指定されたパラメータで複利計算を実行し、各年の金額の配列を取得します。

グラフの描画

1
2
3
4
5
6

plt.plot(np.arange(1, n+1), investment_values, marker='o')
plt.title('Compound Interest over {} Years'.format(n))
plt.xlabel('Years')
plt.ylabel('Investment Value')
plt.grid(True)
plt.show()

• plt.plot() は、与えられたデータをプロットしてグラフを作成します。
  ここでは、各年の金額の配列を横軸に対してプロットしています。
• plt.title() はグラフのタイトルを設定します。
• plt.xlabel() と plt.ylabel() はそれぞれ横軸と縦軸のラベルを設定します。
• plt.grid(True) はグリッド線を表示します。
• plt.show() は、グラフを表示します。

これにより、元本$ (P) $を元にした複利計算の結果が$ 10 $年間の投資期間にわたってグラフ化されます。

結果解説

[実行結果]

このグラフは、投資の元本$ (P) $、年利率$ (r) $、投資期間$ (n) $に基づいて計算された複利効果を示しています。

横軸は時間を表し、投資期間が経過するにつれて増加します。
縦軸は投資の価値を表し、元本$ (P) $から始まり、時間が経つにつれて複利効果により増加していきます。

具体的には、投資期間の各年における投資の価値が点で示されています。
この点は、元本$ (P) $を年利率$ (r) $で複利計算した結果を示しています。
グラフ全体として、時間が経過するにつれて投資の価値が増加していく様子が視覚化されています。

グラフのタイトルには、投資期間が何年かが示されており、横軸には各年が、縦軸には投資の価値がラベル付けされています。
また、グリッド線が描かれており、各年ごとの価値を読み取りやすくしています。

このグラフは、投資が時間の経過とともにどのように成長するかを視覚的に理解するのに役立ちます。

動力学

航空学に関連する方程式の例としては、航空機の上昇および降下の動力学を記述するための運動方程式があります。

典型的な航空機の運動方程式は、重力、揚力、抗力、推力などの力を考慮しています。

以下に、簡単なモデルとして航空機の上昇と降下を表現する方程式を示します。

1. 運動方程式：$ ( F = m \cdot a ) $
2. 揚力（Lift）：$ ( L = \frac{1}{2} \cdot C_L \cdot \rho \cdot V^2 \cdot S ) $
3. 重力（Gravity）：$ ( W = m \cdot g ) $
4. 抗力（Drag）：$ ( D = \frac{1}{2} \cdot C_D \cdot \rho \cdot V^2 \cdot S ) $
5. 推力（Thrust）：$ ( T ) $

これらの方程式を組み合わせて、航空機の上昇および降下の挙動をモデル化します。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57

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

# 航空機の運動方程式を定義
def aircraft_motion_equation(state, t, thrust, drag_coefficient, lift_coefficient, mass, gravity):
    # 状態変数
    altitude, velocity = state
    
    # 密度の簡易モデル（高度による変化を無視）
    rho = 1.225
    
    # 揚力
    lift = 0.5 * lift_coefficient * rho * velocity**2
    
    # 重力
    weight = mass * gravity
    
    # 抗力
    drag = 0.5 * drag_coefficient * rho * velocity**2
    
    # 上昇力
    lift_force = lift - weight
    
    # 加速度
    acceleration = (thrust - drag) / mass
    
    return [velocity, acceleration]

# パラメータの設定
thrust = 50000  # 推力（N）
drag_coefficient = 0.3  # 抗力係数
lift_coefficient = 0.4  # 揚力係数
mass = 5000  # 航空機の質量（kg）
gravity = 9.81  # 重力加速度（m/s^2）

# 初期条件
initial_altitude = 0  # 初期高度（m）
initial_velocity = 0  # 初期速度（m/s）
initial_state = [initial_altitude, initial_velocity]

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

# 運動方程式を解く
result = odeint(aircraft_motion_equation, initial_state, t, args=(thrust, drag_coefficient, lift_coefficient, mass, gravity))
altitude, velocity = result.T

# グラフのプロット
plt.plot(t, altitude, label='Altitude')
plt.plot(t, velocity, label='Velocity')
plt.xlabel('Time (s)')
plt.ylabel('Altitude (m) / Velocity (m/s)')
plt.title('Aircraft Motion')
plt.legend()
plt.grid(True)
plt.show()

このコードは、航空機の運動方程式を解き、時間に対する高度と速度の変化をグラフ化しています。

航空機の上昇および降下の動力学を簡単にモデル化することができます。

[実行結果]

ソースコード解説

このソースコードは、航空機の運動方程式を解いて結果をグラフ化する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
10
11
12
13
14
15
16
17
18
19
20
21
22
23

def aircraft_motion_equation(state, t, thrust, drag_coefficient, lift_coefficient, mass, gravity):
    # 状態変数
    altitude, velocity = state
    
    # 密度の簡易モデル（高度による変化を無視）
    rho = 1.225
    
    # 揚力
    lift = 0.5 * lift_coefficient * rho * velocity**2
    
    # 重力
    weight = mass * gravity
    
    # 抗力
    drag = 0.5 * drag_coefficient * rho * velocity**2
    
    # 上昇力
    lift_force = lift - weight
    
    # 加速度
    acceleration = (thrust - drag) / mass
    
    return [velocity, acceleration]

• 航空機の運動方程式を表す関数です。
  引数として状態変数 state、時間 t、航空機のパラメータ（推力、抗力係数、揚力係数、質量、重力加速度）を取ります。
• 関数内で、高度と速度から揚力、重力、抗力を計算し、運動方程式に基づいて加速度を計算します。
• 計算結果として、速度と加速度のリストを返します。

3. パラメータの設定

1
2
3
4
5

thrust = 50000  # 推力（N）
drag_coefficient = 0.3  # 抗力係数
lift_coefficient = 0.4  # 揚力係数
mass = 5000  # 航空機の質量（kg）
gravity = 9.81  # 重力加速度（m/s^2）

• 航空機のパラメータを設定します。

推力、抗力係数、揚力係数、航空機の質量、重力加速度が含まれます。

4. 初期条件の設定

1
2
3

initial_altitude = 0  # 初期高度（m）
initial_velocity = 0  # 初期速度（m/s）
initial_state = [initial_altitude, initial_velocity]

• 初期高度と初期速度を設定し、状態変数の初期値としてリストに格納します。

5. 時間の範囲の設定

1

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

• 時間の範囲を設定します。
  ここでは$0$から$100$までの範囲を等間隔で$1000$分割しています。

6. 運動方程式の解析

1
2

result = odeint(aircraft_motion_equation, initial_state, t, args=(thrust, drag_coefficient, lift_coefficient, mass, gravity))
altitude, velocity = result.T

• odeint を使用して、航空機の運動方程式を解析します。
  初期状態、時間、およびパラメータを引数として渡します。
• 解析結果から高度と速度を取得します。

7. グラフのプロット

1
2
3
4
5
6
7
8

plt.plot(t, altitude, label='Altitude')
plt.plot(t, velocity, label='Velocity')
plt.xlabel('Time (s)')
plt.ylabel('Altitude (m) / Velocity (m/s)')
plt.title('Aircraft Motion')
plt.legend()
plt.grid(True)
plt.show()

• matplotlib.pyplot を使用して、時間に対する高度と速度の変化をプロットします。
• $X軸$に時間（秒）、$Y軸$に高度（メートル）と速度（メートル/秒）を表示します。
• グラフのタイトルやラベルを設定し、凡例とグリッドを追加しています。

結果解説

[実行結果]

上記のコードで生成されるグラフには、時間と高度、速度の関係が表示されます。

• 横軸は時間（秒）を表し、縦軸には高度（メートル）と速度（メートル/秒）が表示されます。
• 「Altitude」という線は、航空機の高度の変化を示しています。
  航空機の高度が時間とともにどのように変化するかを示しています。
• 「Velocity」という線は、航空機の速度の変化を示しています。
  航空機の速度が時間とともにどのように変化するかを示しています。

グラフの挙動を詳しく説明すると：

• 初期状態では、高度と速度はともにゼロです。
• 時間が経つにつれて、航空機の推力が働くことで高度が上昇し、速度が増加します。
• 一定の高度に達した後、航空機は水平飛行に移行し、高度は一定になりますが、速度は一定にはなりません。
  速度は、推力と抗力のバランスによって決まります。
• 最終的に、航空機が安定した速度に達し、高度も一定になります。

このように、グラフは航空機の運動を時間とともに可視化しています。

感染症の伝播（SIRモデル）

病気に関する方程式の例として、感染症の伝播を表すSIRモデルを取り上げます。

このモデルは、人口を感受性（S）、感染（I）、回復（R）の3つのカテゴリに分類し、これらのカテゴリ間で移行が起こると仮定します。

SIRモデルの方程式は以下の通りです：

$$
\frac{dS}{dt} = -\frac{\beta SI}{N}
$$

$$
\frac{dI}{dt} = \frac{\beta SI}{N} - \gamma I
$$

$$
\frac{dR}{dt} = \gamma I
$$

ここで、$ (S) $は感受性保持者の数、$ (I) $は感染者の数、$ (R) $は回復者の数、$ (N) $は総人口の数、$ (\beta) $は感染率、$ (\gamma) $は回復率です。

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

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

# SIRモデルの微分方程式
def sir_model(y, t, beta, gamma, N):
    S, I, R = y
    dSdt = -beta * S * I / N
    dIdt = beta * S * I / N - gamma * I
    dRdt = gamma * I
    return [dSdt, dIdt, dRdt]

# 初期値
N = 1000  # 総人口
I0 = 1    # 初期感染者数
R0 = 0    # 初期回復者数
S0 = N - I0 - R0  # 初期感受性保持者数
beta = 0.3  # 感染率
gamma = 0.1  # 回復率

# 時間の設定
t = np.linspace(0, 160, 160)

# SIRモデルを解く
solution = odeint(sir_model, [S0, I0, R0], t, args=(beta, gamma, N))
S, I, R = solution.T

# グラフ化
plt.figure(figsize=(10, 6))
plt.plot(t, S, label='Susceptible')
plt.plot(t, I, label='Infected')
plt.plot(t, R, label='Recovered')
plt.xlabel('Time (days)')
plt.ylabel('Population')
plt.title('SIR Model')
plt.legend()
plt.grid(True)
plt.show()

このコードでは、初期の感染者数と回復者数を設定し、微分方程式を解いて各カテゴリの人数を求め、それをグラフ化しています。

[実行結果]

ソースコード解説

このソースコードは、感染症の伝播をモデル化するSIRモデルを使用しています。

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

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

• import numpy as np:
  数値計算を行うためのライブラリであるNumPyをインポートします。
  これにより、配列操作や数学的な計算が可能になります。
• from scipy.integrate import odeint:
  微分方程式を解くためのodeint関数をSciPyからインポートします。
• import matplotlib.pyplot as plt:
  データの可視化のためのMatplotlibをインポートします。

2. SIRモデルの微分方程式の定義:

• def sir_model(y, t, beta, gamma, N):
  SIRモデルの微分方程式を定義します。
  この関数は、現在の状態$y$と時間$t$に対して微分方程式を解きます。
  また、感染率（$beta$）、回復率（$gamma$）、総人口（$N$）のパラメータも受け取ります。

3. 初期値の設定:

• 初期感染者数（$I0$）、初期回復者数（$R0$）などの各種パラメータが設定されています。

4. 時間の設定:

• t = np.linspace(0, 160, 160):
  $0$から$160$までの範囲を$160$個の点に分割して、時間の配列を作成します。
  これにより、微分方程式が時間に沿って解かれます。

5. 微分方程式の解析:

• solution = odeint(sir_model, [S0, I0, R0], t, args=(beta, gamma, N)):
  odeint関数を使用してSIRモデルの微分方程式を解きます。
  初期条件として、感受性保持者数（$S0$）、感染者数（$I0$）、回復者数（$R0$）を与え、時間t、感染率（$beta$）、回復率（$gamma$）、総人口（$N$）を引数として渡します。

6. グラフ化:

• 解析結果をプロットして、感受性保持者、感染者、回復者の人数の時間変化を可視化します。
  それぞれのプロットは異なる色で示され、凡例にはそれぞれのカテゴリ名が表示されます。
• plt.xlabel, plt.ylabel, plt.titleを使用して、軸ラベルやタイトルを追加します。
• plt.legend()を使って凡例を表示し、plt.grid(True)でグリッド線を表示します。
• plt.show()でグラフを表示します。

このコードは、SIRモデルを使って感染症の伝播を数値的に解析し、結果をグラフ化するものです。

結果解説

[実行結果]

このグラフは、SIRモデル（感受性(S)、感染(I)、回復(R)）に基づいて感染症の伝播をモデル化した結果を示しています。

時間軸 (Time):

グラフの横軸は時間を表し、単位は日数です。
感染の伝播が時間の経過とともにどのように変化するかを示します。

人口 (Population):

グラフの縦軸は人口を表します。
感受性保持者、感染者、回復者の数が表示されます。

感受性保持者 (Susceptible):

曲線”Susceptible”は感受性保持者の数を示します。
つまり、まだ感染していない人々の数です。
感染者との接触によって感染のリスクにさらされています。

感染者 (Infected):

曲線”Infected”は感染者の数を示します。
感染者は他の人々に病気を伝染させる可能性があります。
感染者の数が増えると、感染が拡大していることを示します。

回復者 (Recovered):

曲線”Recovered”は回復者の数を示します。
回復者は感染から回復し、免疫を獲得した人々です。
回復者の数が増えると、感染の広がりが鈍化し、感染者数が減少することを示します。

ピーク (Peak):

“Infected”曲線の頂点が感染のピークを表します。
ピーク時に感染者数が最大となり、その後減少し始めます。

このグラフを通じて、感染症の伝播パターンや感染のピーク時期、回復の進行など、感染症の伝播プロセスに関する情報が視覚的に理解できます。

天気に関する方程式

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

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

$$
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 ) $を表しています。

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

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

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

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