確率微分方程式

このページでは、簡単な確率微分方程式の数値計算について考えてみる。

ゆらぎを伴う微分方程式

ランダムなゆらぎ$\xi(t)$を考え、それによって動かされる粒子の位置$x(t)$が微分方程式 $$ \frac{d x(t)}{dt} = \xi(t) \tag{1} $$ で変化するようなモデルを考えてみよう。 ここで、 ゆらぎの平均は0、すなわち$i$番目のサンプルの揺らぎの時系列を$\xi_i(t)$とすれば $$ \left\langle \xi(t) \right\rangle = \lim_{N\to\infty} \frac{1}{N} \sum_{i=1}^N \xi_i(t) = 0 $$ で、自己相関が $$ \left\langle \xi(t) \xi(t+\tau) \right\rangle = \lim_{N\to\infty} \frac{1}{N} \sum_{i=1}^N \xi_i(t) \xi_i(t+\tau) = R \delta(\tau) $$ であるような場合を考える（$\delta(\tau)$はデルタ関数）。 そして、これらの関係は時間の原点には依らずに成り立つ（$\xi(t)$は定常的）と仮定する。

この種の確率的（ランダム）な変数を含む微分方程式は確率微分方程式（stochastic differential equation）と呼ばれている。 確率微分方程式を数値的に解くには、ゆらぎの部分の扱い方に、通常の常微分方程式には無かった注意が必要である。

確率微分方程式の差分化

ここで、時間$t$を$\Delta$毎の区間に $$ [0, \Delta),\ [\Delta, 2\Delta),\ [2\Delta, 3\Delta),\cdots,[k \Delta, (k+1)\Delta),\cdots $$ のように等分割して、 $k$番目の区間の揺らぎの累積（積分）を $$ q_k = \int_{k\Delta}^{(k+1)\Delta} \xi(s) ds $$ と定義する。 すると、そのサンプル平均は $$ \left\langle q_k \right\rangle = \lim_{N \to \infty} \frac{1}{N} \sum_{i=1}^{N} \int_{k\Delta}^{(k+1)\Delta} \xi_i(s) ds = \int_{k\Delta}^{(k+1)\Delta} \left\langle \xi(s) \right\rangle ds = 0 $$ 分散は $$ \begin{eqnarray} \left\langle q_k q_{k+\ell} \right\rangle & = & \lim_{N\to\infty} \frac{1}{N} \sum_{i=1}^{N} \left( \int_{k\Delta}^{(k+1)\Delta} \xi_i(s) ds \right) \left( \int_{(k+\ell)\Delta}^{(k+\ell+1)\Delta} \xi_i(t) dt \right) \\ & = & \lim_{N\to\infty} \frac{1}{N} \sum_{i=1}^{N} \left( \int_{k\Delta}^{(k+1)\Delta} \xi_i(s) ds \right) \left( \int_{k \Delta}^{(k+1)\Delta} \xi_i(t + \ell\Delta) dt \right) \\ & = & \int_{k\Delta}^{(k+1)\Delta} \int_{k \Delta}^{(k+1)\Delta} \left[ \frac{1}{N} \lim_{N\to\infty} \sum_{i=1}^{N} \xi_i(s) \xi_i(t + \ell\Delta) \right] ds dt \\ & = & \int_{k\Delta}^{(k+1)\Delta} \int_{k \Delta}^{(k+1)\Delta} \left\langle \xi(s) \xi(t + \ell\Delta) \right\rangle ds dt \\ & = & \int_{k\Delta}^{(k+1)\Delta} \int_{k \Delta}^{(k+1)\Delta} R \delta(t + \ell\Delta -s) dt ds \\ & = & R \int_{k\Delta}^{(k+1)\Delta} \left\{ \int_{k \Delta + \ell \Delta - s}^{(k+1)\Delta + \ell \Delta - s} \delta(u) du \right\} ds \\ \end{eqnarray} $$ となる。 ここで、デルタ関数の積分 $$ \int_{k \Delta + \ell \Delta - s}^{(k+1)\Delta + \ell \Delta - s} \delta(u) du $$ は、$\ell=0$の場合に限り、$k \Delta \le s \le (k+1)\Delta$の範囲で積分区間が$u=0$を「通過」すること、 および $\epsilon \gt 0$ に対して $\int_{0-\epsilon}^{0+\epsilon} \delta(u) du=1$ を考慮すると、 $$ \left\langle q_k q_{k+\ell} \right\rangle = R \Delta \delta_{0,\ell} $$ が得られる（$\delta_{0,\ell}$はクロネッカーのデルタ）。

オイラー・丸山の近似

