分子動力学法(1) 理論的背景と実装の基礎 / Simulation 05

シミュレーション工学

Transcript

1. 1 85 分子動力学法(1)理論的背景と数値積分 シミュレーション工学 慶應義塾大学大学院理工学研究科基礎理工学専攻物理情報専修 渡辺宙志
2. 2 85 • 解析力学の復習 • 分子動力学法の時間発展アルゴリズム 本講義の目的 分子動力学法とは • 原子や分子にかかる「力」を計算し、位置 を更新していく数値計算手法 • 背景にある理論は解析力学
3. 3 85 𝑚𝑎 = 𝐹 物体の加速度は力に比例する 物体の動きにくさは質量に比例する
4. 4 85 𝑎 = 𝑑𝑣 𝑑𝑡 𝑣 = 𝑑𝑟 𝑑𝑡 加速度とは速度の時間変化率 速度とは位置の時間変化率 𝑑2𝑟 𝑑𝑡2 = 𝐹 𝑚 運動方程式は 位置の二階微分方程式
5. 5 85 O 𝑟 自然長の位置を原点にとる 質量とバネ定数は1とする 𝑑2𝑟 𝑑𝑡2 = −𝑟
6. 6 85 𝑟 = 0 物体が原点にいる これだけでは系の状態が定まらない 止まった状態かもしれない 運動中にちょうど原点を 通過したところかもしれない
7. 7 85 𝑟 = 0,𝑣 = 1 物体が原点にいて、速度が1 𝑑2𝑟 𝑑𝑡2 = −𝑟 二階の常微分方程式で初期条件を二つ与えたので、 この系の運動は完全に定まる
8. 8 85 (𝑣, 𝑟) = (1,0) 𝑣 𝑟 O この系の状態は(𝑣, 𝑟)空間の点を指定することで一意に定まる この空間を位相空間(Phase Space)と呼ぶ
9. 9 85 系の状態は(𝑣, 𝑟)空間で指定できるのだから 運動方程式も(𝑣, 𝑟)で書きたい 𝑑2𝑟 𝑑𝑡2 = −𝑟 𝑑𝑣 𝑑𝑡 = −𝑟 二階の常微分方程式を一階の連立常微分方程式に 𝑑𝑟 𝑑𝑡 = 𝑣
10. 10 85 変数の組(𝑣, 𝑟)の時間微分( ሶ 𝑣, ሶ 𝑟)が (𝑣, 𝑟)の関数として書けている 𝑑𝑣 𝑑𝑡 = 𝑓𝑣 (𝑟, 𝑣) 𝑑𝑟 𝑑𝑡 = 𝑓𝑟 (𝑟, 𝑣) 𝑓𝑟 𝑟, 𝑣 = 𝑣 𝑓𝑣 𝑟, 𝑣 = −𝑟 このような系を力学系(dynamical system)と呼ぶ 一般的に書くと ※運動方程式は力学系の一種 𝑑𝑣 𝑑𝑡 = −𝑟 𝑑𝑟 𝑑𝑡 = 𝑣
11. 11 85 (𝑣, 𝑟) = (1,0) 𝑟 𝑣 O 現在の状態(原点で右向きに運動) 時間微分(右向きに運動し、加速度0) ( ሶ 𝑣, ሶ 𝑟) = (0,1) ( ሶ 𝑣, ሶ 𝑟) = (−𝑟, 𝑣)
12. 12 85 現在の状態(バネが伸び切って停止) (𝑣, 𝑟) = (0,1) 時間微分(左向きに加速し、速度0) ( ሶ 𝑣, ሶ 𝑟) = (−1,0) 𝑟 𝑣 O ( ሶ 𝑣, ሶ 𝑟) = (−𝑟, 𝑣)
13. 13 85 𝑟 𝑣 O 今の手続きを様々な点で繰り返すと・・・ 運動方程式は位相空間にベクトル場を定義する 回転する速度場
14. 14 85 𝑟 𝑣 O (𝑣, 𝑟) = (1,0) 初期条件：位相空間にトレーサーを置くこと 方程式の解：トレーサーの軌跡
15. 15 85 • 運動方程式に従う系の状態を一意に決め ることができる空間を位相空間と呼ぶ • 運動方程式とは、位相空間にベクトル場 (速度場)を定義するものである • 系の状態は、位相空間上の一点で表現で きる • 系の運動は、位相空間上の軌跡で表現さ れ、運動方程式が定義した速度場に従っ て動く点の軌跡である
16. 16 85 3次元空間に𝑁粒子がいる場合、その位相空間は6𝑁次元 (位置が3成分、速度が3成分で粒子あたり6成分) (𝑟𝑥 1, 𝑟𝑦 1, 𝑟𝑧 1, ⋯ , 𝑟𝑥 𝑁, 𝑟𝑦 𝑁, 𝑟𝑧 𝑁, 𝑣𝑥 1, 𝑣𝑦 1, 𝑣𝑧 1, ⋯ , 𝑣𝑥 𝑁, 𝑣𝑦 𝑁, 𝑣𝑧 𝑁) 6𝑁 6𝑁次元の空間にベクトル場を定義するには、6𝑁本の関数が必要 𝑑𝑟𝑥 1 𝑑𝑡 = 𝑓1 (𝑟𝑥 1, 𝑟𝑦 1, 𝑟𝑧 1, ⋯ , 𝑟𝑥 𝑁, 𝑟𝑦 𝑁, 𝑟𝑧 𝑁, 𝑣𝑥 1, 𝑣𝑦 1, 𝑣𝑧 1, ⋯ , 𝑣𝑥 𝑁, 𝑣𝑦 𝑁, 𝑣𝑧 𝑁) 𝑑𝑣𝑧 𝑁 𝑑𝑡 = 𝑓6𝑁 (𝑟𝑥 1, 𝑟𝑦 1, 𝑟𝑧 1, ⋯ , 𝑟𝑥 𝑁, 𝑟𝑦 𝑁, 𝑟𝑧 𝑁, 𝑣𝑥 1, 𝑣𝑦 1, 𝑣𝑧 1, ⋯ , 𝑣𝑥 𝑁, 𝑣𝑦 𝑁, 𝑣𝑧 𝑁) ・・・ 6𝑁
17. 17 85 ℒ(𝑟𝑥 1, 𝑟𝑦 1, 𝑟𝑧 1, ⋯ , 𝑟𝑥 𝑁, 𝑟𝑦 𝑁, 𝑟𝑧 𝑁, 𝑣𝑥 1, 𝑣𝑦 1, 𝑣𝑧 1, ⋯ , 𝑣𝑥 𝑁, 𝑣𝑦 𝑁, 𝑣𝑧 𝑁) ラグランジアンというスカラー関数ひとつから、運動方 程式に必要な6𝑁本の方程式を生み出すことができる 𝑑𝑟𝑥 1 𝑑𝑡 = 𝑓1 (𝑟𝑥 1, 𝑟𝑦 1, 𝑟𝑧 1, ⋯ , 𝑟𝑥 𝑁, 𝑟𝑦 𝑁, 𝑟𝑧 𝑁, 𝑣𝑥 1, 𝑣𝑦 1, 𝑣𝑧 1, ⋯ , 𝑣𝑥 𝑁, 𝑣𝑦 𝑁, 𝑣𝑧 𝑁) 𝑑𝑣𝑧 𝑁 𝑑𝑡 = 𝑓6𝑁 (𝑟𝑥 1, 𝑟𝑦 1, 𝑟𝑧 1, ⋯ , 𝑟𝑥 𝑁, 𝑟𝑦 𝑁, 𝑟𝑧 𝑁, 𝑣𝑥 1, 𝑣𝑦 1, 𝑣𝑧 1, ⋯ , 𝑣𝑥 𝑁, 𝑣𝑦 𝑁, 𝑣𝑧 𝑁) ・・・ 6𝑁
18. 18 85 𝐼 = න 𝑡0 𝑡1 ℒ( ሶ 𝑟, 𝑟)𝑑𝑡 運動はラグランジアンの時間積分を最小にするように決まる(※) ※「最小」ではなく「極小」の方が正確だが、以下では「最小」を使う 𝛿𝐼 = 0 𝑑𝑟𝑥 1 𝑑𝑡 = 𝑓1 (𝑟𝑥 1, 𝑟𝑦 1, 𝑟𝑧 1, ⋯ , 𝑟𝑥 𝑁, 𝑟𝑦 𝑁, 𝑟𝑧 𝑁, 𝑣𝑥 1, 𝑣𝑦 1, 𝑣𝑧 1, ⋯ , 𝑣𝑥 𝑁, 𝑣𝑦 𝑁, 𝑣𝑧 𝑁) 𝑑𝑣𝑧 𝑁 𝑑𝑡 = 𝑓6𝑁 (𝑟𝑥 1, 𝑟𝑦 1, 𝑟𝑧 1, ⋯ , 𝑟𝑥 𝑁, 𝑟𝑦 𝑁, 𝑟𝑧 𝑁, 𝑣𝑥 1, 𝑣𝑦 1, 𝑣𝑧 1, ⋯ , 𝑣𝑥 𝑁, 𝑣𝑦 𝑁, 𝑣𝑧 𝑁) ・・・
19. 19 85 歩きやすい 歩きにくい A B A地点からB地点まで最短時間で到達したい
20. 20 85 歩きやすい 歩きにくい A B A地点からB地点まで最短時間で到達したい 歩きづらさの総和が最小になる経路を選ぶ
21. 21 85 A B 屈折率が小さい 進みやすい 光路長が短い 屈折率が大きい 進みにくい 光路長が長い 光の経路は、光学的距離を最小(極小)にするように決まる 光は空中から水中に入る際に屈折する
22. 22 85 𝐼 = න 𝐶 𝑛 𝑥, 𝑦 𝑑𝑠 点(𝑥, 𝑦)における屈折率を𝑛(𝑥, 𝑦)として、光の経路を𝐶とすると 光路長を最小にするように経路が選ばれる(フェルマーの原理) 屈折率は、その場所における「光の進みづらさ」を表している 光は「進みづらさ」を最小にしたい この積分の意味は？ 積分を最小にする経路とは？
23. 23 85 𝐶 𝑥 𝑦 経路は一次元なので、一つのスカラーパラメータ𝑠で指定できる 𝑠 = 0 𝑠 = 1 (𝑥 𝑠 , 𝑦 𝑠 ) O
24. 24 85 屈折率𝑛(𝑥, 𝑦)は場所(𝑥, 𝑦)の関数 →場所(ベクトル)が入力、スカラーが出力 経路𝐶は、パラメータ𝑠の関数 →スカラーが入力、場所(ベクトル)が出力 𝐼 = න 𝐶 𝑛 𝑥, 𝑦 𝑑𝑠 作用積分𝐼は関数(経路)の汎関数 → 経路が入力、スカラーが出力 ここを決めたい ここを最小に するように これを決める方法が汎関数微分
25. 25 85 𝐼 = න 𝑡0 𝑡1 ℒ( ሶ 𝑟, 𝑟)𝑑𝑡 𝛿ℒ = 𝜕ℒ ሶ 𝜕𝑟 𝛿 ሶ 𝑟 + 𝜕ℒ 𝜕𝑟 𝛿𝑟 = − 𝑑 𝑑𝑡 𝜕ℒ ሶ 𝜕𝑟 𝛿𝑟 + 𝜕ℒ 𝜕𝑟 𝛿𝑟 = − 𝑑 𝑑𝑡 𝜕ℒ ሶ 𝜕𝑟 + 𝜕ℒ 𝜕𝑟 𝛿𝑟 汎関数微分 部分積分 𝛿𝑟でまとめる
26. 26 85 𝐼 = න 𝑡0 𝑡1 ℒ( ሶ 𝑟, 𝑟)𝑑𝑡 が最小 →𝛿𝐼 = 0 𝛿ℒ = − 𝑑 𝑑𝑡 𝜕ℒ ሶ 𝜕𝑟 + 𝜕ℒ 𝜕𝑟 𝛿𝑟 𝛿𝐼 = 0 → =0 − 𝑑 𝑑𝑡 𝜕ℒ ሶ 𝜕𝑟 + 𝜕ℒ 𝜕𝑟 = 0 以上から、 これをオイラー・ラグランジュ方程式とよぶ
27. 27 85 x,y,z方向を区別するのが面倒なので 𝑟𝑥 𝑖 = 𝑞3𝑖 , 𝑟𝑦 𝑖 = 𝑞3𝑖+1 , 𝑟𝑦 𝑖 = 𝑞3𝑖+2 と通し番号をつける − 𝑑 𝑑𝑡 𝜕ℒ 𝜕 ሶ 𝑟𝑖 + 𝜕ℒ 𝜕𝑟𝑖 = 0 𝛿𝐼 = 0 ひとつの多変数スカラー関数の変分から 必要な数の微分方程式が全て得られる ℒ(𝑞1 , ⋯ 𝑞3𝑁 , ሶ 𝑞1 , ⋯ , ሶ 𝑞3𝑁 )
28. 28 85 • ラグランジアンは位相空間における「進み づらさ」を表す関数 • 作用積分とは軌道に沿ったラグランジアン の線積分であり、作用とは進みづらさの総 和を意味する • 多粒子系の運動方程式は連立微分方程式と なるが、その全てが単一のスカラー関数で あるラグランジアンの変分から得られる(オ イラー・ラグランジュ方程式)
29. 29 85 ラグランジアンの第一引数をルジャンドル変換 ℒ( ሶ 𝑞, 𝑞) ℋ(𝑝, 𝑞) 𝑝 = 𝜕ℒ 𝜕 ሶ 𝑞 ℋ = 𝑝𝑞 − ℒ ラグランジアンから ハミルトニアンへ
30. 30 85 双対変換の一種 双対変換とは「何かを入れ替える」変換 二度行うともとに戻る 変換で情報は失われない 双対変換の例 正多面体の「点」と「面」の入れ替え 正四面体←→正四面体 正六面体←→正八面体 正四面体←→正八面体 正十二面体←→正二十面体 ※フーリエ変換なども双対変換
31. 31 85 ルジャンドル変換は(𝑥, 𝑦)空間から(𝑋, 𝑌)への双対変換 ルジャンドル変換には二つの表式がある 接線表式 面線表式 𝑦 = 𝑋𝑥 + 𝑌 𝑦 + 𝑌 = 𝑥𝑋
32. 32 85 O 𝑥 𝑦 𝑦 = 𝑋𝑥 + 𝑌 𝑥 𝑦 = 𝑦(𝑥) 関数𝑦 = 𝑦(𝑥)の各点に接線をひく 接線の傾きを𝑋、y切片を𝑌とする 𝑋 = 𝑑𝑦 𝑑𝑥 𝑌 = 𝑦 − 𝑥𝑋
33. 33 85 𝑋 = 𝑑𝑦 𝑑𝑥 𝑌 = 𝑦 − 𝑥𝑋 𝑥 = − 𝑑𝑌 𝑑𝑋 𝑦 = 𝑌 − 𝑥𝑋 逆変換を考えてやると ここに余計な負符号がつく 最終的に辻褄が合えばどうでも良いが、これを嫌って 𝑌 = 𝑥𝑋 − 𝑦 と定義すると、変換と逆変換が対称になる
34. 34 85 O 𝑋 = 𝑓 𝑥 𝑋 𝑔 𝑋 𝑥 𝑥 = 𝑔 𝑋 𝑋 𝑥 𝑓 𝑥 (𝑥, 𝑋)空間上の曲線を考える xを自由変数にとった場合 Xを自由変数にとった場合
35. 35 85 O 𝑋 𝑥 𝑋 𝑥 𝑦 = න 0 𝑥 𝑓 𝑥 𝑑𝑥 𝑌 = න 0 𝑋 𝑔 𝑋 𝑑𝑋 𝑦 + 𝑌 = 𝑥𝑋 長方形の面積 曲線の下側の面積 曲線の左側の面積 𝑋 = 𝑓 𝑥 𝑥 = 𝑔 𝑋
36. 36 85 𝑦 + 𝑌 = 𝑥𝑋 𝑌 = න 0 𝑋 𝑔 𝑋 𝑑𝑋 𝑋 = 𝑑𝑦 𝑑𝑥 𝑦 = න 0 𝑥 𝑓 𝑥 𝑑𝑥 𝑥 = 𝑑𝑌 𝑑𝑋 𝑋 = 𝑓 𝑥 𝑥 = 𝑔 𝑋 変換公式が自然に対称になる
37. 37 85 ルジャンドル変換は双対変換 双対変換は情報を保存する →ハミルトニアンへの変換で情報は増えない なぜハミルトニアンを考えるか？ 見通しが良くなるから − 𝑑 𝑑𝑡 𝜕ℒ ሶ 𝜕𝑞 + 𝜕ℒ 𝜕𝑞 = 0 𝑑𝑝 𝑑𝑡 = − 𝜕ℋ 𝜕𝑞 𝑑𝑞 𝑑𝑡 = 𝜕ℋ 𝜕𝑝 オイラー・ラグランジュ方程式 ハミルトンの運動方程式
38. 38 85 𝑑𝑝 𝑑𝑡 = − 𝜕ℋ 𝜕𝑞 𝑑𝑞 𝑑𝑡 = 𝜕ℋ 𝜕𝑝 ハミルトンの運動方程式 ℋの時間微分 𝑑ℋ 𝑑𝑡 = 𝜕ℋ 𝜕𝑝 ሶ 𝑝 + 𝜕ℋ 𝜕𝑞 ሶ 𝑞 = − 𝜕ℋ 𝜕𝑝 𝜕ℋ 𝜕𝑞 + 𝜕ℋ 𝜕𝑞 𝜕ℋ 𝜕𝑝 = 0 エネルギーの時間微分がゼロ→エネルギーが保存する この事実の幾何学的な意味を考える
39. 39 85 pとqをまとめて一つのベクトルで表現する Ԧ 𝑧 = 𝑝 𝑞 ハミルトニアンの勾配(gradient)を求める ∇ℋ = 𝜕𝑝 ℋ 𝜕𝑞 ℋ 運動方程式がベクトルの式で表現できる 𝑑 Ԧ 𝑧 𝑑𝑡 = ሶ 𝑝 ሶ 𝑞 = 0 −1 1 0 ∇ℋ
40. 40 85 ∇ℋ = 𝜕𝑝 ℋ 𝜕𝑞 ℋ p q スカラー場の勾配は、最も変化が大きい方向へのベクトル 勾配ベクトルは等高線と直交する エネルギーの等高線 この点での勾配
41. 41 85 𝑑 Ԧ 𝑧 𝑑𝑡 = ሶ 𝑝 ሶ 𝑞 = −𝜕𝑞 ℋ 𝜕𝑝 ℋ = 0 −1 1 0 ∇ℋ 勾配ベクトルを 90度回している p q ∇ℋ 𝑑 Ԧ 𝑧 𝑑𝑡 運動方向は、勾配ベクトルと直交する →等高線に沿って運動する →エネルギーが保存する
42. 42 85 𝑑𝑝 𝑑𝑡 = − 𝜕ℋ 𝜕𝑞 𝑑𝑞 𝑑𝑡 = 𝜕ℋ 𝜕𝑝 ハミルトンの運動方程式 ∇ 𝑑 Ԧ 𝑧 𝑑𝑡 = 𝜕 ሶ 𝑝 𝜕𝑝 + 𝜕 ሶ 𝑞 𝜕𝑞 = 0 速度場の発散(divergence)がゼロ ハミルトンベクトル場の発散がゼロの意味とは？ 𝑑 Ԧ 𝑧 𝑑𝑡 ハミルトニアンが作るベクトル場 →ハミルトンベクトル場
43. 43 85 𝑥 𝑥 𝑥 + ℎ 𝑣(𝑥 + ℎ) 𝑣(𝑥) 速度場𝑣(𝑥)の発散 入ってくる量 出ていく量 𝑣 𝑥 + ℎ − 𝑣 𝑥 ∼ 𝑑𝑣 𝑑𝑥 ℎ 出ていく量 入ってくる量 領域中に残る量
44. 44 85 p q O ∇ 𝑑 Ԧ 𝑧 𝑑𝑡 = 𝜕 ሶ 𝑝 𝜕𝑝 + 𝜕 ሶ 𝑞 𝜕𝑞 p方向の収支 q方向の収支 単位時間、単位面積あたりの収支 速度場の発散は、その点における量の非保存を意味する
45. 45 85 𝑑𝑝 𝑑𝑡 = − 𝜕ℋ 𝜕𝑞 𝑑𝑞 𝑑𝑡 = 𝜕ℋ 𝜕𝑝 ハミルトンの運動方程式 ∇ 𝑑 Ԧ 𝑧 𝑑𝑡 = 𝜕 ሶ 𝑝 𝜕𝑝 + 𝜕 ሶ 𝑞 𝜕𝑞 = 0 速度場の発散(divergence)がゼロ 速度場の発散がゼロ →各点で入ってくる量と出ていく量が等しい →流れに沿って密度が変化しない →ハミルトンベクトル場は非圧縮流
46. 46 85 O 非圧縮流：流れに沿って位相空間体積が変化しない
47. 47 85 • ハミルトンの運動方程式は、位相空間に流れ場 を作る(ハミルトンベクトル場) • ハミルトンベクトル場は、回転する流れを表現 している • ハミルトンベクトル場は、非圧縮流になってい る 時間発展とは、位相空間における回転
48. 48 85 一次元一自由度系 𝑝, 𝑞 における物理量𝐴を考える 点(𝑝(𝑡), 𝑞(𝑡))は運動方程式に従って時間発展する この世界の物理量は全て(𝑝, 𝑞)の関数として表現できる →𝐴は(𝑝, 𝑞)の関数 𝐴(𝑝 𝑡 , 𝑞 𝑡 ) 𝐴(𝑡) と略記 𝑝 𝑞 𝐴 𝐴(𝑡) (𝑝 𝑡 , 𝑞(𝑡)) 系の軌道 物理量の時間発展
49. 49 85 時刻𝑡 + ℎにおける値𝐴(𝑡 + ℎ)のテイラー展開を考える 𝐴 𝑡 + ℎ = 𝐴 𝑡 + ℎ 𝑑𝐴 𝑑𝑡 + ℎ2 2 𝑑2𝐴 𝑑𝑡2 + ⋯ = ෍ 𝑘=0 ℎ𝑘 𝑘! 𝑑𝑘 𝑑𝑡𝑘 𝐴(𝑡) = exp ℎ 𝑑 𝑑𝑡 𝐴(𝑡) ≡ 𝑈 ℎ 𝐴(𝑡) 𝑈 ℎ
50. 50 85 𝑈 ℎ 𝐴 𝑡 = 𝐴(𝑡 + ℎ) 𝐴(𝑡)に𝑈(ℎ)をかけると 時刻がℎだけ進む 𝑈 ℎ = exp ℎ 𝑑 𝑑𝑡 は時間発展演算子
51. 51 85 𝐴の時間微分を考える 𝐴は(𝑝, 𝑞)の関数であったから 𝑑𝐴 𝑑𝑡 = 𝜕𝐴 𝜕𝑝 ሶ 𝑝 + 𝜕𝐴 𝜕𝑞 ሶ 𝑞 = ሶ 𝑝 𝜕 𝜕𝑝 + ሶ 𝑞 𝜕 𝜕𝑞 𝐴 −𝑖𝐿 = −𝑖𝐿𝐴
52. 52 85 −𝑖𝐿𝐴 = 𝑑𝐴 𝑑𝑡 𝐴(𝑡)に−𝑖𝐿をかけると 時間微分が得られる −𝑖𝐿 = ሶ 𝑝 𝜕 𝜕𝑝 + ሶ 𝑞 𝜕 𝜕𝑞 は時間微分演算子 これをリュービル演算子(Liouville operator)と呼ぶ ※ 虚数単位をつけるのは、演算子をエルミートにするため ※ 物理量にかける場合は負符号をつける
53. 53 85 −𝑖𝐿𝐴 = 𝑑𝐴 𝑑𝑡 形式的に積分する 微分して自分自身が出てくるので指数関数 𝐴(𝑡 + ℎ) = exp −𝑖ℎ𝐿 𝐴(𝑡) 𝐴(𝑡 + ℎ) = 𝑈(ℎ)𝐴(𝑡) 時間発展演算子の定義 𝑈(ℎ) = exp −𝑖ℎ𝐿 時間発展演算子はリュービル演算子を指数関数の肩に載せたもの
54. 54 85 𝑑𝐴 𝑑𝑡 = −𝑖𝐿𝐴 任意の物理量𝐴について なので 𝑑 𝑑𝑡 = −𝑖𝐿 𝑈 ℎ = exp ℎ 𝑑 𝑑𝑡 微分演算子をリュービル演算子 で置き換えた
55. 55 85 一次元調和振動子の運動方程式 ሶ 𝑝 = −𝑞 ሶ 𝑞 = 𝑝 Ԧ 𝑧(𝑡) = 𝑝(𝑡) 𝑞(𝑡) とベクトル表示すると 𝑑 𝑑𝑡 Ԧ 𝑧 = 𝐿 Ԧ 𝑧 𝐿 = 0 −1 1 0 リュービル演算子の行列表示 ※計算の煩雑さをさけるため、負符号と虚数単位は省いてある
56. 56 85 𝑑 𝑑𝑡 Ԧ 𝑧 = 𝐿 Ԧ 𝑧 を形式的に解くと Ԧ 𝑧 𝑡 + ℎ = exp(ℎ𝐿) Ԧ 𝑧 𝑡 exp ℎ𝐿 = 𝐼 + ℎ𝐿 + ℎ2 2 𝐿2 + ⋯ + ℎ𝑘 𝑘! 𝐿𝑘 + ⋯ 𝐿 = 0 −1 1 0 𝐿2 = 1 0 0 −1 などから厳密に計算できる 𝑈 ℎ
57. 57 85 𝐿 = 0 −1 1 0 exp ℎ𝐿 = cos ℎ −sin ℎ sin ℎ cos ℎ 回転行列 回転行列の生成子 (微小回転) リュービル演算子は回転の「方向」を表現している 時間発展演算子は有限の回転を表現している
58. 58 85 時間を進める演算子が時間発展演算子 𝑈 ℎ 𝐴 𝑡 = 𝐴(𝑡 + ℎ) 時間で微分する演算子がリュービル演算子 −𝑖𝐿𝐴 𝑡 = 𝑑𝐴 𝑑𝑡 リュービル演算子を指数関数の肩に載せたものが 時間発展演算子 𝑈(ℎ) = exp −𝑖ℎ𝐿 時間発展とは回転であり、有限の回転を表現するのが 時間発展演算子で、無限小の回転を表現するのが リュービル演算子
59. 59 85 微分方程式が与えられている時、現時点での値から将来を予測したい 𝑑𝑥 𝑑𝑡 = 𝑓(𝑥) 𝑥(𝑡) 𝑥(𝑡 + ℎ) 方程式 初期値 未来の値 𝑥(𝑡 + ℎ) を 𝑥(𝑡) と 𝑓(𝑥) で近似する この積分が厳密評価できない 𝑥 𝑡 + ℎ = 𝑥 𝑡 + න 𝑡 𝑡+ℎ 𝑓 𝑥 𝑑𝑡 厳密な表式
60. 60 85 𝑥 𝑡 + ℎ = 𝑥 𝑡 + න 𝑡 𝑡+ℎ 𝑓 𝑥 𝑑𝑡 ∼ 𝑥 𝑡 + 𝑓 𝑥 ℎ 積分区間が短いと思って𝑓(𝑥)を定数で近似する 𝑥 𝑡 + ℎ = 𝑥 𝑡 + 𝑓 𝑥 𝑡 ℎ 知りたい量 知ってる量 未知の量𝑥(𝑡 + ℎ)を、既知の量𝑥 𝑡 と𝑓(𝑥(𝑡))だけで表現できた この数値積分法をオイラー法と呼ぶ
61. 61 85 𝑑𝑝 𝑑𝑡 = − 𝜕ℋ 𝜕𝑞 𝑑𝑞 𝑑𝑡 = 𝜕ℋ 𝜕𝑝 運動方程式 現在の値 𝑝 𝑡 , 𝑞(𝑡) 未来の値 𝑝 𝑡 + ℎ , 𝑞(𝑡 + ℎ) 運動方程式を使って、未来の値を現在の値の関数として表現したい
62. 62 85 𝑝 𝑡 + ℎ = 𝑝 𝑡 + ሶ 𝑝 𝑡 ℎ = 𝑝 𝑡 − 𝜕ℋ 𝑑𝑞 ℎ 𝑞 𝑡 + ℎ = 𝑞 𝑡 + ሶ 𝑞 𝑡 ℎ = 𝑞 𝑡 + 𝜕ℋ 𝑑𝑝 ℎ 運動方程式にオイラー法を適用する 現在の値 現在の微分係数 (傾き) t 現在の値 真の未来の値 近似した未来の値 すごくずれていきそうな気がする
63. 63 85 t = 0 h = 0.05 p = 1.0 q = 0.0 T = 1000 for _ in range(T): print(f"{t} {p} {q}") t = t + h dp = - q * h dq = p * h p += dp q += dq 調和振動子の運動方程式をオイラー法で解くコード 𝑝 𝑡 + ℎ = 𝑝 𝑡 − 𝑞 𝑡 ℎ 𝑞 𝑡 + ℎ = 𝑞 𝑡 + 𝑝 𝑡 ℎ
64. 64 85 オイラー法で計算した調和振動子の軌道 単位円を描くべき軌道が 徐々に広がっている →エネルギーが増えている
65. 65 85 𝑝(𝑡 + ℎ) = 𝑝(𝑡) − 𝑞(𝑡)ℎ 𝑞 𝑡 + ℎ = 𝑞 𝑡 + 𝑝(𝑡)ℎ 調和振動子にオイラー法を適用 行列表示 𝑝(𝑡 + ℎ) 𝑞(𝑡 + ℎ) = 1 −ℎ ℎ 1 𝑝(𝑡) 𝑞(𝑡) 𝑝(𝑡 + ℎ) 𝑞(𝑡 + ℎ) = cos ℎ −sin ℎ sin ℎ cos ℎ 𝑝(𝑡) 𝑞(𝑡) 厳密解 行列式= 1 + ℎ2 > 1 行列式= 1 位相空間の面積非保存がエネルギー非保存の原因
66. 66 85 オイラー法は、テイラー展開の一次まで正しいコード →高次の積分法を構築する 𝑞 𝑡 + ℎ = 𝑞 𝑡 + ሶ 𝑞 𝑡 ℎ + ሷ 𝑞 𝑡 ℎ2 2 ሶ 𝑝 = 𝑓(𝑞), ሶ 𝑞 = 𝑝 以下の運動方程式を考える 位置を二次までテイラー展開 = 𝑞 𝑡 + 𝑝(𝑡)ℎ + 𝑓(𝑡)ℎ2 2 ポテンシャル力
67. 67 85 𝑝 𝑡 + ℎ = 𝑝 𝑡 + ሶ 𝑝 𝑡 ℎ + ሷ 𝑝 𝑡 ℎ2 2 運動量も二次までテイラー展開 = 𝑝 𝑡 + 𝑓 𝑡 ℎ + ሶ 𝑓 𝑡 ℎ2 2 𝑑𝑓 𝑑𝑡 ∼ 𝑓 𝑡 + ℎ − 𝑓(𝑡) ℎ 力の時間微分を差分近似して代入 𝑝 𝑡 + ℎ = 𝑝 𝑡 + 𝑓 𝑡 + ℎ + 𝑓(𝑡) 2 ℎ 先に𝑞(𝑡 + ℎ)が求まっているので値がわかる
68. 68 85 ሶ 𝑝 = −𝑓(𝑞), ሶ 𝑞 = 𝑝 以下の運動方程式について 𝑝 𝑡 + ℎ = 𝑝 𝑡 + 𝑓 𝑡 + ℎ + 𝑓(𝑡) 2 ℎ 𝑞(𝑡 + ℎ) = 𝑞 𝑡 + 𝑝(𝑡)ℎ + 𝑓(𝑡)ℎ2 2 以下の二次の数値積分法を得た これをVelocity Verlet法と呼ぶ ※ 速度ベルレ法や、ベレの方法などとも呼ばれる
69. 69 85 t = 0 h = 0.05 p = 1.0 q = 0.0 T = 1000 for _ in range(T): print(f"{t} {p} {q}") t = t + h ft = -q q += p * h - q * h**2 * 0.5 ft2 = -q p += (ft2 + ft) * h * 0.5 調和振動子の運動方程式をVelocity Verlet法で解くコード 𝑝 𝑡 + ℎ = 𝑝 𝑡 + 𝑓 𝑡 + ℎ + 𝑓(𝑡) 2 ℎ 𝑞(𝑡 + ℎ) = 𝑞 𝑡 + 𝑝(𝑡)ℎ + 𝑓(𝑡)ℎ2 2
70. 70 85 Velocity Verlet法で計算した調和振動子の軌道 軌道が完全な円を描いている
71. 71 85 𝑝 𝑡 + ℎ = 𝑝 𝑡 + 𝑓 𝑡 + ℎ + 𝑓(𝑡) 2 ℎ 𝑞(𝑡 + ℎ) = 𝑞 𝑡 + 𝑝(𝑡)ℎ + 𝑓(𝑡)ℎ2 2 調和振動子にVelocity Verlet法を適用 𝑝(𝑡 + ℎ) 𝑞(𝑡 + ℎ) = 1 − ℎ2/2 −ℎ −ℎ + ℎ3/4 1 − ℎ2/2 𝑝(𝑡) 𝑞(𝑡) 行列表示 𝑝(𝑡 + ℎ) 𝑞(𝑡 + ℎ) = cos ℎ −sin ℎ sin ℎ cos ℎ 𝑝(𝑡) 𝑞(𝑡) 厳密解 行列式= 1 行列式= 1 Velocity Verlet法は時間発展演算子を近似しているが、 行列式は厳密に1
72. 72 85 𝑝 𝑡 + ℎ = 𝑃, 𝑞 𝑡 + ℎ = 𝑄 と書く 数値積分法とは(𝑝, 𝑞)から(𝑃, 𝑄)への写像を作るもの 厳密な時間発展ではエネルギーが保存する 𝑝2 2 + 𝑞2 2 = 𝑃2 2 + 𝑄2 2 Velocity Verletでは「ずれた」エネルギーが厳密に保存する 𝑝2 2 + 1 − ℎ2 4 𝑞2 2 = 𝑃2 2 + 1 − ℎ2 4 𝑄2 2 エネルギーが真の値のまわりを揺らぎながらも保存する
73. 73 85 𝑑 𝑑𝑡 = −𝑖𝐿 = − 𝜕ℋ 𝜕𝑞 𝜕 𝜕𝑝 + 𝜕ℋ 𝜕𝑝 𝜕 𝜕𝑞 リュービル演算子 時間発展演算子 𝑈(ℎ) = exp −𝑖ℎ𝐿 一般に厳密に計算できない リュービル演算子から近似した時間発展演算子を作りたい
74. 74 85 −𝑖𝐿 = − 𝜕ℋ 𝜕𝑞 𝜕 𝜕𝑝 + 𝜕ℋ 𝜕𝑝 𝜕 𝜕𝑞 = −𝑖𝐿𝑉 + −𝑖𝐿𝐾 −𝑖𝐿𝑉 −𝑖𝐿𝐾 リュービル演算子も二つに分解する ハミルトニアンが運動項とポテンシャル項に分離しているとする ℋ 𝑝, 𝑞 = 𝐾 𝑝 + 𝑉(𝑞) = − 𝜕𝑉 𝜕𝑞 𝜕 𝜕𝑝 + 𝜕𝐾 𝜕𝑝 𝜕 𝜕𝑞
75. 75 85 −𝑖𝐿𝑉 を𝑝に二回演算すると消える −𝑖𝐿𝑉 𝑝 = − 𝜕ℋ 𝜕𝑞 𝜕 𝜕𝑝 𝑝 = − 𝜕ℋ 𝜕𝑞 −𝑖𝐿𝑉 2𝑝 = − 𝜕ℋ 𝜕𝑞 𝜕 𝜕𝑝 − 𝜕ℋ 𝜕𝑞 = 0 ℋ 𝑝, 𝑞 = 𝐾 𝑝 + 𝑉(𝑞) ハミルトニアンが運動項とポテンシャル項に分かれている時 qだけの関数 pによる偏微分 −𝑖𝐿𝑉 を𝑞に演算すると消える −𝑖𝐿𝑉 𝑞 = − 𝜕ℋ 𝜕𝑞 𝜕 𝜕𝑝 𝑞 = 0 ※ − 𝑖𝐿𝐾 も同様
76. 76 85 𝑈𝐾 ℎ = exp −𝑖ℎ𝐿𝐾 = 1 + −𝑖ℎ𝐿𝐾 + −𝑖ℎ𝐿𝐾 2 2 + ⋯ ここがゼロなので これ以降全部ゼロ = 1 + −𝑖ℎ𝐿𝐾 𝑈𝑉 ℎ = exp −𝑖ℎ𝐿𝑉 = 1 + (−𝑖ℎ𝐿𝑉 ) 同様に 𝑈(ℎ) = exp −𝑖ℎ𝐿 は厳密に求められないが exp(−𝑖𝐿𝐾 )やexp(−𝑖𝐿𝑉 )は厳密に求めることができる
77. 77 85 一般に演算子𝐴と𝐵について exp 𝐴 + 𝐵 ≠ exp 𝐴 exp 𝐵 しかし、以下が成り立つ(Lie-Trotter公式) exp 𝐴 + 𝐵 = lim 𝑛→∞ exp 𝐴/𝑛 exp(𝐵/𝑛) 𝑛 有限のnで止めることで、以下の近似式を得る exp ℎ(𝐴 + 𝐵) = exp ℎ𝐴 exp ℎ𝐵 + 𝑂 ℎ2 exp ℎ(𝐴 + 𝐵) = exp ℎ𝐴/2 exp ℎ𝐵 exp ℎ𝐴/2 + 𝑂 ℎ3 一次の公式 二次の公式
78. 78 85 −𝑖𝐿 = −𝑖𝐿𝑉 + −𝑖𝐿𝐾 リュービル演算子を二つに分解する 求めたい時間発展演算子をLie-Trotter公式で分解する 𝑈 ℎ = exp −𝑖ℎ𝐿 = exp −𝑖ℎ𝐿𝑉 exp −𝑖ℎ𝐿𝐾 + 𝑂(ℎ2) = exp −𝑖ℎ𝐿𝑉 /2 exp −𝑖ℎ𝐿𝐾 exp −𝑖ℎ𝐿𝑉 /2 + 𝑂(ℎ3) = exp −𝑖ℎ(𝐿𝑉 + 𝐿𝐾 ) ෩ 𝑈1 (ℎ) ෩ 𝑈2 (ℎ) 一次近似された時間発展演算子 二次近似された時間発展演算子
79. 79 85 厳密な時間発展 𝑃 𝑄 = 𝑈 ℎ 𝑝 𝑞 = exp(−𝑖ℎ𝐿) 𝑝 𝑞 一次近似された時間発展 𝑃 𝑄 = ෩ 𝑈1 ℎ 𝑝 𝑞 = exp(−𝑖ℎ𝐿𝐾 )exp(−𝑖ℎ𝐿𝑉 ) 𝑝 𝑞 最初にこれを演算して 次にこれを演算する このように、指数分解公式で作られた時間発展演算子を使う 数値積分法をシンプレクティック積分法と呼ぶ
80. 80 85 𝑈𝑉 (ℎ)を𝑝に演算してみる 𝑈𝑉 ℎ 𝑝 = exp −𝑖ℎ𝐿𝑉 𝑝 = (1 − 𝑖ℎ𝐿𝑉 )𝑝 = 1 − ℎ 𝜕ℋ 𝜕𝑞 𝜕 𝜕𝑝 𝑝 = 𝑝 − ℎ 𝜕𝑉 𝜕𝑞 それぞれ、お互いを無視して時間ℎだけ運動 𝑈𝐾 (ℎ)を𝑞に演算してみる 𝑈𝐾 ℎ 𝑞 = exp −𝑖ℎ𝐿𝐾 𝑞 = 1 − 𝑖ℎ𝐿𝐾 𝑞 = 1 + ℎ 𝜕ℋ 𝜕𝑝 𝜕 𝜕𝑞 𝑞 = 𝑞 + ℎ 𝜕𝐾 𝜕𝑞 𝑝の時間をℎだけ進める 𝑞の時間をℎだけ進める
81. 81 85 調和振動子の場合 𝑈𝑉 ℎ 𝑝 𝑞 = 𝑝 − ℎ𝑞 𝑞 = 1 −ℎ 0 1 𝑝 𝑞 𝑈𝐾 ℎ 𝑝 𝑞 = 𝑝 𝑞 + ℎ𝑝 = 1 0 ℎ 1 𝑝 𝑞 𝑝の時間をℎだけ進める演算子 𝑞の時間をℎだけ進める演算子 二つともかけると 𝑈𝐾 ℎ 𝑈𝑉 ℎ 𝑝 𝑞 = 1 −ℎ ℎ 1 − ℎ2 𝑝 𝑞 行列式が1
82. 82 85 𝑝, 𝑞 から(𝑃, 𝑄)への変換 𝑃 𝑄 = 𝑈 𝑝 𝑞 において 𝑑𝑃𝑑𝑄 = 𝑑𝑝𝑑𝑞 が満たされる時、この変換をシンプレクティック変換と呼ぶ 𝑑𝑃𝑑𝑄 = 𝜕𝑝 𝑃 𝜕𝑞 𝑃 𝜕𝑝 𝑄 𝜕𝑞 𝑄 𝑑𝑝𝑑𝑞 であるから シンプレクティック変換とは、変換のヤコビアンが1であること ※ 時間発展演算子を行列表示したときの行列式が1
83. 83 85 ෩ 𝑈2 (ℎ) = exp −𝑖ℎ𝐿𝑉 /2 exp −𝑖ℎ𝐿𝐾 exp −𝑖ℎ𝐿𝑉 /2 二次の指数分解公式 (1) (2) (3) (1) 𝑝 𝑡 + ℎ/2 = 𝑝(𝑡) + 𝑓 𝑡 2 (2) 𝑄 = 𝑞(𝑡) + 𝑝 𝑡 + ℎ/2 ℎ (3) 𝑃 = 𝑝(𝑡 + ℎ/2) + 𝑓 𝑡 + ℎ ℎ 現在の力がそのまま続いたとしてℎ/2だけ時間を進める (1)で得た速度がそのまま続いたとしてℎだけ時間をすすめる (2)の力がそのまま続いたとしてℎ/2だけ時間を進める
84. 84 85 𝑝 𝑡 + ℎ/2 = 𝑝(𝑡) + 𝑓 𝑡 2 ℎ 𝑄 = 𝑞(𝑡) + 𝑝 𝑡 + ℎ/2 ℎ 𝑃 = 𝑝(𝑡 + ℎ/2) + 𝑓 𝑡 + ℎ ℎ 𝑝(𝑡 + ℎ/2)を消去して整理するとVelocity Verlet公式を得る 𝑃 = 𝑝 𝑡 + 𝑓 𝑡 + ℎ + 𝑓(𝑡) 2 ℎ 𝑄 = 𝑞 𝑡 + 𝑝(𝑡)ℎ + 𝑓(𝑡)ℎ2 2 Velocity Verlet法は、シンプレクティック積分
85. 85 85 • 数値積分とは、現在時刻𝑡の点(𝑝, 𝑞)から、時刻𝑡 + ℎの 点(𝑃, 𝑄)を得る写像である • 写像を表す変換のヤコビアンが1であるとき、位相空 間の体積が保存される。このような写像をシンプレク ティック写像と呼ぶ • 変換がシンプレクティック写像であるような数値積分 法をシンプレクティック積分と呼ぶ • シンプレクティック積分は、Lie-Trotterの指数分解公 式から構築できる (𝑝, 𝑞) (𝑃, 𝑄) 𝑑𝑃𝑑𝑄 = 𝑑𝑝𝑑𝑞