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

Slide 1

Slide 1 text

1 Langevin法と確率微分方程式慶應義塾大学理工学部物理情報工学科渡辺 2021/07/02

Slide 2

Slide 2 text

2 確率微分方程式の厳密な取り扱いはかなり面倒ですここで紹介する確率微分方程式の解釈は、「とりあえずそう思っておく」程度と認識してください

Slide 3

Slide 3 text

3 分子動力学法(Molecular Dynamics Method, MD)とは狭義にはハミルトンの運動方程式を数値積分する手法 ሶ 𝑝 = − 𝜕𝐻 𝜕𝑞 ሶ 𝑞 = 𝜕𝐻 𝜕𝑝 𝑑𝐻 𝑑𝑡 = 𝜕𝐻 𝜕𝑝 ሶ 𝑝 + 𝜕𝐻 𝜕𝑞 ሶ 𝑞 = 0 より、この時間発展ではエネルギーが保存される (NVEアンサンブル) 我々が興味があるのは、物理量のエネルギー依存性ではなく温度依存性温度を制御したい

Slide 4

Slide 4 text

4 ハミルトンの運動方程式に、温度制御項を付け加える ሶ 𝑝 = − 𝜕𝐻 𝜕𝑞 − 𝛾 𝜕𝐻 𝜕𝑝 + 2𝐷 ෠ 𝑅 ሶ 𝑞 = 𝜕𝐻 𝜕𝑝 散逸項揺動項散逸項：エネルギーを減らそうとする揺動項：エネルギーを注入する散逸と揺動が釣り合ったところで温度が一定となる → Langevin法温度制御項

Slide 5

Slide 5 text

5 ሶ 𝑝 = − 𝜕𝐻 𝜕𝑞 − 𝛾 𝜕𝐻 𝜕𝑝 + 2𝐷 ෠ 𝑅 ሶ 𝑞 = 𝜕𝐻 𝜕𝑝 揺動項揺動項のRは時間に依存する確率的な変数であり、以下を満たす ෠ 𝑅(𝑡) = 0 ෠ 𝑅(𝑡1 ) ෠ 𝑅(𝑡2 ) = 𝛿(𝑡1- 𝑡2) 平均ゼロ無相関(白色雑音) このような確率変数を含む微分方程式を確率微分方程式と呼ぶ特に、ブラウン運動を記述する確率微分方程式をLangevin方程式と呼ぶ

Slide 6

Slide 6 text

6 ሶ 𝑝 = − 𝜕𝐻 𝜕𝑞 − 𝛾 𝜕𝐻 𝜕𝑝 + 2𝐷 ෠ 𝑅 ሶ 𝑞 = 𝜕𝐻 𝜕𝑝 上記の時間発展に従う系の定常分布はカノニカル分布となる 𝑓(𝑝, 𝑞)~exp −𝛽𝐻 摩擦係数と拡散係数の比が温度を与える(揺動散逸定理) 𝛽 = 1 𝑘𝑇 = 𝛾/𝐷 (NVTアンサンブル)

Slide 7

Slide 7 text

7 ሶ 𝑝 = − 𝜕𝐻 𝜕𝑞 − 𝛾 𝜕𝐻 𝜕𝑝 + 2𝐷 ෠ 𝑅 ሶ 𝑞 = 𝜕𝐻 𝜕𝑝 Langevin法は、温度制御項のついた運動方程式を数値積分する手法揺動項 Rは確率変数 • 確率変数を積分するにはどうすればよいか？ • 運動方程式と分布関数の関係はどうなっているか？

Slide 8

Slide 8 text

8 積分区間を適当に分割し、以下の和をとる(Riemann和) න 𝑥𝑠 𝑥𝑒 𝑓 𝑥 𝑑𝑥 関数f(x)の定積分を求めたい O x f(x) O x f(x) ≒

Slide 9

Slide 9 text

9 分割の最大幅をδとする任意の正数εに対してあるδが存在し、 ෍ 𝑖 𝑛−1 𝑓 𝑥𝑖 𝑥𝑖+1 − 𝑥𝑖 − 𝐴 < 𝜀 を満たすAを定義できる場合、f(x)はRiemann可積分であり、 Aをf(x)の定積分、もしくはRiemann積分と呼ぶ 𝛿 = max(𝑥𝑖+1 − 𝑥𝑖 ) න 𝑥𝑠 𝑥𝑒 𝑓 𝑥 𝑑𝑡 ≡ 𝐴 簡単にいえば、短冊の和の極限がRiemann積分我々が通常使う「積分」とはRiemann積分のこと