式(1)の両辺を$k\Delta$から$(k+1)\Delta$まで積分すると $$ x(k\Delta + \Delta) - x(k\Delta) = q_k $$ となるが、$x(k\Delta)$をあらためて$x_k$と書けば $$ x_{k+1} = x_k + q_k $$ である。 ここで、分散$R$を持つ乱数を$\xi$とすると、$q_k$の分散はその$\Delta$倍（標準偏差は$\sqrt{\Delta}$倍）であるから、 $$ x_{k+1} = x_k + \sqrt{\Delta} \xi $$ という差分近似によって$\Delta$だけ先の$x(t)$を数値的に見積もることが可能となる（オイラー・丸山法）。

同様に $$ \frac{dx(t)}{dt} = f(x(t)) + \xi(t) $$ に対しては、 $$ x_{k+1} = x_k + f(x_k) \Delta + \sqrt{\Delta} \xi $$ となり、確率的な揺らぎの項は他と時間差分$\Delta$の次数が異なってくる。

　解説: 白色ノイズ

$\Delta$毎に変動し、その値が平均0分散$R$で分布するような無限に続くランダムな時系列 $v(t_i)$を考える。 ここで、$i$は何番目の区分か、を表すことにしよう。 そして、$i$番目と$j$番目（$i \neq j$)は無相関、すなわち $$ \left\langle v(t_i) v(t_j) \right\rangle = R \delta_{i,j} $$ とする。

次に、時系列信号の時間のスケールを1/2、時間あたりのエネルギーを2倍にするような変換を施してみる。 つまり、波形を細かくしても、「時間刻み」あたりのエネルギーは同じになるように調整するわけである。 そうするためには、グラフの横軸を1/2に、縦軸を$\sqrt{2}$倍に変形すれば良い（下図）。

こうした操作を繰り返すことで、その自己相関が $$ \left\langle v(t) v(t') \right\rangle = R \delta(t-t') $$ で、ランダムなスパイクが連続するような波形が最終的に得られ、白色雑音（ホワイトノイズ）と呼ばれている。

ウィーナー・ヒンチンの定理により、 時系列信号の自己相関関数をフーリエ変換すると信号のパワースペクトル（振動数あたりのエネルギー密度）が得られるので、これを白色雑音に適用すると、 $$ S(\omega) = \int_{-\infty}^{\infty} R \delta(\tau) e^{-i \omega \tau} d\tau = R $$ となり、全ての角振動数について、エネルギー密度が一定（$R$）となる。 そうすると、全ての振動数についてエネルギーを積分すると無限大となってしまうことになるが、 現実の問題では、髙振動数側のどこかに必ずカットオフが存在するため、信号のエネルギーが発散するようなことは無い。 全ての振動数で強度が一定となることから、光のスペクトルになぞらえられたのが「白色」の由来である。

　練習： 確率微分方程式

オイラー・丸山の方法で微分方程式を数値的に解くコードの例を以下に示した。

これを雛形にして、$x(0)=0$を初期条件とした複数の経路を重ね書きするコードを作成してみなさい。

# coding: utf-8

import math
import random
import matplotlib.pyplot as plt

R=1
t=0
dt=0.01
x=0

X=[]
T=[]

while t<10:
    X.append(x)
    T.append(t)
    x = x + math.sqrt(R*dt) * random.gauss(0,1)
    t += dt

plt.plot(T,X, color='red', linewidth=1.0)
plt.grid(True)
plt.xlabel('t')
plt.ylabel('x')
plt.show()

　解説:ブラウン運動のモデル

ブラウン運動を表現する確率微分方程式

