Julia言語によるLangevin方程式の数値計算

横浜市立大学大学院
生命ナノシステム科学研究科 前期博士課程 物質システム科学専攻
量子物理化学研究室

大野 周平

2021年10月23日(土) @ 第8回 ぶつりがく徒のつどい

はじめに

本講演ではLangevin方程式の歴史背景と数値解法について概説する.

計算にはJulia言語, 描写にはPlots.jlを用いた. 全てのソースコードは講義ノートにまとめてある.
講義ノートはIJulia.jlによるJupyter Notebookでまとめており, 計算や描写はノート上で行った.
このスライドはMarp for VS Codeを用いて作成した.

資料

資料は全て https://github.com/ohno/butsudoi2021 にアップロードしてあります.

質問・指摘等は https://github.com/ohno/butsudoi2021/discussions にお願いします.

Julia言語

Julia言語とは2018年にバージョン1.0がリリースされた Pythonのように書けて, Cのように動く 次世代のプログラミング言語である. JuliaはPythonの書きやすさとCの速さを両立しており, さらにFortranに比類する強力な配列操作と行列演算が可能である. 単体での性能もさることながら, Python, C, Fortranの豊富なライブラリを呼び出すなど, 他言語との糊の役割を果たすグルー言語としても優秀である. 以上に説明した通り, Julia言語は速さ, 書きやすさ, 連携性など, 計算科学における標準語として活躍できるだけの資質を備えている.

また, Julia言語の疑似乱数にはdSFMTが採用されており, 特別な準備をせずに良質な疑似乱数を用いてモンテカルロ計算を行うことができる. 乱数を使う分野だと最強！

個人的なJulia歴

年 エピソード
2014 第4回世界トップレベル研究拠点プログラム(WPI)合同シンポジウムにて, どこかの
おじいさん（たぶんすごい研究者）から「最近はJuliaという速い言語があるらしい」
と聞いていた.（開発スタート2009年, 初回リリース2012年なのでホットな時期）
2015～ 当時はExcel VBAで数値計算をしていた. 以降, Python, Octave, R, C, JavaScriptなどをかじるが, 圧倒的な書きやすさと速さを誇るFortranに魅了される.
2018 Ver.1.0リリースをきっかけに使い始める. 印象は「使いやすいFortran」.
2020 Julia言語で量子モンテカルロ法のプログラムを開発. スパコン上で運用.
2021 研究室のFORTRANをJuliaから呼び出すAPIを開発し, スパコン上で運用中.

Julia in Physics 2021 Online

2021年9月3日(金) 13:30～

Juliaの入門講座や物理学における使用例の紹介があった. スライドや講演動画がアップロードされているので, Juliaに入門したい方におすすめ.

サンプルコード

# `plot()`を使用するためには事前に`using Plots`を宣言しておく必要がある.
# using Plots

# 型の宣言は不要である. zeros()は値が全て0の配列を生成する
step = 20000
X = zeros(step)
Y = zeros(step)

# For文は以下のように書く,
# randn()は平均0分散1の正規分布に従う疑似乱数.
for i in 2:step
    X[i] = X[i-1] + randn()
    Y[i] = Y[i-1] + randn()
end

# 描写は`plot(X,Y)`のように2軸の配列を渡す.
plot(X, Y, label="", xlabel="x", ylabel="y")

# 保存は`savefig(ファイル名)`
savefig("sample.svg")
plot!()

Plots.jlの入門はこちらのノートを参照.

目標

拡散方程式のMSDではなく, Langevin方程式のMSDを数値的に再現すること！

名称 拡散方程式 Langevin方程式
方程式 ∂P(x,t)∂t=D∇2P(x,t)\frac{\partial P(\pmb{x},t)}{\partial t}=D\nabla^2 P(\pmb{x},t) mdV(t)dt=−γV(t)+R(t)m\frac{\mathrm{d}\pmb{V}(t)}{\mathrm{d}t} = -\gamma \pmb{V}(t) + \pmb{R}(t)
MSD ⟨∣X(t)∣2⟩=2nDt\langle|\pmb{X}(t)|^2\rangle = 2nDt ⟨∣X(t)∣2⟩=2nDumγ3(γmt+e−γmt−1)\langle|\pmb{X}(t)|^2\rangle = 2n\frac{D_u m}{\gamma^3} \left(\frac{\gamma}{m} t+e^{-\frac{\gamma}{m} t}-1 \right)