Slide 10

Slide 10 text

10 普通に積分できる場合、普通に数値積分できる ሶ 𝑥 = 𝑓(𝑥) を数値積分したい 𝑥 𝑡 + ℎ = 𝑥 𝑡 + ሶ 𝑥 𝑡 ℎ + 𝑂(ℎ2) 数値積分とは𝑥 𝑡 と ሶ 𝑥 𝑡 から 𝑥 𝑡 + ℎ を求めること න 𝑡 𝑡+ℎ ሶ 𝑥𝑑𝑡 = 𝑥 𝑡 + ℎ − 𝑥(𝑡) ←これは厳密に成り立つが一般に積分は厳密に求められない Taylor展開して一次までとる無視する 𝑥 𝑡 と ሶ 𝑥 𝑡 から 𝑥 𝑡 + ℎ が近似的に求まった(一次のEuler法) Runge-Kutta法のような高次の積分や、Symplectic積分のように性質の良い積分など、数値積分法には様々なものがある

Slide 11

Slide 11 text

11 Riemann可積分である場合、積分値は分割の方法に依存しない 𝑥𝑘 𝑥𝑘+1 𝑓(𝑥𝑘 ) 𝑥𝑘 𝑥𝑘+1 𝑓(𝑥𝑘+1 ) 𝑥𝑘 𝑥𝑘+1 𝑓(𝑥𝑘+1 ) + 𝑓(𝑥𝑘 ) 2 →分割区間と、短冊の「高さ」の参照点をどこにとっても良い参照点を左端にとる参照点を右端にとる平均値を使う

Slide 12

Slide 12 text

12 ሶ 𝑥 = ෠ 𝑅 例えばこんな微分方程式を積分したい 𝑥 𝑡 − 𝑥 0 = න 0 𝑡 ෠ 𝑅𝑑𝑡 形式的に積分してみる න 0 𝑡 ෠ 𝑅𝑑𝑡 = lim 𝑁→∞ ෍ 𝑘=0 𝑁−1 ෠ 𝑅 𝑡𝑘 ℎ これを評価するためにみたいなものを考えたいが、右辺はよくわからない

Slide 13

Slide 13 text

13 න 0 𝑡 ෠ 𝑅𝑑𝑡 = 𝑥 𝑡 − 𝑥(0) こっちはよくわからないのでこっちから考える ሶ 𝑥 = ෠ 𝑅 連続的なランダムウォークの代わりに離散的なランダムウォークを考える 𝑥𝑘+1 = 𝑥𝑘 + ℎ ෠ 𝑅𝑘

Slide 14

Slide 14 text

14 公平なコインを投げて、表なら右へ、裏なら左へ一歩進む一歩の長さをh、一歩にかかる時間を𝜏としておく h x kステップ目の位置を𝑥𝑘 とすると、以下の式に従う 𝑥𝑘+1 = 𝑥𝑘 + ℎ ෠ 𝑅𝑘 ただし、 ෠ 𝑅𝑘 は確率1/2で1か-1を取る確率変数で、異なるステップ間は無相関 ෠ 𝑅𝑘 = 0 ෠ 𝑅𝑘 ෠ 𝑅𝑙 = 𝛿𝑘𝑙 ෠ 𝑅(𝑡) = 0 ෠ 𝑅(𝑡1 ) ෠ 𝑅(𝑡2 ) = 𝛿(𝑡1- 𝑡2) cf. 連続版の確率変数

Slide 15

Slide 15 text

15 n回コインを投げた時、表が出る回数 ෠ 𝑋 𝑃 ෠ 𝑋 = 𝑘 = 𝑛 𝑘 1 2𝑛 これは確率1/2、試行回数n回の二項分布に従う ෠ 𝑋 = 𝑛/2 ෠ 𝑋2 = 𝑛/4 1次と2次のモーメントはすぐに計算できる

Slide 16

Slide 16 text

16 n回コインを投げてX枚が表だったのなら、右にX歩、左にn-X歩進んでいる 𝑥𝑛 = 2 ෠ 𝑋 − 𝑛 ℎ 𝑥𝑛 = 0 𝑥𝑛 2 = 𝑛ℎ2 平均分散二項分布は試行回数が多い時、同じ平均、分散を持つ正規分布で近似できる 𝑃 𝑥𝑛 = 𝑥 = 1 2𝜋𝑛ℎ2 exp − 𝑥2 2𝑛ℎ2

Slide 17

Slide 17 text

