Langevin法と確率微分方程式 / Langevin Equation

Langevin方程式の数値解法、Fokker-Planck方程式、伊藤/Stratonovich積分について

Transcript

1. 1 Langevin法と確率微分方程式 慶應義塾大学理工学部物理情報工学科 渡辺 2021/07/02
2. 2 確率微分方程式の厳密な取り扱いはかなり面倒です ここで紹介する確率微分方程式の解釈は、「とりあえずそう思っておく」 程度と認識してください
3. 3 分子動力学法(Molecular Dynamics Method, MD)とは 狭義にはハミルトンの運動方程式を数値積分する手法 ሶ 𝑝 = − 𝜕𝐻 𝜕𝑞 ሶ 𝑞 = 𝜕𝐻 𝜕𝑝 𝑑𝐻 𝑑𝑡 = 𝜕𝐻 𝜕𝑝 ሶ 𝑝 + 𝜕𝐻 𝜕𝑞 ሶ 𝑞 = 0 より、この時間発展ではエネルギーが保存される (NVEアンサンブル) 我々が興味があるのは、物理量のエネルギー依存性ではなく温度依存性 温度を制御したい
4. 4 ハミルトンの運動方程式に、温度制御項を付け加える ሶ 𝑝 = − 𝜕𝐻 𝜕𝑞 − 𝛾 𝜕𝐻 𝜕𝑝 + 2𝐷 ෠ 𝑅 ሶ 𝑞 = 𝜕𝐻 𝜕𝑝 散逸項 揺動項 散逸項：エネルギーを減らそうとする 揺動項：エネルギーを注入する 散逸と揺動が釣り合ったところで温度が一定となる → Langevin法 温度制御項
5. 5 ሶ 𝑝 = − 𝜕𝐻 𝜕𝑞 − 𝛾 𝜕𝐻 𝜕𝑝 + 2𝐷 ෠ 𝑅 ሶ 𝑞 = 𝜕𝐻 𝜕𝑝 揺動項 揺動項のRは時間に依存する確率的な変数であり、以下を満たす ෠ 𝑅(𝑡) = 0 ෠ 𝑅(𝑡1 ) ෠ 𝑅(𝑡2 ) = 𝛿(𝑡1- 𝑡2) 平均ゼロ 無相関(白色雑音) このような確率変数を含む微分方程式を確率微分方程式と呼ぶ 特に、ブラウン運動を記述する確率微分方程式をLangevin方程式と呼ぶ
6. 6 ሶ 𝑝 = − 𝜕𝐻 𝜕𝑞 − 𝛾 𝜕𝐻 𝜕𝑝 + 2𝐷 ෠ 𝑅 ሶ 𝑞 = 𝜕𝐻 𝜕𝑝 上記の時間発展に従う系の定常分布はカノニカル分布となる 𝑓(𝑝, 𝑞)~exp −𝛽𝐻 摩擦係数と拡散係数の比が温度を与える(揺動散逸定理) 𝛽 = 1 𝑘𝑇 = 𝛾/𝐷 (NVTアンサンブル)
7. 7 ሶ 𝑝 = − 𝜕𝐻 𝜕𝑞 − 𝛾 𝜕𝐻 𝜕𝑝 + 2𝐷 ෠ 𝑅 ሶ 𝑞 = 𝜕𝐻 𝜕𝑝 Langevin法は、温度制御項のついた運動方程式を 数値積分する手法 揺動項 Rは確率変数 • 確率変数を積分するにはどうすればよいか？ • 運動方程式と分布関数の関係はどうなっているか？
8. 8 積分区間を適当に分割し、以下の和をとる(Riemann和) න 𝑥𝑠 𝑥𝑒 𝑓 𝑥 𝑑𝑥 関数f(x)の定積分 を求めたい O x f(x) O x f(x) ≒
9. 9 分割の最大幅をδとする 任意の正数εに対してあるδが存在し、 ෍ 𝑖 𝑛−1 𝑓 𝑥𝑖 𝑥𝑖+1 − 𝑥𝑖 − 𝐴 < 𝜀 を満たすAを定義できる場合、f(x)はRiemann可積分であり、 Aをf(x)の定積分、もしくはRiemann積分と呼ぶ 𝛿 = max(𝑥𝑖+1 − 𝑥𝑖 ) න 𝑥𝑠 𝑥𝑒 𝑓 𝑥 𝑑𝑡 ≡ 𝐴 簡単にいえば、短冊の和の極限がRiemann積分 我々が通常使う「積分」とはRiemann積分のこと
10. 10 普通に積分できる場合、普通に数値積分できる ሶ 𝑥 = 𝑓(𝑥) を数値積分したい 𝑥 𝑡 + ℎ = 𝑥 𝑡 + ሶ 𝑥 𝑡 ℎ + 𝑂(ℎ2) 数値積分とは𝑥 𝑡 と ሶ 𝑥 𝑡 から 𝑥 𝑡 + ℎ を求めること න 𝑡 𝑡+ℎ ሶ 𝑥𝑑𝑡 = 𝑥 𝑡 + ℎ − 𝑥(𝑡) ←これは厳密に成り立つが 一般に積分は厳密に求められない Taylor展開して一次までとる 無視する 𝑥 𝑡 と ሶ 𝑥 𝑡 から 𝑥 𝑡 + ℎ が近似的に求まった(一次のEuler法) Runge-Kutta法のような高次の積分や、Symplectic積分のよう に性質の良い積分など、数値積分法には様々なものがある
11. 11 Riemann可積分である場合、積分値は分割の方法に依存しない 𝑥𝑘 𝑥𝑘+1 𝑓(𝑥𝑘 ) 𝑥𝑘 𝑥𝑘+1 𝑓(𝑥𝑘+1 ) 𝑥𝑘 𝑥𝑘+1 𝑓(𝑥𝑘+1 ) + 𝑓(𝑥𝑘 ) 2 →分割区間と、短冊の「高さ」の参照点をどこにとっても良い 参照点を左端にとる 参照点を右端にとる 平均値を使う
12. 12 ሶ 𝑥 = ෠ 𝑅 例えばこんな微分方程式を積分したい 𝑥 𝑡 − 𝑥 0 = න 0 𝑡 ෠ 𝑅𝑑𝑡 形式的に積分してみる න 0 𝑡 ෠ 𝑅𝑑𝑡 = lim 𝑁→∞ ෍ 𝑘=0 𝑁−1 ෠ 𝑅 𝑡𝑘 ℎ これを評価するために みたいなものを考えたいが、右辺はよくわからない
13. 13 න 0 𝑡 ෠ 𝑅𝑑𝑡 = 𝑥 𝑡 − 𝑥(0) こっちはよくわからないので こっちから考える ሶ 𝑥 = ෠ 𝑅 連続的なランダムウォークの代わりに 離散的なランダムウォークを考える 𝑥𝑘+1 = 𝑥𝑘 + ℎ ෠ 𝑅𝑘
14. 14 公平なコインを投げて、表なら右へ、裏なら左へ一歩進む 一歩の長さをh、一歩にかかる時間を𝜏としておく h x kステップ目の位置を𝑥𝑘 とすると、以下の式に従う 𝑥𝑘+1 = 𝑥𝑘 + ℎ ෠ 𝑅𝑘 ただし、 ෠ 𝑅𝑘 は確率1/2で1か-1を取る確率変数で、異なるステップ間は無相関 ෠ 𝑅𝑘 = 0 ෠ 𝑅𝑘 ෠ 𝑅𝑙 = 𝛿𝑘𝑙 ෠ 𝑅(𝑡) = 0 ෠ 𝑅(𝑡1 ) ෠ 𝑅(𝑡2 ) = 𝛿(𝑡1- 𝑡2) cf. 連続版の確率変数
15. 15 n回コインを投げた時、表が出る回数 ෠ 𝑋 𝑃 ෠ 𝑋 = 𝑘 = 𝑛 𝑘 1 2𝑛 これは確率1/2、試行回数n回の二項分布に従う ෠ 𝑋 = 𝑛/2 ෠ 𝑋2 = 𝑛/4 1次と2次のモーメントはすぐに計算できる
16. 16 n回コインを投げてX枚が表だったのなら、右にX歩、左にn-X歩進んでいる 𝑥𝑛 = 2 ෠ 𝑋 − 𝑛 ℎ 𝑥𝑛 = 0 𝑥𝑛 2 = 𝑛ℎ2 平均 分散 二項分布は試行回数が多い時、同じ平均、分散を持つ正規分布で近似できる 𝑃 𝑥𝑛 = 𝑥 = 1 2𝜋𝑛ℎ2 exp − 𝑥2 2𝑛ℎ2
17. 17 𝑥𝑛 2 = 𝑛ℎ2 から、 𝑛ℎ2 = 𝑡 𝑛𝜏 = 𝑡 (空間の連続極限) (時間の連続極限) を固定したままn無限大の極限をとる 𝑃 𝑥 < 𝑥 𝑡 < 𝑥 + 𝑑𝑥 = 1 2𝜋𝑡 exp − 𝑥2 2𝑡 𝑑𝑥 𝑃 𝑥𝑛 = 𝑥 = 1 2𝜋𝑛ℎ2 exp − 𝑥2 2𝑛ℎ2
18. 18 𝑥𝑘+1 = 𝑥𝑘 + ℎ ෠ 𝑅𝑘 𝑥𝑛 = ෍ 𝑘=0 𝑛 ෠ 𝑅𝑘 ℎ lim 𝑛→∞ 𝑥𝑛 = 𝑥(𝑡) ሶ 𝑥 = ෠ 𝑅 𝑥 𝑡 = න 0 𝑡 ෠ 𝑅𝑑𝑡 lim 𝑛→∞ ෍ 𝑘=0 𝑛 ෠ 𝑅𝑘 ℎ = න 0 𝑡 ෠ 𝑅𝑑𝑡 離散 連続 極限をとった時の対応 𝑥 𝑡 = න 0 𝑡 ෠ 𝑅𝑑𝑡 こっちが平均0、分散tの正規分布に 従うんだから こっちもそうでしょう
19. 19 𝑃 𝑥 < න 0 𝑡 ෠ 𝑅𝑑𝑡 < 𝑥 + 𝑑𝑥 = 1 2𝜋𝑡 exp − 𝑥2 2𝑡 𝑑𝑥 最終的に得られたもの න 0 𝑡 ෠ 𝑅𝑑𝑡 平均0、分散tの正規分布に従う確率変数 ※ 通常はWiener過程として「そういうものだ」と定義してしまう は とみなすことができる(※) 確率微分方程式の数値解法を構築できる
20. 20 例えばこんなLangevin方程式を解きたい ሶ 𝑣 = −𝛾𝑣 + 2𝐷 ෠ 𝑅 形式的に積分する 𝑣 𝑡 + ℎ − 𝑣 𝑡 = න 𝑡 𝑡+ℎ −𝛾𝑣 + 2𝐷 ෠ 𝑅 𝑑𝑡 時間刻みhが小さいと思って近似する 𝑣 𝑡 + ℎ − 𝑣 𝑡 ~ − 𝛾𝑣ℎ + 2𝐷 න 𝑡 𝑡+ℎ ෠ 𝑅𝑑𝑡 = −𝛾𝑣ℎ + 𝑤(0,2𝐷ℎ) 平均0、分散2Dhの 正規分布に従う乱数 これをEuler-Maruyama法と呼ぶ 𝑣 −𝛾𝑣
21. 21 ሶ 𝑣 = −𝛾𝑣 + 2𝐷 ෠ 𝑅 import numpy as np N = 100000 # ステップ数 h = 0.01 # 時間刻み v = 0 data = [] for _ in range(N): v = v - v * h + np.random.normal(0,np.sqrt(2.0*h)) data.append(v) コードで書くとこんな感じ (𝛾 = 𝐷 = 1) 𝑣(𝑡 + ℎ) = 𝑣(𝑡) − 𝛾𝑣ℎ + 𝑤(0,2𝐷ℎ) Euler-Maruyama法を適用
22. 22 ሶ 𝑣 = −𝛾𝑣 + 2𝐷 ෠ 𝑅 上記のLangevin方程式に従う系の定常分布はこうなる 𝑓 𝑣 = exp −𝛽𝑣2/2 2𝜋/𝛽 ↓確かになっている Langevin方程式と分布関数の時間発展の関係は？
23. 23 ሶ 𝑥 = 𝑣(𝑥) 場所によって走りやすさが異なるコースでのマラソンを考える →速度が場所に依存する この微分方程式は「どの場所でどのくらいの速度が出せるか」を表している →ランナー視点(Lagrange描像) x 𝑣(𝑥1 ) 𝑥1 𝑥2 𝑥3 𝑣(𝑥2 ) 𝑣(𝑥3 )
24. 24 ある時刻、ある場所でどのくらいランナーがいるかを考える 𝑓(𝑥, 𝑡) この分布関数は、時刻t、場所xのおけるランナーの密度を表している →観客視点(Euler描像) x ሶ 𝑥 = 𝑣(𝑥) 運動方程式と分布関数の時間発展の関係を知りたい 𝜕𝑓 𝜕𝑡 =? 各ランナーの位置と速度が わかっている時 ある場所で観測したランナーの数の 増減を知りたい
25. 25 𝑥 − ℎ 𝑥 + ℎ 𝑥 𝑓(𝑥, 𝑡) 𝑓(𝑥 − ℎ, 𝑡) 𝑣(𝑥 + ℎ) 𝑣(𝑥) 単位時間あたりに この区間の人数の 増減に着目 𝑓 𝑥 − ℎ 𝑣(𝑥) 出ていく人数 𝑓 𝑥 𝑣(𝑥 + ℎ) 人数の変化 Δ𝑓 𝑥, 𝑡 = 𝑓 𝑥 − ℎ 𝑣 𝑥 − 𝑓 𝑥 𝑣(𝑥 + ℎ) ~ 𝑓 − 𝜕𝑓 𝜕𝑥 ℎ 𝑣 − 𝑓 𝑣 + 𝜕𝑣 𝜕𝑥 ℎ = − 𝜕 𝜕𝑥 𝑣𝑓 ℎ 入ってくる人数 ሶ 𝑥 = 𝑣(𝑥) 𝜕𝑓 𝜕𝑡 = − 𝜕 𝜕𝑥 𝑣𝑓 𝑥 離散化して考える
26. 26 ሶ 𝑥 = ෠ 𝑅 拡散項の場合 𝑥 − ℎ 𝑥 + ℎ 𝑥 𝑥 𝑓(𝑥, 𝑡) 𝑓(𝑥 − ℎ, 𝑡) 𝑓(𝑥 + ℎ, 𝑡) 単位時間あたりに、αの割合だけ隣に移動する 𝛼 𝑓 𝑥 − ℎ + 𝑓(𝑥 + ℎ) 出ていく人数 入ってくる人数 2𝛼𝑓(𝑥) 人数の変化 Δ𝑓 𝑥, 𝑡 = 𝛼 𝑓 𝑥 − ℎ + 𝑓(𝑥 + ℎ) − 2𝛼𝑓(𝑥) ~𝛼 𝜕2𝑓 𝜕𝑥2 ℎ2 𝜕𝑓 𝜕𝑡 = 1 2 𝜕2𝑓 𝜕𝑥2 ሶ 𝑥 = ෠ 𝑅 ※なぜ1/2がつくかはここでは説明しません
27. 27 ሶ 𝑣 = −𝛾𝑣 + 2𝐷 ෠ 𝑅 Langevin方程式 対応する分布関数の時間発展 (Fokker-Planck方程式, FPE) 𝜕𝑓 𝜕𝑡 = − 𝜕 𝜕𝑣 −𝛾𝑣𝑓 + 𝐷 𝜕2𝑓 𝜕𝑣2 𝜕𝑓 𝜕𝑡 = 0 定常状態なら なので 𝜕 𝜕𝑣 𝛾𝑣𝑓 + 𝐷 𝜕2𝑓 𝜕𝑣2 = 0 𝜕𝑓 𝜕𝑣 = − 𝛾 𝐷 𝑣𝑓 整理すると これを解くと 𝑓 𝑣 = exp −𝛽𝑣2/2 2𝜋/𝛽 𝐻 = 1 2 𝑣2 に対するカノニカル分布になっている これは
28. 28 ሶ 𝑝 = − 𝜕𝐻 𝜕𝑞 − 𝛾 𝜕𝐻 𝜕𝑝 + 2𝐷 ෠ 𝑅 ሶ 𝑞 = 𝜕𝐻 𝜕𝑝 散逸項 揺動項 与えられたハミルトニアンHに対して以下の運動方程式を解くことで 定常状態として指定温度のカノニカル分布 を実現する数値計算手法 𝑓(𝑝, 𝑞)~exp −𝛽𝐻
29. 29 ሶ 𝑣 = −𝛾𝑣 + 2𝐷 ෠ 𝑅 以下のLangevin方程式では、拡散係数が定数だった 確率変数の係数が定数である場合を相加ノイズ(additive noise)と呼ぶ ሶ 𝑣 = −𝛾𝑣 + 𝐷(𝑣) ෠ 𝑅 一般には、拡散係数は速度に依存する 係数に状態依存性がある場合を相乗ノイズ(multiplicative noise)と呼ぶ 𝑣(𝑡)に時間依存性があるので、当然𝐷(𝑣) にも時間依存性が生じる
30. 30 ሶ 𝑥 = 𝑓(𝑥) ෠ 𝑅(𝑡) 相乗的なノイズを持つ確率微分方程式を考える න 𝑡 𝑡+ℎ 𝑓 𝑥 ෠ 𝑅 𝑡 𝑑𝑡 ~𝑓 𝑥𝑡 𝑤 小さい時間刻みについて、こんな近似をしたくなる 𝑡 𝑡 + ℎ 𝑓 𝑥𝑡 𝑤 これは、短冊の参照点として左端を採用したことに対応する この定義を伊藤積分、 この解釈を伊藤解釈と呼ぶ
31. 31 ሶ 𝑥 = 𝑓(𝑥) ෠ 𝑅(𝑡) න 𝑡 𝑡+ℎ 𝑓 𝑥 ෠ 𝑅 𝑡 𝑑𝑡 ~ 𝑓 𝑥𝑡+ℎ + 𝑓 𝑥𝑡 2 𝑤 左端ではなく、両側の値の平均値を使ってみる 𝑡 𝑡 + ℎ 𝑓 𝑥𝑡 𝑤 𝑓 𝑥𝑡+ℎ 𝑤 𝑓 𝑥𝑡+ℎ + 𝑓 𝑥𝑡 2 この定義をStratonovich積分、 この解釈をStratonovich解釈と呼ぶ
32. 32 න 𝑡 𝑡+ℎ 𝑓 𝑥 ෠ 𝑅 𝑡 𝑑𝑡 ~ 𝑓 𝑥𝑡+ℎ + 𝑓 𝑥𝑡 2 𝑤 න 𝑡 𝑡+ℎ 𝑓 𝑥 ෠ 𝑅 𝑡 𝑑𝑡 ~𝑓 𝑥𝑡 𝑤 伊藤積分 Stratonovich積分 Riemann可積分の場合は両者の極限は一致するが、 乗法的ノイズを持つ系では両者は一致しない どの解釈を採用するかで結果が変わる
33. 33 ሶ 𝑣 = −𝛾𝑣 + 2𝐷 ෠ 𝑅 Langevin方程式 𝜕𝑓 𝜕𝑡 = − 𝜕 𝜕𝑣 −𝛾𝑣𝑓 + 𝜕2𝐷𝑓 𝜕𝑣2 伊藤解釈のFPE Stratonovich解釈のFPE 𝜕𝑓 𝜕𝑡 = − 𝜕 𝜕𝑣 −𝛾𝑣𝑓 + 𝜕 𝜕𝑣 𝐷 𝜕 𝜕𝑣 𝐷𝑓 Dが定数(相加ノイズ)なら両者は一致する 𝜕𝑓 𝜕𝑡 = − 𝜕 𝜕𝑣 −𝛾𝑣𝑓 + 𝐷 𝜕2𝑓 𝜕𝑣2 両者は時間発展も定常状態も異なる
34. 34 ሶ 𝑥 = 𝑓(𝑥) ෠ 𝑅(𝑡) 𝑥𝑡+ℎ = 𝑥𝑡 + 𝑓 𝑥𝑡+ℎ + 𝑓 𝑥𝑡 2 𝑤 Stratonovich解釈では 未知の値が右辺に表れてしまう 伊藤解釈では 𝑥𝑡+ℎ = 𝑥𝑡 + 𝑓 𝑥𝑡 𝑤 伊藤解釈で未来の値を予想し、それを使って改めてStratonovich解釈で計算する (Two-Step法) 𝑥𝑡+ℎ = 𝑥𝑡 + 𝑓 𝑥𝑡+ℎ 𝐼 + 𝑓 𝑥𝑡 2 𝑤 𝑥𝑡+ℎ 𝐼 = 𝑥𝑡 + 𝑓 𝑥𝑡 𝑤 同じ値を使う必要が あることに注意
35. 35 ሶ 𝑥 = −𝑥3 + 2𝑥 ෠ 𝑅1 + 2 ෠ 𝑅2 以下の相乗的ノイズを持つLangevin方程式を考える R1とR2は独立な白色雑音 R2はエルゴード性を満たすために追加 伊藤解釈 𝜕𝑓 𝜕𝑡 = − 𝜕 𝜕𝑡 −𝑥3𝑓 − 𝜕 𝜕𝑡 𝑥2𝑓 − 𝜕𝑓 𝜕𝑡 𝑓I ~ exp −𝑥2 1 + 𝑥2 Fokker-Planck方程式 定常状態 ෠ 𝑅1 由来 ෠ 𝑅2 由来 Stratonovich解釈 𝜕𝑓 𝜕𝑡 = − 𝜕 𝜕𝑡 −𝑥3𝑓 − 𝑥 𝜕 𝜕𝑡 𝑥𝑓 − 𝜕𝑓 𝜕𝑡 ෠ 𝑅1 由来 ෠ 𝑅2 由来 𝑓S ~ exp −𝑥2 相乗ノイズ 相加ノイズ
36. 36 import numpy as np N = 100000 # ステップ数 h = 0.01 # 時間刻み x = 0 data = [] for _ in range(N): w1 = np.random.normal(0,np.sqrt(2*h)) # R1に対応する白色雑音 w2 = np.random.normal(0,np.sqrt(2*h)) # R2に対応する白色雑音 x_i = x + (-x * x * x) * h + x * w1 + w2 x = x + (-x * x * x) * h + (x+x_i)*0.5*w1 + w2 data.append(x) ሶ 𝑥 = −𝑥3 + 2𝑥 ෠ 𝑅1 + 2 ෠ 𝑅2 Two-step法の実装 似た計算を二度行うが、両者で「同じノイズ」を使う必要がある
37. 37 𝑓I ~ exp −𝑥2 1 + 𝑥2 𝑓S ~ exp −𝑥2 伊藤解釈 Stratonovich解釈 同じ確率微分方程式に対して、異なる解釈を採用すると、異なる定常分布を得る
38. 38 • Langevin法は、ハミルトンの運動方程式に散逸項と揺動項を 加えて指定温度のカノニカル分布を実現する数値計算手法 • 白色雑音の積分はガウス分布に従う確率変数と解釈する • 運動方程式に対応する分布関数の時間発展方程式を得る (Fokker-Planck方程式) • 運動方程式はLagrange描像 • Fokker-Planck方程式はEuler描像 • 確率微分方程式の解釈には任意性がある • 広く使われているのは伊藤解釈とStratonovich解釈 • 相加的ノイズなら両者は一致 • 相乗的ノイズの場合は異なる結果を与える