初めてこのグラフを見た人の反応

  • 異常拡散では？
    → 厳密なBrown運動のMSDです
  • こんなの見たことない
  • 出典は？
  • 合ってるの？
    → Kitahara(1997)の式(4.59)
    　Zwanzig(2001)の式(1.34)
    　このノートの3.4節を読んで！
  • 積分が難しい！
    → Juliaで検証してみよう！

第１部

Brown(1828), Einstein(1905), Perrin(1908), Langevin(1908)の一連の研究について簡単に紹介する.

Brown(1828)

  • Brownの功績はBrown運動を物理学の世界へと持ち込んだこと.
  • Brownは液体中の微粒子のランダムな運動が生命現象である可能性を否定した.
  • Brown運動はBrownが最初に発見したわけではないことは本人も明言している.

never considered by me as wholly original　―　Brown(1828) p.171


Wikipedia『ブラウン運動にまつわる誤解』でも取り上げられているように,
Brownの観察した微粒子はよく花粉と間違えられる. 例えば,

Brown運動は水に浮かぶ花粉の不規則運動
―　L.カラザス, S.E.シュレーブ『ブラウン運動と確率積分』(2012) 第2章 2.1序論

Einsteinまでの80年間

アインシュタインの理論が成功したのは, ブラウン粒子の速度に関する部分を表向き回避した点にあると言ってもよい.
―　米沢登美子『物理学 One Point - 27 ブラウン運動』(共立出版, 1986) p.50

Einstein(1905)

Einsteinの三大業績の一つ. 「奇跡の年」の2本目の論文. この論文では粒子の平均二乗変位からアボガドロ定数NNを算出する理論を発表した. 最後の式

N=tλx2⋅RT3πkPN = \frac{t}{{\lambda_x}^2} \cdot \frac{RT}{3\pi kP}

から, 根平均二乗変位λx{\lambda_x}を実験的に求めれば良いことがわかる. Einsteinは導出の過程で拡散方程式とその根平均二乗変位（RMSD）

∂f∂t=D∂2f∂x2\frac{\partial f}{\partial t}=D\frac{\partial^2 f}{\partial x^2}

λx=2Dt\lambda_x = \sqrt{2Dt}

を用いた. 詳細は 米沢登美子『ブラウン運動』(共立出版, 1986) p.56を参照.

Perrin(1908)

最終行にNA=6.7×1023N_\mathrm{A}=6.7\times10^{23}が得られたと書いてある. 当時はまだ原子や分子の実在性に疑問を持つ学者も多く, 革新的な結果だった.

ぺランは, 一連の実験により,（最初の結果を公表したのは1908年の論文）それまで誰にも達成できなかった正確さで, アインシュタインの予測のほとんどすべてを裏付けた.
―　ジョン・スタチェル編, 青木薫訳『アインシュタイン論文選「奇跡の年」の5論文』(ちくま学芸文庫, 2011) p.201

Perrinは妥当か？

2021年現在の最新のアボガドロ定数
は2018年CODATA推奨値

NA=6.022 140 76×1023mol−1N_\mathrm{A} = 6.022~140~76 \times 10^{23} \mathrm{mol}^{-1}

Becker(2001)によると, 当時としてはおかしくないことがわかる. Perrinが1926年にノーベル物理学賞を授賞する頃までは分子運動による測定方法が主流だった
→ 拡散係数を過小評価
→ アボガドロ定数を過大評価？

Langevin(1908)

今日の本題. EinsteinとSmoluchowskiの論文を元に, 粒子が粘性抵抗とランダムな力(揺動力)を受けるときの運動方程式としてLangevin方程式

md2xdt2=−6πμadxdt+Xm\frac{\mathrm{d}^2 x}{\mathrm{d}t^2} = -6\pi\mu a \frac{\mathrm{d} x}{\mathrm{d}t} + X

を提案した. また, Langevin方程式から平均二乗変位を導出し, Einsteinの論文と同じ結果を得た.

Δx2ˉ=x2ˉ−x02ˉ=RTN13πμaτ\bar{\Delta_x^2} = \bar{x^2}-\bar{x^2_0} = \frac{RT}{N}\frac{1}{3\pi\mu a}\tau