17 𝑥𝑛 2 = 𝑛ℎ2 から、 𝑛ℎ2 = 𝑡 𝑛𝜏 = 𝑡 (空間の連続極限) (時間の連続極限) を固定したままn無限大の極限をとる 𝑃 𝑥 < 𝑥 𝑡 < 𝑥 + 𝑑𝑥 = 1 2𝜋𝑡 exp − 𝑥2 2𝑡 𝑑𝑥 𝑃 𝑥𝑛 = 𝑥 = 1 2𝜋𝑛ℎ2 exp − 𝑥2 2𝑛ℎ2

Slide 18

Slide 18 text

18 𝑥𝑘+1 = 𝑥𝑘 + ℎ ෠ 𝑅𝑘 𝑥𝑛 = ෍ 𝑘=0 𝑛 ෠ 𝑅𝑘 ℎ lim 𝑛→∞ 𝑥𝑛 = 𝑥(𝑡) ሶ 𝑥 = ෠ 𝑅 𝑥 𝑡 = න 0 𝑡 ෠ 𝑅𝑑𝑡 lim 𝑛→∞ ෍ 𝑘=0 𝑛 ෠ 𝑅𝑘 ℎ = න 0 𝑡 ෠ 𝑅𝑑𝑡 離散連続極限をとった時の対応 𝑥 𝑡 = න 0 𝑡 ෠ 𝑅𝑑𝑡 こっちが平均0、分散tの正規分布に従うんだからこっちもそうでしょう

Slide 19

Slide 19 text

19 𝑃 𝑥 < න 0 𝑡 ෠ 𝑅𝑑𝑡 < 𝑥 + 𝑑𝑥 = 1 2𝜋𝑡 exp − 𝑥2 2𝑡 𝑑𝑥 最終的に得られたもの න 0 𝑡 ෠ 𝑅𝑑𝑡 平均0、分散tの正規分布に従う確率変数 ※ 通常はWiener過程として「そういうものだ」と定義してしまうはとみなすことができる(※) 確率微分方程式の数値解法を構築できる

Slide 20

Slide 20 text

20 例えばこんなLangevin方程式を解きたい ሶ 𝑣 = −𝛾𝑣 + 2𝐷 ෠ 𝑅 形式的に積分する 𝑣 𝑡 + ℎ − 𝑣 𝑡 = න 𝑡 𝑡+ℎ −𝛾𝑣 + 2𝐷 ෠ 𝑅 𝑑𝑡 時間刻みhが小さいと思って近似する 𝑣 𝑡 + ℎ − 𝑣 𝑡 ~ − 𝛾𝑣ℎ + 2𝐷 න 𝑡 𝑡+ℎ ෠ 𝑅𝑑𝑡 = −𝛾𝑣ℎ + 𝑤(0,2𝐷ℎ) 平均0、分散2Dhの正規分布に従う乱数これをEuler-Maruyama法と呼ぶ 𝑣 −𝛾𝑣

Slide 21

Slide 21 text

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法を適用

Slide 22

Slide 22 text

22 ሶ 𝑣 = −𝛾𝑣 + 2𝐷 ෠ 𝑅 上記のLangevin方程式に従う系の定常分布はこうなる 𝑓 𝑣 = exp −𝛽𝑣2/2 2𝜋/𝛽 ↓確かになっている Langevin方程式と分布関数の時間発展の関係は？

Slide 23

Slide 23 text

23 ሶ 𝑥 = 𝑣(𝑥) 場所によって走りやすさが異なるコースでのマラソンを考える →速度が場所に依存するこの微分方程式は「どの場所でどのくらいの速度が出せるか」を表している →ランナー視点(Lagrange描像) x 𝑣(𝑥1 ) 𝑥1 𝑥2 𝑥3 𝑣(𝑥2 ) 𝑣(𝑥3 )

Slide 24

Slide 24 text

24 ある時刻、ある場所でどのくらいランナーがいるかを考える 𝑓(𝑥, 𝑡) この分布関数は、時刻t、場所xのおけるランナーの密度を表している →観客視点(Euler描像) x ሶ 𝑥 = 𝑣(𝑥) 運動方程式と分布関数の時間発展の関係を知りたい 𝜕𝑓 𝜕𝑡 =? 各ランナーの位置と速度がわかっている時ある場所で観測したランナーの数の増減を知りたい

Slide 25

Slide 25 text

