フォッカー-プランク方程式

フォッカー-プランク方程式

フォッカー-プランク方程式は、確率過程の時間発展を記述する偏微分方程式です。

例えば、ブラウン運動などの確率過程の確率密度関数の時間変化を記述します。

1次元のフォッカー-プランク方程式は以下のように表されます。

$$
∂P(x,t)/∂t = -∂/∂x[D(1)(x,t)P(x,t)] + (1/2)∂^2/∂x^2[D(2)(x,t)P(x,t)]
$$

ここで、$P(x,t)$は確率密度関数、$D(1)(x,t)$は漂流係数、$D(2)(x,t)$は拡散係数です。

この方程式を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
40
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

# 定数
D1 = 0.1  # 漂流係数
D2 = 0.01  # 拡散係数
X_MAX = 10.0  # x範囲
T_MAX = 10.0  # 時間範囲
N_X = 100  # x方向の分割数
N_T = 1000  # 時間ステップ数

# 配列の初期化
dx = X_MAX / (N_X - 1)
dt = T_MAX / (N_T - 1)
x = np.linspace(0, X_MAX, N_X)
t = np.linspace(0, T_MAX, N_T)
P = np.zeros((N_T, N_X))

# 初期条件
P[0] = np.exp(-(x - 3)**2 / 0.5)  # ガウス分布

# フォッカー-プランク方程式の解法
for i in range(N_T - 1):
    for j in range(1, N_X - 1):
        P[i+1, j] = P[i, j] \
                    - dt * (D1 * (P[i, j+1] - P[i, j-1]) / (2 * dx)) \
                    + (D2 * dt / (2 * dx**2)) * (P[i, j+1] - 2 * P[i, j] + P[i, j-1])
    P[i+1, 0] = P[i+1, 1]  # 端の条件
    P[i+1, -1] = P[i+1, -2]  # 端の条件

# 結果のプロット
fig, ax = plt.subplots(figsize=(8, 6))
for i in range(0, N_T, 100):
    ax.plot(x, P[i], label=f't = {t[i]:.2f}')
ax.set_xlabel('x')
ax.set_ylabel('P(x,t)')
ax.set_title('Fokker-Planck Equation')
ax.legend()
plt.show()

このコードでは、フォッカー-プランク方程式を有限差分法で近似的に解いています。
初期条件としてガウス分布を与え、時間発展させています。

出力される図は、確率密度関数$P(x,t)$の時間変化を示しています。
時間が経つにつれて、分布が拡散していく様子が確認できます。

[実行結果]

ソースコード解説

ここでは、フォッカー-プランク方程式の数値解を計算し、その結果を可視化するPythonコードについて説明します。

1. ライブラリのインポートと設定

1
2
3
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

NumPyとMatplotlibライブラリをインポートし、Matplotlibのインラインモードを有効にしています。

2. 定数の設定

1
2
3
4
5
6
D1 = 0.1  # 漂流係数
D2 = 0.01  # 拡散係数
X_MAX = 10.0  # x範囲
T_MAX = 10.0  # 時間範囲
N_X = 100  # x方向の分割数
N_T = 1000  # 時間ステップ数

フォッカー-プランク方程式の係数$(D1, D2)$、計算範囲$(X_MAX, T_MAX)$、空間と時間の分割数$(N_X, N_T)$を設定しています。

3. 配列の初期化

1
2
3
4
5
dx = X_MAX / (N_X - 1)
dt = T_MAX / (N_T - 1)
x = np.linspace(0, X_MAX, N_X)
t = np.linspace(0, T_MAX, N_T)
P = np.zeros((N_T, N_X))

空間と時間の刻み幅$(dx, dt)$、座標と時間の配列$(x, t)$、確率密度関数の配列$P(N_T x N_X)$を初期化しています。

4. 初期条件の設定

1
P[0] = np.exp(-(x - 3)**2 / 0.5)  # ガウス分布

初期の確率密度関数にガウス分布(平均値$3$、分散$0.5$)を設定しています。

5. フォッカー-プランク方程式の解法

1
2
3
4
5
6
7
for i in range(N_T - 1):
    for j in range(1, N_X - 1):
        P[i+1, j] = P[i, j] \
                    - dt * (D1 * (P[i, j+1] - P[i, j-1]) / (2 * dx)) \
                    + (D2 * dt / (2 * dx**2)) * (P[i, j+1] - 2 * P[i, j] + P[i, j-1])
    P[i+1, 0] = P[i+1, 1]  # 端の条件
    P[i+1, -1] = P[i+1, -2]  # 端の条件

フォッカー-プランク方程式を有限差分法で近似的に解いています。

時間ループとx方向のループを使って、各時間ステップと空間格子点における確率密度関数の値を計算しています。
また、端の条件も設定しています。

6. 結果のプロット

1
2
3
4
5
6
7
8
fig, ax = plt.subplots(figsize=(8, 6))
for i in range(0, N_T, 100):
    ax.plot(x, P[i], label=f't = {t[i]:.2f}')
ax.set_xlabel('x')
ax.set_ylabel('P(x,t)') 
ax.set_title('Fokker-Planck Equation')
ax.legend()
plt.show()

最後に、計算結果をプロットしています。
時間ステップごとに確率密度関数$P(x,t)$をプロットし、各曲線に対応する時刻tを凡例に表示しています。
また、軸ラベルとタイトルを設定し、グラフを表示しています。

このコードでは、フォッカー-プランク方程式の数値解を計算し、確率密度関数の時間発展を可視化しています。
初期条件にガウス分布を与え、時間の経過とともに確率密度がどのように拡散していくかを示しています。

結果解説

[実行結果]

出力されるグラフは、確率密度関数 $P(x,t)$の時間発展を示しています。

  • 横軸は$x$で、確率変数の値を表します。
  • 縦軸は$P(x,t)$で、その時刻$t$における確率密度を表します。

グラフには複数の曲線が描かれています。
それぞれの曲線は異なる時刻$t$における確率密度関数$P(x,t)$を表しています。

  • 初期条件($t=0$)の曲線は、平均値が$3$、分散が$0.5$のガウス分布になっています。
  • 時間が経つにつれ、曲線が広がっていく様子が見て取れます。
    これは確率密度関数が拡散していることを表しています。
  • 端($x=0$および$x=X_MAX$)の条件を満たすため、曲線の両端は$0$に近づくように調整されています。

凡例には、各曲線に対応する時刻tが表示されています。
例えば、最初の曲線は$t=0.00$、次の曲線は$t=1.00$というように時間が経過しています。

このグラフから、次のような情報が読み取れます:

  1. 初期の確率分布の形状
  2. 時間の経過とともに確率分布がどのように変化するか
  3. 確率密度がどのように拡散し、平らになっていくか
  4. 端の条件がどのように満たされているか

フォッカー-プランク方程式は、このように確率過程の時間発展を記述する役割があり、このグラフはその数値解の可視化結果になっています。