まとめ

  • Brown(1828)はBrown運動が生命現象である可能性を否定.
  • Einstein(1905)は粒子の平均二乗変位からアボガドロ定数を算出する理論を提唱.
  • Perrin(1908)でNA=6.7×1023N_A=6.7\times10^{23}が得られ, 原子・分子の存在が実証された.
  • Langevin(1908) はLangevin方程式を提案し, Einsteinと同じ結果を導いた.

拡散方程式とLangevin方程式のMSDは異なる, という認識に基づけば, Einstein(1905)とLangevin(1908)が同じ結論にたどり着いたことは非常に面白い！

第２部

確率過程の概念, Eular-Maruyama法による種々の確率微分方程式の数値解法について, Julia言語での実装例を交えて解説する.

確率過程

確率過程（stochastic process）とは, 時間に依存して変化する確率変数（stochastic variable）のことである. 時間がt=1,2,3,4,⋯t=1,2,3,4,\cdotsのように離散的な値をとる確率過程を離散確率過程（discrete stochastic process）といい, 連続的な値をとる確率過程を連続確率過程（containuous stochastic process）という.

いい加減な定義ではあるものの, 確率変数がX:Ω→RnX:\Omega\rightarrow\mathbb{R^n}であり, 確率過程がX:Ω×T→RnX:\Omega\times T\rightarrow\mathbb{R^n}であるから, 概ね正しいと思われる.

ランダムウォークモデル

離散確率過程の例：ランダムウォークモデル

コインを投げて表ならRi=+1R_i=+1, 裏ならRi=−1R_i=-1とする. コインをii回目に投げたときの合計得点をXi=∑j=1iRjX_i=\sum_{j=1}^{i}R_jとする.

i=1i=1のとき, Xi=+1,−1X_i=+1,-1のどちらか
i=2i=2のとき, Xi=+2,0,−2X_i=+2,0,-2のどれか
i=3i=3のとき, Xi=+3,+1,0,−1,−3X_i=+3,+1,0,-1,-3

つまり, XiX_iの確率分布が時間iiに依存する.
なお, RiR_iの確率分布は時間に依らない.

Wiener過程

連続確率過程の例：Wiener過程W(t)W(t)

⟨R(t)⟩=0⟨R(t1)R(t2)⟩=2Duδ(t1−t2)⟨W(t)⟩=0⟨W(t1)W(t2)⟩=2Dumin(t1,t2)W(t)∼N(0,t)\begin{aligned} \langle R(t) \rangle &= 0\\ \langle R(t_1)R(t_2) \rangle &= 2D_u \delta(t_1-t_2)\\ \\ \langle W(t) \rangle &= 0\\ \langle W(t_1)W(t_2) \rangle &= 2D_u {\rm min}(t_1,t_2)\\ W(t) &\sim N(0,t) \end{aligned}

白色雑音R(t)R(t)はWiener過程W(t)W(t)の時間微分のようなものであり, ちょうどコインの裏表RiR_iと合計得点数XiX_iの関係によく似たものである.

白色雑音を取り巻く問題

結論から言うと, このような白色雑音は, 数学的には確率過程のカテゴリーに入れて議論することができない

白色雑音を数学的に厳密に取り扱うには, 確率超過程と呼ばれる概念を導入しなければならない
―　兼清泰明『確率微分方程式とその応用』(森北出版, 2017) p.100

Langevin方程式を扱う物理学の文脈では, 白色雑音は揺動力として導入され, その性質を出発点とした相関関数の計算によって議論が進められる. 一方, 数学では白色雑音の扱いを避けて, Wiener過程を議論の出発点に採用することが多いようである.

このノートでも, 白色雑音の性質は不要なので使わない. また, Wiener過程についてもW(t)−W(s)∼N(0,2Du(t−s))W(t)-W(s) \sim N(0,2D_u(t-s))という性質を使うのみである.

Wiener過程の最も重要な性質

Wiener過程W(t)W(t)の最も重要な性質は

W(t)−W(s)∼N(0,2Du(t−s)) W(t)-W(s) \sim N(0,2D_u(t-s))

である. 兼清(2017) p.103ではこの性質をWiener過程W(t)W(t)の定義に含めるが, 三井(2004) p.139では定理として導出している.

いずれにせよ, この性質によってΔt\Delta tあたりの変化量ΔWi\Delta W_iが正規分布N(0,2DuΔt)N(0,2D_u\Delta t)に従う疑似乱数として計算できる.