25 𝑥 − ℎ 𝑥 + ℎ 𝑥 𝑓(𝑥, 𝑡) 𝑓(𝑥 − ℎ, 𝑡) 𝑣(𝑥 + ℎ) 𝑣(𝑥) 単位時間あたりにこの区間の人数の増減に着目 𝑓 𝑥 − ℎ 𝑣(𝑥) 出ていく人数 𝑓 𝑥 𝑣(𝑥 + ℎ) 人数の変化 Δ𝑓 𝑥, 𝑡 = 𝑓 𝑥 − ℎ 𝑣 𝑥 − 𝑓 𝑥 𝑣(𝑥 + ℎ) ~ 𝑓 − 𝜕𝑓 𝜕𝑥 ℎ 𝑣 − 𝑓 𝑣 + 𝜕𝑣 𝜕𝑥 ℎ = − 𝜕 𝜕𝑥 𝑣𝑓 ℎ 入ってくる人数 ሶ 𝑥 = 𝑣(𝑥) 𝜕𝑓 𝜕𝑡 = − 𝜕 𝜕𝑥 𝑣𝑓 𝑥 離散化して考える

Slide 26

Slide 26 text

26 ሶ 𝑥 = ෠ 𝑅 拡散項の場合 𝑥 − ℎ 𝑥 + ℎ 𝑥 𝑥 𝑓(𝑥, 𝑡) 𝑓(𝑥 − ℎ, 𝑡) 𝑓(𝑥 + ℎ, 𝑡) 単位時間あたりに、αの割合だけ隣に移動する 𝛼 𝑓 𝑥 − ℎ + 𝑓(𝑥 + ℎ) 出ていく人数入ってくる人数 2𝛼𝑓(𝑥) 人数の変化 Δ𝑓 𝑥, 𝑡 = 𝛼 𝑓 𝑥 − ℎ + 𝑓(𝑥 + ℎ) − 2𝛼𝑓(𝑥) ~𝛼 𝜕2𝑓 𝜕𝑥2 ℎ2 𝜕𝑓 𝜕𝑡 = 1 2 𝜕2𝑓 𝜕𝑥2 ሶ 𝑥 = ෠ 𝑅 ※なぜ1/2がつくかはここでは説明しません

Slide 27

Slide 27 text

27 ሶ 𝑣 = −𝛾𝑣 + 2𝐷 ෠ 𝑅 Langevin方程式対応する分布関数の時間発展 (Fokker-Planck方程式, FPE) 𝜕𝑓 𝜕𝑡 = − 𝜕 𝜕𝑣 −𝛾𝑣𝑓 + 𝐷 𝜕2𝑓 𝜕𝑣2 𝜕𝑓 𝜕𝑡 = 0 定常状態ならなので 𝜕 𝜕𝑣 𝛾𝑣𝑓 + 𝐷 𝜕2𝑓 𝜕𝑣2 = 0 𝜕𝑓 𝜕𝑣 = − 𝛾 𝐷 𝑣𝑓 整理するとこれを解くと 𝑓 𝑣 = exp −𝛽𝑣2/2 2𝜋/𝛽 𝐻 = 1 2 𝑣2 に対するカノニカル分布になっているこれは

Slide 28

Slide 28 text

28 ሶ 𝑝 = − 𝜕𝐻 𝜕𝑞 − 𝛾 𝜕𝐻 𝜕𝑝 + 2𝐷 ෠ 𝑅 ሶ 𝑞 = 𝜕𝐻 𝜕𝑝 散逸項揺動項与えられたハミルトニアンHに対して以下の運動方程式を解くことで定常状態として指定温度のカノニカル分布を実現する数値計算手法 𝑓(𝑝, 𝑞)~exp −𝛽𝐻

Slide 29

Slide 29 text

29 ሶ 𝑣 = −𝛾𝑣 + 2𝐷 ෠ 𝑅 以下のLangevin方程式では、拡散係数が定数だった確率変数の係数が定数である場合を相加ノイズ(additive noise)と呼ぶ ሶ 𝑣 = −𝛾𝑣 + 𝐷(𝑣) ෠ 𝑅 一般には、拡散係数は速度に依存する係数に状態依存性がある場合を相乗ノイズ(multiplicative noise)と呼ぶ 𝑣(𝑡)に時間依存性があるので、当然𝐷(𝑣) にも時間依存性が生じる

Slide 30

Slide 30 text

30 ሶ 𝑥 = 𝑓(𝑥) ෠ 𝑅(𝑡) 相乗的なノイズを持つ確率微分方程式を考える න 𝑡 𝑡+ℎ 𝑓 𝑥 ෠ 𝑅 𝑡 𝑑𝑡 ~𝑓 𝑥𝑡 𝑤 小さい時間刻みについて、こんな近似をしたくなる 𝑡 𝑡 + ℎ 𝑓 𝑥𝑡 𝑤 これは、短冊の参照点として左端を採用したことに対応するこの定義を伊藤積分、この解釈を伊藤解釈と呼ぶ