花粉やコロイドなどの微粒子がブラウン運動する様子は、 $$ \begin{eqnarray} \dot{x}(t) & = & v(t) \\ m \dot{v}(t) & = & -\gamma v(t) + \xi(t) \tag{2} \end{eqnarray} $$ のように表される(ドットは$t$による微分）。ここで、$m$は微粒子の質量、$\gamma$は流体からの受ける粘性抵抗の程度（摩擦）を表す係数、 $\xi(t)$は周囲の分子から粒子が受ける揺動力を表す。

ブラウン粒子の平均二乗変位

ここで、時刻$t=0$で原点$x=0$にあった粒子が、どのように周囲に移動するかを考えてみる。

変位の自乗、すなわち$x^2$の変化に着目して、その二階微分を $$ \frac{d (x^2)}{d t} = 2 v^2 + 2x \dot{v} $$ から計算してみる。全体に$m$をかけてから右辺の$\dot{v}$のところに運動方程式(2)を代入すると $$ m \frac{d (x^2)}{d t} = 2 m v^2 + 2 x (-\gamma v + \xi) $$ となるが、$\frac{d x^2}{dt} = 2 xv$であることを思い出すと、各変数を二乗した量についての微分方程式 $$ m \frac{d (x^2)}{d t} = 2 m v^2 -\gamma \frac{d (x^2)}{dt} + 2 x \xi $$ を得る。

ここで、初期条件のみ同じにして、同様の沢山の粒子についての平均を取ってみる。ある量 $\cdot$ の平均を $\langle \cdot \rangle$で表すことにして、 上記の微分方程式全体の平均を取ってみると $$ m \frac{d \langle x^2 \rangle}{d t} = 2 m \langle v^2 \rangle -\gamma \frac{d \langle x^2 \rangle}{dt} $$ ただし、揺動力が粒子の位置とは無相関であり、$\xi$の平均は0であるから、それらの積の平均は $\langle x \xi \rangle = \langle x \rangle \langle \xi \rangle= 0$とした。

コロイド等の微粒子にもエネルギー等分配の法則が成り立つと仮定すると $$ \frac{1}{2} m \langle v^2 \rangle = \frac{1}{2} k_B T $$ であるから、これを代入し、さらに$y(t) = \langle x(t)^2 \rangle$と記号を置き直すことで、 $$ m \frac{d^2 y}{d t^2} + \gamma \frac{d y}{d t} = 2 k_B T $$ を得る。

簡単な微分方程式なのですぐに手で解けるものの、このページの趣旨に鑑みて、SymPyで解かせてみると

from sympy import Function, Symbol,dsolve, Derivative, latex

t = Symbol('t')
m = Symbol('m',positive=True)
gamma = Symbol('gamma',positive=True)
kBT = Symbol('kBT',positive=True)

y = Function('y')

result = dsolve(m * Derivative(y(t), t, t) + gamma * Derivative(y(t),t) - 2*kBT, y(t))
print(result)
# print(latex(result))

Eq(y(t), C1 + C2*exp(-gamma*t/m) + 2*kBT*t/gamma)

を得る。初期条件$y(0)=0$に加え、十分に時間が経過すると指数関数の項は0と見なせるから、長時間的な振る舞いは $$ y(t) = \langle (x(t))^2 \rangle = \frac{2 k_B T}{\gamma} t $$ となる。

拡散定数$D$を、時間が$\tau$だけ経過する間の位置の平均自乗変位を用いて $$ D \equiv \frac{\langle x^2 \rangle}{2 \tau} $$ と定義すると、ブラウン運動の二乗変位の結果と合わせて $$ D = \frac{k_B T}{\gamma} $$ と表せる（ ここで、粘性抵抗にストークスの式$\gamma = 6 \pi \eta a$（$\eta$:粘性係数、$a$:粒子の半径）を用いたのがストークス・アインシュタインの式）。

速度分布と揺動力

一方、速度の変動を表す式(2)を $$ \frac{d v(t)}{dt} = -\frac{\gamma}{m} v(t) + \frac{1}{m} \xi $$ と書き直すと、これは速度空間のポテンシャル $$ U(v) = - \frac{\gamma}{2 m} v^2 $$ 中の運動と見なせるが、そのときの定常分布は、こちらのページにあるとおり（ただし、$x$を$v$と読み替えて） $$ p(v) \sim e^{-\frac{2}{R} U(v)} \tag{3} $$ となる。ここで$R$は揺動項（この場合は$\xi/m$）の自己相関で、 $$ R = \langle \xi(t)^2 \rangle / m^2 $$ に対応する。

熱平衡状態にある物理的な系を考えているので、微粒子のエネルギーを$E$とすると、その確率分布はボルツマン分布 $$ p(E) \sim e^{-\frac{E}{k_B T}} $$ に従う。 この場合、粒子の運動エネルギーのみを考慮すれば良いので $E = \frac{1}{2} m v^2$ であるから $$ p(E(v)) \sim e^{-\frac{m v^2}{2 k_B T}} \tag{4} $$ である。

式(3)と(4)は本来同じ現象を違う見方をしているのだから、両者を見比べることで、式(2)の揺動項の自己相関と物理量との関係 $$ \langle \xi(t)^2 \rangle = 2 \gamma k_B T $$ を得る。 このように、周囲の分子の熱運動によって微粒子に作用する実効的な揺動力は、温度（周囲の分子の運動エネルギー）に比例するだけでなく、 （熱平衡であるための要請から）摩擦の程度とも連動している。

　練習：ブラウン運動のシミュレーション

確率微分方程式(2)に基づいて、ブラウン運動のシミュレーションを行ってみなさい。

ここでは、温度 $T=300{\textrm K}$、粒子の半径$a=1 \mu {\textrm m}$、水の粘性係数 $\eta = 0.001 \textrm{Pa} \cdot {\textrm s}$ とする。