このΔWi\Delta W_iが今日の計算のカギ！

確率微分方程式

微分方程式の解は微分可能でなければならない. Brown運動の軌道は明らかに微分可能ではないから, Brown運動の運動方程式(微分方程式)を考えるのはおかしい.

そのようなわけで, Langevin方程式は通常の微分方程式としては解釈できない. ではどのように解釈すればよいかというと, 確率積分を含む方程式の別表記と考えるのである. というわけで, SDEは微分方程式というよりはむしろ積分方程式である.

確率微分方程式には確率積分としてIto積分を採用する立場とStratonovich積分を採用する立場の2つの立場があるが, このノートでは全てIto積分を採用する.

Ito方程式

(3.2.3)の第2の積分をIto積分と解釈したとき, これらの方程式をItoのLangevin方程式あるいはItoの方程式という.
　―　堀淳一『ランジュバン方程式』(岩波書店, 2015) p.61

dX(t)dt=f(X(t),t)+g(X(t),t)R(t)(3.2.1)dX(t)=f(X(t),t)dt+g(X(t),t)dW(t)(3.2.2)X(t)−X(t0)=∫t0t1f(X(s),s)ds+∫t0t1g(X(s),s)dW(s)(3.2.3)\begin{aligned} \frac{{\rm d}X(t)}{{\rm d}t} = f(X(t),t) + g(X(t),t)R(t) &&(3.2.1)\\ \mathrm{d}X(t) = f(X(t),t)\mathrm{d}t + g(X(t),t)\mathrm{d}W(t) &&(3.2.2)\\ X(t)-X(t_0) = \int_{t_0}^{t_1} f(X(s),s){\rm d}s + \int_{t_0}^{t_1} g(X(s),s){\rm d}W(s) &&(3.2.3) \end{aligned}

(3.2.3)の右辺の第2項の積分は, Ito積分を採用することで決まった値を持つようになり,(3.2.3)が意味を持つようになる. (3.2.1)と(3.2.2)は形式的にdW(t)=R(t)dt\mathrm{d}W(t)=R(t)\mathrm{d}tを認めたときの(3.2.3)の別表記であり, 同じ意味を持つものとする. これによって, (3.2.1)のような通常の微分方程式と同じ表記が使えるようになる.

Fokker-Planck方程式

確率論と解析学の間には実り多い相互関係がある. その研究は少なくともKolmogorov(1931)にさかのぼる. (中略) 楕円型, 放物型偏微分方程式の多くの問題の解は確率汎関数の期待値として表現することができる.
―　L.カラザス, S.E.シュレーブ『ブラウン運動と確率積分』(2012) 第4章 4.1序論

このノートでは, 上記のEular-Maruyama法による数値解の妥当性をFokker-Planck方程式によって評価する. Ito方程式

dX(t)dt=f(X(t),t)+g(X(t),t)R(t) \frac{{\rm d}X(t)}{{\rm d}t} = f(X(t),t) + g(X(t),t)R(t)

の解X(t)X(t)の確率密度関数P(x,t)P(x,t)が従う偏微分方程式としてFokker-Planck方程式

∂P(x,t)∂t=∂2[Du(g(x,t))2P(x,t)]∂x2−∂[f(x,t)P(x,t)]∂x \frac{\partial P(x,t)}{\partial t} = \frac{\partial^2 [D_u (g(x,t))^2 P(x,t)]}{\partial x^2}-\frac{\partial [f(x,t)P(x,t)]}{\partial x}

が知られている.

以下に具体例をまとめる. 今回はD=DuD=D_uとしてよい.