Slide 31

Slide 31 text

31 ሶ 𝑥 = 𝑓(𝑥) ෠ 𝑅(𝑡) න 𝑡 𝑡+ℎ 𝑓 𝑥 ෠ 𝑅 𝑡 𝑑𝑡 ~ 𝑓 𝑥𝑡+ℎ + 𝑓 𝑥𝑡 2 𝑤 左端ではなく、両側の値の平均値を使ってみる 𝑡 𝑡 + ℎ 𝑓 𝑥𝑡 𝑤 𝑓 𝑥𝑡+ℎ 𝑤 𝑓 𝑥𝑡+ℎ + 𝑓 𝑥𝑡 2 この定義をStratonovich積分、この解釈をStratonovich解釈と呼ぶ

Slide 32

Slide 32 text

32 න 𝑡 𝑡+ℎ 𝑓 𝑥 ෠ 𝑅 𝑡 𝑑𝑡 ~ 𝑓 𝑥𝑡+ℎ + 𝑓 𝑥𝑡 2 𝑤 න 𝑡 𝑡+ℎ 𝑓 𝑥 ෠ 𝑅 𝑡 𝑑𝑡 ~𝑓 𝑥𝑡 𝑤 伊藤積分 Stratonovich積分 Riemann可積分の場合は両者の極限は一致するが、乗法的ノイズを持つ系では両者は一致しないどの解釈を採用するかで結果が変わる

Slide 33

Slide 33 text

33 ሶ 𝑣 = −𝛾𝑣 + 2𝐷 ෠ 𝑅 Langevin方程式 𝜕𝑓 𝜕𝑡 = − 𝜕 𝜕𝑣 −𝛾𝑣𝑓 + 𝜕2𝐷𝑓 𝜕𝑣2 伊藤解釈のFPE Stratonovich解釈のFPE 𝜕𝑓 𝜕𝑡 = − 𝜕 𝜕𝑣 −𝛾𝑣𝑓 + 𝜕 𝜕𝑣 𝐷 𝜕 𝜕𝑣 𝐷𝑓 Dが定数(相加ノイズ)なら両者は一致する 𝜕𝑓 𝜕𝑡 = − 𝜕 𝜕𝑣 −𝛾𝑣𝑓 + 𝐷 𝜕2𝑓 𝜕𝑣2 両者は時間発展も定常状態も異なる

Slide 34

Slide 34 text

34 ሶ 𝑥 = 𝑓(𝑥) ෠ 𝑅(𝑡) 𝑥𝑡+ℎ = 𝑥𝑡 + 𝑓 𝑥𝑡+ℎ + 𝑓 𝑥𝑡 2 𝑤 Stratonovich解釈では未知の値が右辺に表れてしまう伊藤解釈では 𝑥𝑡+ℎ = 𝑥𝑡 + 𝑓 𝑥𝑡 𝑤 伊藤解釈で未来の値を予想し、それを使って改めてStratonovich解釈で計算する (Two-Step法) 𝑥𝑡+ℎ = 𝑥𝑡 + 𝑓 𝑥𝑡+ℎ 𝐼 + 𝑓 𝑥𝑡 2 𝑤 𝑥𝑡+ℎ 𝐼 = 𝑥𝑡 + 𝑓 𝑥𝑡 𝑤 同じ値を使う必要があることに注意

Slide 35

Slide 35 text

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 相乗ノイズ相加ノイズ

Slide 36

Slide 36 text

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法の実装似た計算を二度行うが、両者で「同じノイズ」を使う必要がある

Slide 37

Slide 37 text

37 𝑓I ~ exp −𝑥2 1 + 𝑥2 𝑓S ~ exp −𝑥2 伊藤解釈 Stratonovich解釈同じ確率微分方程式に対して、異なる解釈を採用すると、異なる定常分布を得る

Slide 38

Slide 38 text

38 • Langevin法は、ハミルトンの運動方程式に散逸項と揺動項を加えて指定温度のカノニカル分布を実現する数値計算手法 • 白色雑音の積分はガウス分布に従う確率変数と解釈する • 運動方程式に対応する分布関数の時間発展方程式を得る (Fokker-Planck方程式) • 運動方程式はLagrange描像 • Fokker-Planck方程式はEuler描像 • 確率微分方程式の解釈には任意性がある • 広く使われているのは伊藤解釈とStratonovich解釈 • 相加的ノイズなら両者は一致 • 相乗的ノイズの場合は異なる結果を与える