コップの中にインクを一滴垂らすと、インクは徐々に広がっていきます。このとき、任意の時刻 $t$ で「インクの分子がある位置 $x$ にいる確率」はどのように変化するでしょうか。個々の分子の軌道はブラウン運動のようにランダムですが、多数の分子の「確率密度」は決定論的な偏微分方程式に従います。この偏微分方程式がフォッカー・プランク方程式（Fokker-Planck equation）です。
フォッカー・プランク方程式は、確率微分方程式（SDE）が記述する個々の粒子の「ミクロな視点」を、確率密度関数の時間発展という「マクロな視点」に変換する橋渡しの役割を果たします。
- 統計力学: ランジュバン方程式からボルツマン分布への導出に使われます。系が熱平衡に達する過程を定量的に記述できます
- 金融工学: オプション価格を求めるブラック・ショールズ方程式は、フォッカー・プランク方程式の双対（コルモゴロフ後退方程式）に他なりません
- 生物物理学: イオンチャネルの開閉やタンパク質のフォールディングなど、ポテンシャル中のブラウン運動の記述に広く使われます
- 制御工学: 不確実性を含む系の確率分布の伝搬を予測するために不可欠です
本記事では、フォッカー・プランク方程式をSDEから体系的に導出し、代表的な場合の解を求め、Pythonによる数値解で理論を検証します。
本記事の内容
- SDEからフォッカー・プランク方程式への導出
- 自由拡散（ブラウン運動）のフォッカー・プランク方程式と解
- ドリフト付き拡散のフォッカー・プランク方程式
- OU過程の定常分布と緩和過程
- ダブルウェルポテンシャル中の拡散
- Pythonによる数値解法と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 伊藤の公式の応用 — 伊藤の公式とSDEの解法
- 伊藤積分の定義と性質 — 伊藤積分の基礎
- 偏微分方程式 — 拡散方程式の基本
ミクロからマクロへ — 二つの視点
個々の粒子の視点（SDE）
確率微分方程式は、個々の粒子の軌道を記述します。1次元の一般的なSDEは次の形をしています。
$$ \begin{equation} dX_t = \mu(X_t, t)\,dt + \sigma(X_t, t)\,dB_t \end{equation} $$
ここで $\mu(x, t)$ はドリフト係数（決定論的な力）、$\sigma(x, t)$ は拡散係数（ノイズの強度）です。このSDEを解くと、各実現値 $\omega$ に対して一本のサンプルパス $X_t(\omega)$ が得られます。
しかし、個々のパスはランダムであり、「典型的にどう振る舞うか」を理解するには統計的な記述が必要です。
確率密度の視点（フォッカー・プランク方程式）
SDEに従う粒子が多数ある（あるいは一つの粒子の実験を多数回繰り返す）とき、時刻 $t$ で粒子が位置 $x$ の近傍にいる確率は確率密度関数 $p(x, t)$ で記述されます。この確率密度関数がどのように時間発展するかを決める方程式がフォッカー・プランク方程式です。
イメージとしては、SDEが「一匹の酔っ払いの千鳥足」を記述するのに対し、フォッカー・プランク方程式は「大勢の酔っ払いが街に散らばっていく様子」を記述します。個々の酔っ払いの動きはランダムですが、群衆全体の密度分布は決定論的な法則に従うのです。
では、SDEからフォッカー・プランク方程式を導出しましょう。
SDEからフォッカー・プランク方程式の導出
導出の方針
導出の基本的なアイデアは、任意の滑らかな試験関数 $\phi(x)$ に対して $E[\phi(X_t)]$ の時間微分を二つの方法で計算し、等置することです。
一方では伊藤の公式を使い、他方では確率密度関数を使って期待値を書きます。
伊藤の公式による計算
$\phi(x)$ が十分滑らかな関数であるとき、伊藤の公式を $\phi(X_t)$ に適用すると、
$$ d\phi(X_t) = \phi'(X_t)\,dX_t + \frac{1}{2}\phi”(X_t)(dX_t)^2 $$
$dX_t = \mu\,dt + \sigma\,dB_t$ と $(dX_t)^2 = \sigma^2\,dt$ を代入すると、
$$ d\phi(X_t) = \left[\mu\phi'(X_t) + \frac{1}{2}\sigma^2\phi”(X_t)\right]dt + \sigma\phi'(X_t)\,dB_t $$
両辺の期待値をとります。伊藤積分の期待値はゼロ（マルチンゲール性）なので、
$$ \frac{d}{dt}E[\phi(X_t)] = E\left[\mu(X_t, t)\phi'(X_t) + \frac{1}{2}\sigma^2(X_t, t)\phi”(X_t)\right] $$
確率密度関数による表現
期待値を確率密度関数 $p(x, t)$ を使って書き換えると、
$$ \frac{d}{dt}\int_{-\infty}^{\infty}\phi(x)p(x, t)\,dx = \int_{-\infty}^{\infty}\left[\mu(x, t)\phi'(x) + \frac{1}{2}\sigma^2(x, t)\phi”(x)\right]p(x, t)\,dx $$
左辺では微分と積分を交換し（正当性は適切な条件のもとで保証される）、
$$ \int_{-\infty}^{\infty}\phi(x)\frac{\partial p}{\partial t}(x, t)\,dx = \int_{-\infty}^{\infty}\left[\mu(x, t)\phi'(x) + \frac{1}{2}\sigma^2(x, t)\phi”(x)\right]p(x, t)\,dx $$
右辺を部分積分で変形します。まず第1項について、
$$ \int \mu\phi’p\,dx = -\int \phi\frac{\partial}{\partial x}(\mu p)\,dx $$
境界での寄与はゼロ（$\phi$ がコンパクトサポートを持つか、$p$ が十分速く減衰する場合）とします。
第2項も2回の部分積分を行います。
$$ \int \frac{1}{2}\sigma^2\phi”p\,dx = \int \phi\frac{\partial^2}{\partial x^2}\left(\frac{1}{2}\sigma^2 p\right)dx $$
以上をまとめると、
$$ \int \phi(x)\frac{\partial p}{\partial t}\,dx = \int \phi(x)\left[-\frac{\partial}{\partial x}(\mu p) + \frac{\partial^2}{\partial x^2}\left(\frac{1}{2}\sigma^2 p\right)\right]dx $$
$\phi$ は任意の試験関数なので、被積分関数が恒等的に等しくなければなりません。
$$ \begin{equation} \frac{\partial p}{\partial t} = -\frac{\partial}{\partial x}\left[\mu(x, t)\,p\right] + \frac{\partial^2}{\partial x^2}\left[\frac{1}{2}\sigma^2(x, t)\,p\right] \end{equation} $$
これがフォッカー・プランク方程式（前進コルモゴロフ方程式とも呼ばれる）です。
方程式の各項の物理的意味
フォッカー・プランク方程式の右辺は二つの項からなります。
第1項: $-\frac{\partial}{\partial x}(\mu p)$ — ドリフト項（advection term）。確率密度をドリフト $\mu$ の方向に輸送する効果を表します。確率の流れ（probability current）$J_{\text{drift}} = \mu p$ の発散のマイナスです。
第2項: $\frac{\partial^2}{\partial x^2}\left(\frac{1}{2}\sigma^2 p\right)$ — 拡散項（diffusion term）。ノイズによって確率密度が広がる（拡散する）効果を表します。拡散係数 $D(x) = \frac{1}{2}\sigma^2(x)$ が定数の場合は通常の拡散方程式 $D\frac{\partial^2 p}{\partial x^2}$ に帰着します。
この方程式は確率の保存則（連続の方程式）の形をしています。確率の流れを $J = \mu p – \frac{\partial}{\partial x}\left(\frac{1}{2}\sigma^2 p\right)$ と定義すると、
$$ \frac{\partial p}{\partial t} + \frac{\partial J}{\partial x} = 0 $$
と書け、確率が生成も消滅もしないことを表しています。
導出の核心は、伊藤の公式の伊藤補正項 $\frac{1}{2}\sigma^2\phi”$ が、フォッカー・プランク方程式の拡散項に対応していることです。通常の微積分（ストラトノビッチ型）を使うと、異なる形のフォッカー・プランク方程式が得られます。
具体的な問題でフォッカー・プランク方程式を解いてみましょう。まずは最も基本的な自由拡散（ブラウン運動）の場合です。
自由拡散のフォッカー・プランク方程式
問題の設定
最も単純な場合として、ドリフトのない純粋な拡散過程を考えます。
$$ dX_t = \sigma\,dB_t, \quad X_0 = x_0 $$
ここで $\sigma$ は定数です。$\mu = 0$, $\sigma = \text{const}$ のフォッカー・プランク方程式は、
$$ \begin{equation} \frac{\partial p}{\partial t} = D\frac{\partial^2 p}{\partial x^2}, \quad D = \frac{\sigma^2}{2} \end{equation} $$
これは物理学で拡散方程式（heat equation）として知られる偏微分方程式です。
解と物理的意味
初期条件 $p(x, 0) = \delta(x – x_0)$（時刻0で粒子が $x_0$ にいる）のもとでの解は、フーリエ変換を使って次のように求まります。
$$ \begin{equation} p(x, t) = \frac{1}{\sqrt{4\pi Dt}}\exp\left(-\frac{(x – x_0)^2}{4Dt}\right) \end{equation} $$
これは平均 $x_0$、分散 $2Dt = \sigma^2 t$ の正規分布です。時間が経つにつれて分布は広がり（$\sigma_x = \sqrt{2Dt}$ が増加し）、最終的にはすべての空間に一様に広がっていきます。
この解の物理的意味を整理しましょう。
- 期待値: $E[X_t] = x_0$ — 粒子の平均位置は動かない（ドリフトがないため）
- 分散: $\text{Var}(X_t) = 2Dt$ — 分散は時間に比例して増大する
- 二乗平均変位: $\langle(X_t – x_0)^2\rangle = 2Dt$ — アインシュタインの関係式
最後の関係はアインシュタインが1905年のブラウン運動の論文で導いた結果であり、ブラウン運動の拡散係数と粒子の揺らぎを結びつける歴史的に重要な公式です。
自由拡散では分布は際限なく広がっていきますが、外力（ドリフト）がある場合はどうでしょうか。次にドリフト付きの場合を見てみましょう。
OU過程のフォッカー・プランク方程式と定常分布
フォッカー・プランク方程式の具体的な形
伊藤の公式の応用で解いたオルンシュタイン・ウーレンベック過程は、
$$ dX_t = -\theta(X_t – \mu_0)\,dt + \sigma\,dB_t $$
のSDEに従います。ここで $\theta > 0$ は回帰速度、$\mu_0$ は長期平均です。$\mu(x) = -\theta(x – \mu_0)$, $\sigma = \text{const}$ を代入すると、フォッカー・プランク方程式は次のようになります。
$$ \begin{equation} \frac{\partial p}{\partial t} = \theta\frac{\partial}{\partial x}[(x – \mu_0)p] + D\frac{\partial^2 p}{\partial x^2} \end{equation} $$
ここで $D = \sigma^2/2$ です。右辺の第1項 $\theta\frac{\partial}{\partial x}[(x – \mu_0)p]$ は、確率密度を $\mu_0$ に向かって引き戻す復元力の効果を表しています。
定常分布の導出
十分長い時間が経過すると、系は定常状態に達し、$\frac{\partial p}{\partial t} = 0$ となります。定常分布 $p_{\text{ss}}(x)$ は次の常微分方程式を満たします。
$$ 0 = \theta\frac{d}{dx}[(x – \mu_0)p_{\text{ss}}] + D\frac{d^2 p_{\text{ss}}}{dx^2} $$
一度 $x$ で積分すると（境界条件 $J \to 0$ as $x \to \pm\infty$）、
$$ 0 = \theta(x – \mu_0)p_{\text{ss}} + D\frac{dp_{\text{ss}}}{dx} $$
変数分離して解きます。
$$ \frac{dp_{\text{ss}}}{p_{\text{ss}}} = -\frac{\theta}{D}(x – \mu_0)\,dx $$
両辺を積分すると、
$$ \ln p_{\text{ss}} = -\frac{\theta}{2D}(x – \mu_0)^2 + C $$
指数をとり、正規化定数を求めると、
$$ \begin{equation} p_{\text{ss}}(x) = \sqrt{\frac{\theta}{2\pi D}}\exp\left(-\frac{\theta}{2D}(x – \mu_0)^2\right) = \frac{1}{\sqrt{2\pi\sigma_{\text{ss}}^2}}\exp\left(-\frac{(x – \mu_0)^2}{2\sigma_{\text{ss}}^2}\right) \end{equation} $$
ここで $\sigma_{\text{ss}}^2 = D/\theta = \sigma^2/(2\theta)$ です。これは伊藤の公式の応用でSDEから直接導いた定常分散と一致しています。
ボルツマン分布との関係
定常分布をポテンシャルの観点から見ると、深い物理的意味が浮かび上がります。OU過程のドリフトは、二次ポテンシャル $V(x) = \frac{\theta}{2}(x – \mu_0)^2$ の勾配から生じる力 $\mu(x) = -V'(x) = -\theta(x – \mu_0)$ に対応しています。
定常分布は次のように書けます。
$$ p_{\text{ss}}(x) \propto \exp\left(-\frac{V(x)}{D}\right) = \exp\left(-\frac{2V(x)}{\sigma^2}\right) $$
これは統計力学のボルツマン分布 $p \propto e^{-V/(k_B T)}$ と同じ形をしています。拡散係数 $D = \sigma^2/2$ が温度 $k_B T$ の役割を果たしているのです。この対応関係は、アインシュタインの揺動散逸定理（fluctuation-dissipation theorem）の表れであり、系の散逸（ドリフトによる減衰）と揺動（ノイズ）が熱平衡で釣り合うことを示しています。
一般のポテンシャルでの定常分布
一般のポテンシャル $V(x)$ 中の過剰減衰（overdamped）ランジュバン方程式、
$$ dX_t = -V'(X_t)\,dt + \sigma\,dB_t $$
の定常分布も同様にボルツマン分布で与えられます。
$$ \begin{equation} p_{\text{ss}}(x) = \frac{1}{Z}\exp\left(-\frac{2V(x)}{\sigma^2}\right), \quad Z = \int_{-\infty}^{\infty}\exp\left(-\frac{2V(x)}{\sigma^2}\right)dx \end{equation} $$
$Z$ は分配関数（正規化定数）です。この結果は、ポテンシャルが二次でなくても成り立ちます。
ここまでの理論をPythonで数値的に検証しましょう。フォッカー・プランク方程式を数値的に解き、SDEのモンテカルロシミュレーションと比較します。
Pythonによる数値解法と可視化
自由拡散の時間発展
まず、最も基本的な自由拡散のフォッカー・プランク方程式を差分法で解き、解析解と比較します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
# --- パラメータ ---
sigma = 1.0
D = sigma**2 / 2 # 拡散係数
x0 = 0.0 # 初期位置
T = 2.0
# --- 空間・時間の格子 ---
L = 6.0 # 空間の範囲 [-L, L]
nx = 400 # 空間の分割数
nt = 10000 # 時間の分割数
dx = 2 * L / nx
dt_fp = T / nt
x = np.linspace(-L, L, nx)
# 安定性条件の確認: dt < dx^2 / (2D)
assert dt_fp < dx**2 / (2 * D), "CFL condition not satisfied"
# --- 初期条件: ガウス近似のデルタ関数 ---
p = np.exp(-x**2 / (2 * 0.01)) / np.sqrt(2 * np.pi * 0.01)
p /= np.trapz(p, x) # 正規化
# --- 差分法による時間発展 ---
snapshots = {}
snapshot_times = [0.0, 0.1, 0.3, 0.5, 1.0, 2.0]
snapshots[0.0] = p.copy()
current_time = 0.0
for step in range(nt):
# 拡散項（中心差分）
d2p_dx2 = (np.roll(p, -1) - 2 * p + np.roll(p, 1)) / dx**2
# 境界条件: ゼロ
d2p_dx2[0] = 0
d2p_dx2[-1] = 0
# 時間更新
p = p + D * d2p_dx2 * dt_fp
current_time += dt_fp
# スナップショット保存
for ts in snapshot_times:
if abs(current_time - ts) < dt_fp / 2 and ts not in snapshots:
snapshots[ts] = p.copy()
fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))
# 左図: 数値解の時間発展
ax = axes[0]
colors = plt.cm.viridis(np.linspace(0, 1, len(snapshot_times)))
for i, ts in enumerate(snapshot_times):
if ts in snapshots:
label = f'$t = {ts}$'
ax.plot(x, snapshots[ts], color=colors[i], linewidth=1.5, label=label)
ax.set_xlabel('$x$', fontsize=12)
ax.set_ylabel('$p(x, t)$', fontsize=12)
ax.set_title('Free diffusion: FP equation (numerical)', fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_xlim(-5, 5)
# 右図: 数値解 vs 解析解
ax = axes[1]
for ts in [0.1, 0.5, 1.0, 2.0]:
if ts in snapshots:
ax.plot(x, snapshots[ts], '-', linewidth=1.5, label=f'Numerical $t={ts}$')
# 解析解
p_exact = norm.pdf(x, x0, np.sqrt(2 * D * ts))
ax.plot(x, p_exact, '--', linewidth=1.5, alpha=0.8)
ax.set_xlabel('$x$', fontsize=12)
ax.set_ylabel('$p(x, t)$', fontsize=12)
ax.set_title('Numerical (solid) vs analytical (dashed)', fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_xlim(-5, 5)
plt.tight_layout()
plt.savefig('fp_free_diffusion.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから、以下のことが確認できます。
-
左図: 確率密度が時間とともに広がっている。初期の鋭いピーク（デルタ関数の近似）が、時間とともにガウス分布として幅を広げていきます。ピークの高さは減少しますが、分布の下の面積（全確率）は常に1に保たれています。これは確率の保存則の帰結です
-
右図: 数値解と解析解が正確に一致している。実線（差分法による数値解）と破線（ガウス分布の解析解）が完全に重なっており、数値的手法が正しく機能していることが確認できます。分散 $\sigma_x^2 = 2Dt$ に比例して幅が広がることも読み取れます
OU過程の定常分布への緩和
次に、OU過程のフォッカー・プランク方程式を数値的に解き、定常分布への緩和過程を可視化します。同時にモンテカルロシミュレーションとも比較します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
np.random.seed(42)
# --- パラメータ ---
theta = 2.0
mu0 = 1.0
sigma = 0.8
D = sigma**2 / 2
T = 3.0
# --- 空間・時間の格子 ---
L = 5.0
nx = 500
nt = 50000
dx = 2 * L / nx
dt_fp = T / nt
x = np.linspace(-L + mu0, L + mu0, nx) # 長期平均を中心に配置
# --- 初期条件: x0 = 3.0 に集中 ---
x0 = 3.0
p = np.exp(-(x - x0)**2 / (2 * 0.02)) / np.sqrt(2 * np.pi * 0.02)
p /= np.trapz(p, x)
# --- FP方程式の数値解 ---
snapshots_fp = {}
snapshot_times = [0.0, 0.2, 0.5, 1.0, 2.0, 3.0]
snapshots_fp[0.0] = p.copy()
for step in range(nt):
# ドリフト項: theta * d/dx[(x - mu0) * p]
drift_flux = theta * (x - mu0) * p
d_drift = (np.roll(drift_flux, -1) - np.roll(drift_flux, 1)) / (2 * dx)
# 拡散項
d2p = (np.roll(p, -1) - 2 * p + np.roll(p, 1)) / dx**2
# 境界処理
d_drift[0] = 0; d_drift[-1] = 0
d2p[0] = 0; d2p[-1] = 0
# 更新
p = p + (d_drift + D * d2p) * dt_fp
p = np.maximum(p, 0) # 負の値を防止
# スナップショット
current_time = (step + 1) * dt_fp
for ts in snapshot_times:
if abs(current_time - ts) < dt_fp / 2 and ts not in snapshots_fp:
snapshots_fp[ts] = p.copy()
# --- モンテカルロシミュレーション ---
n_mc = 50000
dt_mc = 0.001
X_mc = np.full(n_mc, x0)
mc_snapshots = {}
current_time_mc = 0.0
for ts in sorted(snapshot_times):
while current_time_mc < ts - dt_mc / 2:
dB = np.random.normal(0, np.sqrt(dt_mc), n_mc)
X_mc = X_mc - theta * (X_mc - mu0) * dt_mc + sigma * dB
current_time_mc += dt_mc
mc_snapshots[ts] = X_mc.copy()
# --- 定常分布 ---
sigma_ss = sigma / np.sqrt(2 * theta)
p_stationary = norm.pdf(x, mu0, sigma_ss)
fig, axes = plt.subplots(2, 3, figsize=(16, 10))
for i, ts in enumerate(snapshot_times):
ax = axes[i // 3, i % 3]
# FP数値解
if ts in snapshots_fp:
ax.plot(x, snapshots_fp[ts], 'b-', linewidth=2, label='FP equation')
# MC ヒストグラム
if ts in mc_snapshots:
ax.hist(mc_snapshots[ts], bins=80, density=True, alpha=0.4,
color='steelblue', edgecolor='none', label='Monte Carlo')
# 定常分布
ax.plot(x, p_stationary, 'r--', linewidth=1.5, alpha=0.7, label='Stationary')
ax.set_title(f'$t = {ts}$', fontsize=13)
ax.set_xlabel('$x$', fontsize=11)
ax.set_ylabel('$p(x, t)$', fontsize=11)
ax.set_xlim(mu0 - 3, mu0 + 3)
ax.set_ylim(0, 5)
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
plt.suptitle(rf'OU process: $\theta={theta}$, $\mu_0={mu0}$, $\sigma={sigma}$, $x_0={x0}$',
fontsize=14, y=1.01)
plt.tight_layout()
plt.savefig('fp_ou_relaxation.png', dpi=150, bbox_inches='tight')
plt.show()
この6枚のスナップショットから、OU過程の確率密度の緩和過程が詳細に確認できます。
-
$t = 0$: 初期状態。確率密度は $x_0 = 3.0$ に集中した鋭いピークです。定常分布（赤破線）は $\mu_0 = 1.0$ を中心としたガウス分布であり、初期分布とは大きくかけ離れています
-
$t = 0.2$: 拡散と回帰の開始。ピークが広がり始めると同時に、長期平均 $\mu_0$ の方向に移動しています。ドリフト（復元力）と拡散（ノイズ）の両方の効果が見えます
-
$t = 0.5$: 遷移期。分布はかなり広がり、中心が $\mu_0$ に近づいています。フォッカー・プランク方程式の数値解（青実線）とモンテカルロシミュレーション（ヒストグラム）が良く一致しています
-
$t = 1.0 \sim 3.0$: 定常状態への収束。分布が定常分布（赤破線）にほぼ一致しています。回帰速度 $\theta = 2$ の場合、特性時間は $1/\theta = 0.5$ なので、$t = 2$ では $e^{-2\theta \cdot 2} = e^{-8} \approx 0.0003$ となり、初期条件の記憶はほぼ完全に失われています
-
FP方程式の数値解とモンテカルロシミュレーションが全時刻で一致しており、SDEとフォッカー・プランク方程式が同じ物理を異なる視点から記述していることが数値的に確認できます
ダブルウェルポテンシャル中の拡散
最後に、より物理的に興味深い例として、ダブルウェルポテンシャル $V(x) = (x^2 – 1)^2$ 中のブラウン運動を考えます。このポテンシャルは $x = \pm 1$ に二つの極小値を持ち、$x = 0$ に極大値（エネルギー障壁）を持ちます。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
# --- ダブルウェルポテンシャル ---
def V(x):
return (x**2 - 1)**2
def dV(x):
return 4 * x * (x**2 - 1)
# --- パラメータ ---
sigma = 0.7
D = sigma**2 / 2
T = 10.0
# --- モンテカルロシミュレーション ---
n_mc = 100000
dt_mc = 0.001
n_steps = int(T / dt_mc)
# 初期条件: x = -1 の井戸に集中
X = np.full(n_mc, -1.0)
snapshot_times_dw = [0.0, 0.5, 1.0, 2.0, 5.0, 10.0]
mc_snap_dw = {0.0: X.copy()}
for step in range(n_steps):
dB = np.random.normal(0, np.sqrt(dt_mc), n_mc)
X = X - dV(X) * dt_mc + sigma * dB
ct = (step + 1) * dt_mc
for ts in snapshot_times_dw:
if abs(ct - ts) < dt_mc / 2 and ts not in mc_snap_dw:
mc_snap_dw[ts] = X.copy()
# --- 定常分布（ボルツマン分布） ---
x_range = np.linspace(-2.5, 2.5, 500)
p_boltz = np.exp(-2 * V(x_range) / sigma**2)
p_boltz /= np.trapz(p_boltz, x_range)
fig, axes = plt.subplots(2, 3, figsize=(16, 10))
for i, ts in enumerate(snapshot_times_dw):
ax = axes[i // 3, i % 3]
if ts in mc_snap_dw:
ax.hist(mc_snap_dw[ts], bins=100, density=True, alpha=0.6,
color='steelblue', edgecolor='none', label='Monte Carlo')
ax.plot(x_range, p_boltz, 'r-', linewidth=2, label='Boltzmann', alpha=0.8)
# ポテンシャルを背景に表示
ax2 = ax.twinx()
ax2.plot(x_range, V(x_range), 'k--', linewidth=1, alpha=0.3)
ax2.set_ylabel('$V(x)$', fontsize=10, alpha=0.3)
ax2.set_ylim(0, 3)
ax.set_title(f'$t = {ts}$', fontsize=13)
ax.set_xlabel('$x$', fontsize=11)
ax.set_ylabel('$p(x, t)$', fontsize=11)
ax.set_xlim(-2.5, 2.5)
ax.legend(fontsize=9, loc='upper left')
ax.grid(True, alpha=0.3)
plt.suptitle(rf'Double-well potential: $V(x) = (x^2-1)^2$, $\sigma = {sigma}$', fontsize=14, y=1.01)
plt.tight_layout()
plt.savefig('fp_double_well.png', dpi=150, bbox_inches='tight')
plt.show()
ダブルウェルポテンシャルの結果から、以下の物理的に重要な現象が観察できます。
-
初期状態（$t = 0$）: すべての粒子が左の井戸 $x = -1$ に集中しています。右の井戸 $x = +1$ は完全に空です
-
短時間（$t = 0.5 \sim 1.0$）: 粒子は左の井戸内で拡散しますが、$x = 0$ のエネルギー障壁を越えて右の井戸に到達する粒子はほとんどいません。分布は左の井戸の局所的な平衡（局所ボルツマン分布）に近づきます
-
中間時間（$t = 2.0 \sim 5.0$）: 熱揺らぎによってエネルギー障壁を越える粒子が徐々に増え、右の井戸にも粒子が蓄積し始めます。これはクラマースの脱出率（Kramers escape rate）によって支配される確率的なバリア越え現象です
-
長時間（$t = 10.0$）: 分布がボルツマン分布（赤い実線）に収束しています。左右の井戸に等しく粒子が分布し、系は熱平衡に達しています。ポテンシャルの対称性から、二つのピークの高さは等しくなります
-
ノイズ強度 $\sigma$ の役割: $\sigma$ が小さいとバリア越えに指数的に長い時間がかかり、大きいと素早く平衡に達します。これはアレニウスの法則 $\tau \sim \exp(2\Delta V/\sigma^2)$ として定量化されます
コルモゴロフ方程式との関係
コルモゴロフ前進方程式と後退方程式
フォッカー・プランク方程式には「双対」となるもう一つの偏微分方程式が存在します。
コルモゴロフ前進方程式（forward Kolmogorov equation）は、フォッカー・プランク方程式そのものであり、確率密度 $p(x, t | x_0, t_0)$ の $x$ と $t$ に関する方程式です。
$$ \frac{\partial p}{\partial t} = -\frac{\partial}{\partial x}(\mu p) + \frac{1}{2}\frac{\partial^2}{\partial x^2}(\sigma^2 p) $$
コルモゴロフ後退方程式（backward Kolmogorov equation）は、遷移確率密度の初期位置 $x_0$ と初期時刻 $t_0$ に関する方程式です。
$$ \begin{equation} \frac{\partial p}{\partial t_0} + \mu(x_0)\frac{\partial p}{\partial x_0} + \frac{1}{2}\sigma^2(x_0)\frac{\partial^2 p}{\partial x_0^2} = 0 \end{equation} $$
「前進」と「後退」の名前は、時間の方向に由来します。前進方程式は時間を未来に向かって進める方程式であり、後退方程式は時間を過去に向かって遡る方程式です。
後退方程式は金融工学で特に重要です。ブラック・ショールズ方程式は実はコルモゴロフ後退方程式の一種であり、「満期時刻から遡って現在のオプション価格を求める」という問題設定に自然に対応しています。
ファインマン・カッツの公式
SDEと偏微分方程式の間には、ファインマン・カッツの公式（Feynman-Kac formula）というより一般的な対応関係が存在します。
SDE $dX_t = \mu\,dt + \sigma\,dB_t$ に対して、関数 $u(x, t) = E[\phi(X_T) | X_t = x]$（終端条件 $\phi$ のもとでの条件付き期待値）は、次の偏微分方程式を満たします。
$$ \frac{\partial u}{\partial t} + \mu\frac{\partial u}{\partial x} + \frac{1}{2}\sigma^2\frac{\partial^2 u}{\partial x^2} = 0, \quad u(x, T) = \phi(x) $$
この公式は、確率的な問題（SDE、期待値計算）と解析的な問題（PDE）を結びつける強力な橋渡しの役割を果たします。モンテカルロ法（確率的手法）と差分法（解析的手法）のどちらでも同じ問題を解けることの理論的基盤がファインマン・カッツの公式なのです。
まとめ
本記事では、フォッカー・プランク方程式の導出、解法、数値検証について解説しました。
- フォッカー・プランク方程式 $\frac{\partial p}{\partial t} = -\frac{\partial}{\partial x}(\mu p) + \frac{\partial^2}{\partial x^2}(\frac{\sigma^2}{2}p)$ は、SDEに対応する確率密度の時間発展を記述する偏微分方程式です。導出には伊藤の公式が中心的な役割を果たします
- 自由拡散の場合はガウス分布が解であり、分散は $2Dt$ に比例して増大します（アインシュタインの関係式）
- OU過程のフォッカー・プランク方程式は定常分布 $N(\mu_0, \sigma^2/(2\theta))$ に収束します。これはボルツマン分布の一般化であり、揺動散逸定理の表れです
- ダブルウェルポテンシャル中の拡散では、エネルギー障壁の越え方が熱揺らぎの強さに依存し、クラマースの脱出率で特徴づけられます
- コルモゴロフ後退方程式はフォッカー・プランク方程式の双対であり、ファインマン・カッツの公式を通じて確率論と偏微分方程式が結ばれます
次のステップとして、以下の記事も参考にしてください。
- マルチンゲール収束定理と応用 — マルチンゲール理論の深い結果
- 伊藤の公式の応用 — SDEの求解テクニック
- 確率過程の基礎 — 確率過程の基本的な枠組み