条件 Ito方程式 Fokker-Planck方程式
f(X(t),t)=0,g(X(t),t)=1f(X(t),t)=0, g(X(t),t)=1 dX(t)dt=R(t){\frac{{\rm d}X(t)}{{\rm d}t}=R(t)}
Wiener過程		 ∂P(x,t)∂t=D∂2P(x,t)∂x2\frac{\partial P(x,t)}{\partial t}=D\frac{\partial^2 P(x,t)}{\partial x^2}
拡散方程式
f(X(t),t)=c,g(X(t),t)=1f(X(t),t)=c, g(X(t),t)=1 dX(t)dt=c+R(t){\frac{{\rm d}X(t)}{{\rm d}t}=c + R(t)}
ドリフト項のあるWiener過程		 ∂P(x,t)∂t=D∂2P(x,t)∂x2−c∂P(x,t)∂x\frac{\partial P(x,t)}{\partial t}=D\frac{\partial^2 P(x,t)}{\partial x^2}-c\frac{\partial P(x,t)}{\partial x}
移流拡散方程式
f(X(t),t)=−γX(t),g(X(t),t)=1f(X(t),t)=-\gamma X(t), g(X(t),t)=1 dX(t)dt=−γX(t)+R(t){\frac{{\rm d}X(t)}{{\rm d}t}=-\gamma X(t) + R(t)}
Ornstein–Uhlenbeck過程		 ∂P(x,t)∂t=D∂2P(x,t)∂x2−∂(−γXP(x,t))∂x\frac{\partial P(x,t)}{\partial t}=D\frac{\partial^2 P(x,t)}{\partial x^2}-\frac{\partial (-\gamma XP(x,t))}{\partial x}
名前なし

初期条件P(x,0)=δ(x−x0)P(x,0) = \delta(x-x_0), 境界条件P(±∞,t)=0P(\pm\infty,t)=0のときの解はそれぞれ以下の通りである.

P(x,t)=14πDtexp⁡(−(x−x0)24Dt)P(x,t)=14πDtexp⁡(−(x−x0−ct)24Dt)P(x,t)=γ2πD(1−e−2γt)exp⁡(−γ(x−x0e−γt)22D(1−e−2γt))\begin{aligned} P(x,t) = \frac{1}{\sqrt{4\pi Dt}} \exp \left( -\frac{(x-x_0)^2}{4Dt} \right) \\ P(x,t) = \frac{1}{\sqrt{4\pi Dt}} \exp \left( -\frac{(x-x_0-ct)^2}{4Dt} \right) \\ P(x,t) = \sqrt{\frac{\gamma}{2\pi D(1-e^{-2\gamma t})}} \exp \left( -\frac{\gamma(x-x_0 e^{-\gamma t})^2}{2D(1-e^{-2\gamma t})} \right) \end{aligned}

Eular-Maruyama法

Eular-Maruyama法は確率初期値問題

dX(t)dt=f(X(t),t)+g(X(t),t)R(t),t∈[t0,T],X(t0)=X0 \frac{{\rm d}X(t)}{{\rm d}t} = f(X(t),t) + g(X(t),t)R(t), t\in[t_0,T], X(t_0)=X_0

を次のように離散化する方法である.

Xi+1≃Xi+f(Xi,t)Δt+g(Xi,t)ΔWi X_{i+1} \simeq X_i + f(X_i, t)\Delta t + g(X_i, t)\Delta W_i

このとき, Δt\Delta tは時間刻み幅で, 0.1とか0.01などを選ぶ. ΔWi\Delta W_iは正規分布N(0,2DuΔt)N(0,2D_u\Delta t)に従う疑似乱数である. Eular-Maruyama法はEular法の一般化であり, g(X(t),t)=0g(X(t),t)=0の場合はEular法と一致する. Juliaでの実装：

for i in 1:step
    t += dt
    X += f.(X,t)*dt + g.(X,t).*sqrt(2*Du*dt).*randn(N)
end

数値計算：Wiener過程

f(X(t),t)=0,g(X(t),t)=1f(X(t),t)=0, g(X(t),t)=1のとき
dX(t)dt=R(t)dX(t)=dW(t)X(t)−X(t0)=∫t0t1dW(s)\frac{{\rm d}X(t)}{{\rm d}t} = R(t)\\{\rm d}X(t) = {\rm d}W(t)\\X(t)-X(t_0) = \int_{t_0}^{t_1} {\rm d}W(s)
Eular-Maruyama法による離散化
Xi+1≃Xi+ΔWiX_{i+1} \simeq X_i + \Delta W_i
Fokker-Planck方程式（拡散方程式）
∂P(x,t)∂t=D∂2P(x,t)∂x2\frac{\partial P(x,t)}{\partial t}=D\frac{\partial^2 P(x,t)}{\partial x^2}
対応する初期条件
X(0)=X0P(x,0)=δ(x−x0)X(0) = X_0\\P(x,0) = \delta(x-x_0)
境界条件P(±∞,t)=0P(\pm\infty,t)=0のときの解
P(x,t)=14πDtexp⁡(−(x−x0)24Dt)P(x,t) = \frac{1}{\sqrt{4\pi Dt}} \exp \left( -\frac{(x-x_0)^2}{4Dt} \right)

数値計算：ドリフト項のあるWiener過程

f(X(t),t)=c,g(X(t),t)=1f(X(t),t)=c, g(X(t),t)=1のとき
dX(t)dt=c+R(t)dX(t)=cdt+dW(t)X(t)−X(t0)=∫t0t1cds+∫t0t1dW(s)\frac{{\rm d}X(t)}{{\rm d}t} = c + R(t)\\{\rm d}X(t) = c{\rm d}t + {\rm d}W(t)\\X(t)-X(t_0) = \int_{t_0}^{t_1} c{\rm d}s + \int_{t_0}^{t_1} {\rm d}W(s)
Eular-Maruyama法による離散化
Xi+1≃Xi+cΔt+ΔWiX_{i+1} \simeq X_i + c\Delta t + \Delta W_i
Fokker-Planck方程式（移流拡散方程式）
∂P(x,t)∂t=D∂2P(x,t)∂x2−c∂P(x,t)∂x\frac{\partial P(x,t)}{\partial t} = D\frac{\partial^2 P(x,t)}{\partial x^2} - c \frac{\partial P(x,t)}{\partial x}
対応する初期条件
X(0)=X0P(x,0)=δ(x−x0)X(0) = X_0\\P(x,0) = \delta(x-x_0)
境界条件P(±∞,t)=0P(\pm\infty,t)=0のときの解
P(x,t)=14πDtexp⁡(−(x−x0−ct)24Dt)P(x,t) = \frac{1}{\sqrt{4\pi Dt}} \exp \left( -\frac{(x-x_0-ct)^2}{4Dt} \right)

数値計算：Ornstein–Uhlenbeck過程

f(X(t),t)=−γX(t),g(X(t),t)=1f(X(t),t)=-\gamma X(t), g(X(t),t)=1のとき
dX(t)dt=−γX(t)+R(t)dX(t)=−γX(t)dt+dW(t)X(t)−X(t0)=∫t0t1−kX(s)ds+∫t0t1dW(s)\frac{{\rm d}X(t)}{{\rm d}t} = -\gamma X(t) + R(t)\\{\rm d}X(t) = -\gamma X(t){\rm d}t + {\rm d}W(t)\\X(t)-X(t_0) = \int_{t_0}^{t_1} -kX(s){\rm d}s + \int_{t_0}^{t_1} {\rm d}W(s)
Eular-Maruyama法による離散化
Xi+1≃Xi+cΔt+ΔWiX_{i+1} \simeq X_i + c\Delta t + \Delta W_i
Fokker-Planck方程式
∂P(x,t)∂t=D∂2P(x,t)∂x2−∂(−γXP(x,t))∂x\frac{\partial P(x,t)}{\partial t} = D\frac{\partial^2 P(x,t)}{\partial x^2} - \frac{\partial (-\gamma XP(x,t))}{\partial x}
対応する初期条件
X(0)=X0P(x,0)=δ(x−x0)X(0) = X_0\\P(x,0) = \delta(x-x_0)
境界条件P(±∞,t)=0P(\pm\infty,t)=0のときの解
P(x,t)=γ2πD(1−e−2γt)exp⁡(−γ(x−x0e−γt)22D(1−e−2γt))P(x,t) = \sqrt{\frac{\gamma}{2\pi D(1-e^{-2\gamma t})}} \exp \left( -\frac{\gamma(x-x_0 e^{-\gamma t})^2}{2D(1-e^{-2\gamma t})} \right)

まとめ

  • Wiener過程の増分ΔWi\Delta W_iの計算法を理解した！
  • 確率微分方程式をEular-Maruyama法で計算できた！
  • Fokker-Plack方程式との対応関係を理解した！

第３部

第2部の内容を元に, Langevin方程式の数値計算を行う. 北原和夫『非平衡系の統計力学』(1997)から平均二乗変位の式を援用し, 計算結果の評価を行う.

Langevin方程式

Langevin(1908)で与えられたLangevin方程式は次式である.

md2xdt2=−6πμadxdt+Xm\frac{\mathrm{d}^2 x}{\mathrm{d}t^2} = -6\pi\mu a \frac{\mathrm{d} x}{\mathrm{d}t} + X

これは, Stokesの式γ=6πμa\gamma=6\pi\mu aと速度V(t):=dX(t)/dtV(t):=\mathrm{d}X(t)/\mathrm{d}tの定義を使って次のように書き直せる.

dX(t)dt=V(t)mdV(t)dt=−γV(t)+R(t)\begin{aligned} \frac{\mathrm{d} X(t)}{\mathrm{d}t} &= V(t) \\ m\frac{\mathrm{d} V(t)}{\mathrm{d}t} &= -\gamma V(t) + R(t) \end{aligned}

速度V(t)V(t)の方程式はOrnstein–Uhlenbeck過程として解釈できる.

Langevin方程式の離散化

速度V(t)V(t)の方程式にはEular-Maruyama法を, 変位X(t)X(t)の方程式には通常のEular法を使う.

Xi+1≃Xi+ViΔtVi+1≃Vi−γmΔt+1mΔWi\begin{aligned} X_{i+1} \simeq X_i + V_i\Delta t\\ V_{i+1} \simeq V_i - \frac{\gamma}{m}\Delta t + \frac{1}{m}\Delta W_i \end{aligned}

このようにしてBrown運動の軌道を計算できる. しかし, 何を以てこの計算の妥当性を評価するかという問題がある. Ornstein–Uhlenbeck過程のP(v,t)P(v,t)は既に明らかになっているが, P(x,t)P(x,t)を求めるのは難しそうである. その代わりに, 平均二乗変位を用いて計算の妥当性を評価する.

平均二乗変位

北原和夫『岩波基礎物理シリーズ 8 非平衡系の統計力学』(岩波書店, 1997)はp.90の式(4.38)

mdudt=−mγu+R(t)         (4.38)m\frac{\mathrm{d}\boldsymbol{u}}{\mathrm{d}t} = -m\gamma\boldsymbol{u} + \boldsymbol{R}(t)~~~~~~~~~(4.38)

を出発点とし, p.94の式(4.59)の平均二乗変位を導出している.

⟨∣x(t)∣2⟩=6kBTmγ2(γt+e−γt−1)         (4.59)\langle|\boldsymbol{x}(t)|^2\rangle = \frac{6k_\mathrm{B}T}{m\gamma^2}(\gamma t+e^{-\gamma t}-1)~~~~~~~~~(4.59)

この式(4.59)を援用するが, 出発点となる方程式が異なるので, そのままでは検証には使えない.

北原の罠

では粘性抵抗を−mγu-m\gamma\boldsymbol{u}としており, 質量に比例するものとしている. 出発点が異なるため, 平均二乗変位の式(4.59)もそのままではLangevin(1908)の検証には使えない.

γ′=mγ\gamma'=m\gamma, Du=mγkBT=γ′kBTD_u=m\gamma k_\mathrm{B}T=\gamma' k_\mathrm{B}Tとし, さらにnn次元の場合(6→2n6\rightarrow2n)を考えると

⟨∣X(t)∣2⟩=2nDumγ′3(γ′mt+e−γ′mt−1) \langle|\boldsymbol{X}(t)|^2\rangle = 2n\frac{D_u m}{\gamma'^3} \left(\frac{\gamma'}{m} t+e^{-\frac{\gamma'}{m} t}-1 \right)

となる. これを用いて, Langevin方程式の数値計算の結果を評価する.

また, 熱平衡状態を初期値として仮定していることにも注意しなければならない.

数値計算：Langevin方程式

上：初速0の場合
下：初期条件を熱平衡状態にした場合

方程式
dX(t)dt=V(t)mdV(t)dt=−γV(t)+R(t)\frac{\mathrm{d} X(t)}{\mathrm{d}t}= V(t)\\m\frac{\mathrm{d} V(t)}{\mathrm{d}t}= -\gamma V(t) + R(t)
離散化
Xi+1≃Xi+ViΔtVi+1≃Vi−γmΔt+1mΔWiX_{i+1} \simeq X_i + V_i\Delta t\\V_{i+1} \simeq V_i - \frac{\gamma}{m}\Delta t + \frac{1}{m}\Delta W_i
MSD
⟨∣X(t)∣2⟩=2nDumγ′3(γ′mt+e−γ′mt−1)\langle|\boldsymbol{X}(t)|^2\rangle = 2n\frac{D_u m}{\gamma'^3} \left(\frac{\gamma'}{m} t+e^{-\frac{\gamma'}{m} t}-1 \right)

拡散係数の過小評価

Langevin方程式のMSD

⟨∣X(t)∣2⟩=2nDumγ3(γmt+e−γmt−1) \langle|\boldsymbol{X}(t)|^2\rangle = 2n\frac{D_u m}{\gamma^3} \left(\frac{\gamma}{m} t+e^{-\frac{\gamma}{m} t}-1 \right)

は長時間極限t→∞t\rightarrow\inftyではe−γt→0e^{-\gamma t}\rightarrow0となり, 定数項も相対的には無視できることから, 第1項が支配的となり,

⟨∣X(t)∣2⟩≃2nDumγ3γmt \langle|\boldsymbol{X}(t)|^2\rangle \simeq 2n\frac{D_u m}{\gamma^3} \frac{\gamma}{m} t

と考えることができる. さらに, これを拡散方程式のMSD⟨∣X(t)∣2⟩=2nDt\langle|\boldsymbol{X}(t)|^2\rangle = 2nDtと見比べ, Du=γkBTD_u=\gamma k_\mathrm{B}Tを代入すると, よく目にするEinsteinの関係式（揺動散逸定理）

D=Duγ2=kBTγ D = \frac{D_u}{\gamma^2} = \frac{k_\mathrm{B}T}{\gamma}

が得られる. すなわち, どのような粒子であっても十分に長い時間間隔で観測すれば, 拡散係数Du/γ2D_u/\gamma^2の拡散方程式に従う. 拡散係数を過小評価することが自然である.

揺動力の強さDuD_uと拡散係数DD

後に説明するが, 拡散方程式(f(X(t),t)=0,g(X(t),t)=1f(X(t),t)=0, g(X(t),t)=1のIto方程式)とLangevin方程式の導く平均二乗変位⟨∣X(t)∣2⟩\langle|X(t)|^2\rangleは異なるものである. 拡散係数を求める際によく目にする⟨∣X(t)∣2⟩=2nDt\langle|X(t)|^2\rangle = 2nDtとLangevin方程式から導かれる⟨∣X(t)∣2⟩=2nDumγ3(γmt+e−γmt−1)\langle|\boldsymbol{X}(t)|^2\rangle = 2n\frac{D_u m}{\gamma^3} \left(\frac{\gamma}{m} t+e^{-\frac{\gamma}{m} t}-1 \right)は, t→∞t\rightarrow\inftyやγm≃0\frac{\gamma}{m}\simeq0などの特定の条件下でのみEinsteinの関係式

D=Duγ2=kBTγ D = \frac{D_u}{\gamma^2} = \frac{k_\mathrm{B}T}{\gamma}

で結ばれる. Fokker-Planck方程式の説明だけを見ると揺動力の強さDuD_uと拡散係数DDは一致するように見えるが, Langrvin方程式を見据えると, 揺動力の強さDuD_uと拡散係数DDは明確に区別することが教育的であると思われる.

まとめ

  • Langevin方程式のMSDを数値的に検証できた！
  • Langevin方程式の数値計算には落とし穴が多かった！
  • MD計算では拡散係数を過小評価することが自然だとわかった！

最後に

慣性を含んだBrown運動のシミュレーション（2階微分方程式としてのLangevin方程式の計算）はあまり多くないと思われる. これに加えて, 相互作用のあるBrown運動の計算や, 場の中でのBrown運動の計算には新規性のあるテーマがあると思われるので, そういったBrown運動の計算やノイズの関わる数理モデルの構築に興味があればぜひ共同研究をしたい.

参考文献：第１部

参考文献：第２部

参考文献：第３部

Ito方程式(3.2.1)～(3.2.3)をBrown運動の運動方程式(Langevin方程式)として解釈しようとすると, $X(t)$は速度だと思わなくてはいけないので$V(t)$の記号に変えた方がよい. また, $f(X(t),t)=-\gamma X(t), g(X(t),t)=1$の時が狭義でのLangevin方程式なので, かなり一般化されている. そのため, 「一般化Langevin方程式」とでも呼びたいところだが, これはまた別の方程式を指す名称なので紛らわしい. (3.2.1)～(3.2.3)をItoのLangevin方程式と呼ぶのは誤解を招くであろうから, このノートでは単にIto方程式と呼ぶ.