Upgrade to Pro — share decks privately, control downloads, hide ads and more …

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

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

シミュレーション工学

kaityo256
PRO

May 08, 2023
Tweet

More Decks by kaityo256

Other Decks in Education

Transcript

  1. 1
    85
    分子動力学法(1)理論的背景と数値積分
    シミュレーション工学
    慶應義塾大学大学院理工学研究科基礎理工学専攻物理情報専修
    渡辺宙志

    View Slide

  2. 2
    85
    • 解析力学の復習
    • 分子動力学法の時間発展アルゴリズム
    本講義の目的
    分子動力学法とは
    • 原子や分子にかかる「力」を計算し、位置
    を更新していく数値計算手法
    • 背景にある理論は解析力学

    View Slide

  3. 3
    85
    𝑚𝑎 = 𝐹
    物体の加速度は力に比例する
    物体の動きにくさは質量に比例する

    View Slide

  4. 4
    85
    𝑎 =
    𝑑𝑣
    𝑑𝑡
    𝑣 =
    𝑑𝑟
    𝑑𝑡
    加速度とは速度の時間変化率
    速度とは位置の時間変化率
    𝑑2𝑟
    𝑑𝑡2
    =
    𝐹
    𝑚
    運動方程式は
    位置の二階微分方程式

    View Slide

  5. 5
    85
    O
    𝑟
    自然長の位置を原点にとる
    質量とバネ定数は1とする
    𝑑2𝑟
    𝑑𝑡2
    = −𝑟

    View Slide

  6. 6
    85
    𝑟 = 0 物体が原点にいる
    これだけでは系の状態が定まらない
    止まった状態かもしれない
    運動中にちょうど原点を
    通過したところかもしれない

    View Slide

  7. 7
    85
    𝑟 = 0,𝑣 = 1 物体が原点にいて、速度が1
    𝑑2𝑟
    𝑑𝑡2
    = −𝑟
    二階の常微分方程式で初期条件を二つ与えたので、
    この系の運動は完全に定まる

    View Slide

  8. 8
    85
    (𝑣, 𝑟) = (1,0)
    𝑣
    𝑟
    O
    この系の状態は(𝑣, 𝑟)空間の点を指定することで一意に定まる
    この空間を位相空間(Phase Space)と呼ぶ

    View Slide

  9. 9
    85
    系の状態は(𝑣, 𝑟)空間で指定できるのだから
    運動方程式も(𝑣, 𝑟)で書きたい
    𝑑2𝑟
    𝑑𝑡2
    = −𝑟
    𝑑𝑣
    𝑑𝑡
    = −𝑟
    二階の常微分方程式を一階の連立常微分方程式に
    𝑑𝑟
    𝑑𝑡
    = 𝑣

    View Slide

  10. 10
    85
    変数の組(𝑣, 𝑟)の時間微分( ሶ
    𝑣, ሶ
    𝑟)が
    (𝑣, 𝑟)の関数として書けている
    𝑑𝑣
    𝑑𝑡
    = 𝑓𝑣
    (𝑟, 𝑣)
    𝑑𝑟
    𝑑𝑡
    = 𝑓𝑟
    (𝑟, 𝑣) 𝑓𝑟
    𝑟, 𝑣 = 𝑣
    𝑓𝑣
    𝑟, 𝑣 = −𝑟
    このような系を力学系(dynamical system)と呼ぶ
    一般的に書くと
    ※運動方程式は力学系の一種
    𝑑𝑣
    𝑑𝑡
    = −𝑟
    𝑑𝑟
    𝑑𝑡
    = 𝑣

    View Slide

  11. 11
    85
    (𝑣, 𝑟) = (1,0)
    𝑟
    𝑣
    O
    現在の状態(原点で右向きに運動)
    時間微分(右向きに運動し、加速度0)
    ( ሶ
    𝑣, ሶ
    𝑟) = (0,1) ( ሶ
    𝑣, ሶ
    𝑟) = (−𝑟, 𝑣)

    View Slide

  12. 12
    85
    現在の状態(バネが伸び切って停止)
    (𝑣, 𝑟) = (0,1)
    時間微分(左向きに加速し、速度0)
    ( ሶ
    𝑣, ሶ
    𝑟) = (−1,0)
    𝑟
    𝑣
    O
    ( ሶ
    𝑣, ሶ
    𝑟) = (−𝑟, 𝑣)

    View Slide

  13. 13
    85
    𝑟
    𝑣
    O
    今の手続きを様々な点で繰り返すと・・・
    運動方程式は位相空間にベクトル場を定義する
    回転する速度場

    View Slide

  14. 14
    85
    𝑟
    𝑣
    O
    (𝑣, 𝑟) = (1,0)
    初期条件:位相空間にトレーサーを置くこと
    方程式の解:トレーサーの軌跡

    View Slide

  15. 15
    85
    • 運動方程式に従う系の状態を一意に決め
    ることができる空間を位相空間と呼ぶ
    • 運動方程式とは、位相空間にベクトル場
    (速度場)を定義するものである
    • 系の状態は、位相空間上の一点で表現で
    きる
    • 系の運動は、位相空間上の軌跡で表現さ
    れ、運動方程式が定義した速度場に従っ
    て動く点の軌跡である

    View Slide

  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𝑁

    View Slide

  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𝑁

    View Slide

  18. 18
    85
    𝐼 = න
    𝑡0
    𝑡1
    ℒ( ሶ
    𝑟, 𝑟)𝑑𝑡
    運動はラグランジアンの時間積分を最小にするように決まる(※)
    ※「最小」ではなく「極小」の方が正確だが、以下では「最小」を使う
    𝛿𝐼 = 0
    𝑑𝑟𝑥
    1
    𝑑𝑡
    = 𝑓1
    (𝑟𝑥
    1, 𝑟𝑦
    1, 𝑟𝑧
    1, ⋯ , 𝑟𝑥
    𝑁, 𝑟𝑦
    𝑁, 𝑟𝑧
    𝑁, 𝑣𝑥
    1, 𝑣𝑦
    1, 𝑣𝑧
    1, ⋯ , 𝑣𝑥
    𝑁, 𝑣𝑦
    𝑁, 𝑣𝑧
    𝑁)
    𝑑𝑣𝑧
    𝑁
    𝑑𝑡
    = 𝑓6𝑁
    (𝑟𝑥
    1, 𝑟𝑦
    1, 𝑟𝑧
    1, ⋯ , 𝑟𝑥
    𝑁, 𝑟𝑦
    𝑁, 𝑟𝑧
    𝑁, 𝑣𝑥
    1, 𝑣𝑦
    1, 𝑣𝑧
    1, ⋯ , 𝑣𝑥
    𝑁, 𝑣𝑦
    𝑁, 𝑣𝑧
    𝑁)
    ・・・

    View Slide

  19. 19
    85
    歩きやすい
    歩きにくい
    A
    B
    A地点からB地点まで最短時間で到達したい

    View Slide

  20. 20
    85
    歩きやすい
    歩きにくい
    A
    B
    A地点からB地点まで最短時間で到達したい
    歩きづらさの総和が最小になる経路を選ぶ

    View Slide

  21. 21
    85
    A
    B
    屈折率が小さい
    進みやすい
    光路長が短い
    屈折率が大きい
    進みにくい
    光路長が長い
    光の経路は、光学的距離を最小(極小)にするように決まる
    光は空中から水中に入る際に屈折する

    View Slide

  22. 22
    85
    𝐼 = න
    𝐶
    𝑛 𝑥, 𝑦 𝑑𝑠
    点(𝑥, 𝑦)における屈折率を𝑛(𝑥, 𝑦)として、光の経路を𝐶とすると
    光路長を最小にするように経路が選ばれる(フェルマーの原理)
    屈折率は、その場所における「光の進みづらさ」を表している
    光は「進みづらさ」を最小にしたい
    この積分の意味は?
    積分を最小にする経路とは?

    View Slide

  23. 23
    85
    𝐶
    𝑥
    𝑦
    経路は一次元なので、一つのスカラーパラメータ𝑠で指定できる
    𝑠 = 0
    𝑠 = 1
    (𝑥 𝑠 , 𝑦 𝑠 )
    O

    View Slide

  24. 24
    85
    屈折率𝑛(𝑥, 𝑦)は場所(𝑥, 𝑦)の関数
    →場所(ベクトル)が入力、スカラーが出力
    経路𝐶は、パラメータ𝑠の関数
    →スカラーが入力、場所(ベクトル)が出力
    𝐼 = න
    𝐶
    𝑛 𝑥, 𝑦 𝑑𝑠
    作用積分𝐼は関数(経路)の汎関数
    → 経路が入力、スカラーが出力
    ここを決めたい
    ここを最小に
    するように
    これを決める方法が汎関数微分

    View Slide

  25. 25
    85
    𝐼 = න
    𝑡0
    𝑡1
    ℒ( ሶ
    𝑟, 𝑟)𝑑𝑡
    𝛿ℒ =
    𝜕ℒ

    𝜕𝑟
    𝛿 ሶ
    𝑟 +
    𝜕ℒ
    𝜕𝑟
    𝛿𝑟
    = −
    𝑑
    𝑑𝑡
    𝜕ℒ

    𝜕𝑟
    𝛿𝑟 +
    𝜕ℒ
    𝜕𝑟
    𝛿𝑟
    = −
    𝑑
    𝑑𝑡
    𝜕ℒ

    𝜕𝑟
    +
    𝜕ℒ
    𝜕𝑟
    𝛿𝑟
    汎関数微分
    部分積分
    𝛿𝑟でまとめる

    View Slide

  26. 26
    85
    𝐼 = න
    𝑡0
    𝑡1
    ℒ( ሶ
    𝑟, 𝑟)𝑑𝑡 が最小 →𝛿𝐼 = 0
    𝛿ℒ = −
    𝑑
    𝑑𝑡
    𝜕ℒ

    𝜕𝑟
    +
    𝜕ℒ
    𝜕𝑟
    𝛿𝑟
    𝛿𝐼 = 0 →
    =0

    𝑑
    𝑑𝑡
    𝜕ℒ

    𝜕𝑟
    +
    𝜕ℒ
    𝜕𝑟
    = 0
    以上から、
    これをオイラー・ラグランジュ方程式とよぶ

    View Slide

  27. 27
    85
    x,y,z方向を区別するのが面倒なので
    𝑟𝑥
    𝑖 = 𝑞3𝑖
    , 𝑟𝑦
    𝑖 = 𝑞3𝑖+1
    , 𝑟𝑦
    𝑖 = 𝑞3𝑖+2
    と通し番号をつける

    𝑑
    𝑑𝑡
    𝜕ℒ
    𝜕 ሶ
    𝑟𝑖
    +
    𝜕ℒ
    𝜕𝑟𝑖
    = 0
    𝛿𝐼 = 0
    ひとつの多変数スカラー関数の変分から
    必要な数の微分方程式が全て得られる
    ℒ(𝑞1
    , ⋯ 𝑞3𝑁
    , ሶ
    𝑞1
    , ⋯ , ሶ
    𝑞3𝑁
    )

    View Slide

  28. 28
    85
    • ラグランジアンは位相空間における「進み
    づらさ」を表す関数
    • 作用積分とは軌道に沿ったラグランジアン
    の線積分であり、作用とは進みづらさの総
    和を意味する
    • 多粒子系の運動方程式は連立微分方程式と
    なるが、その全てが単一のスカラー関数で
    あるラグランジアンの変分から得られる(オ
    イラー・ラグランジュ方程式)

    View Slide

  29. 29
    85
    ラグランジアンの第一引数をルジャンドル変換
    ℒ( ሶ
    𝑞, 𝑞) ℋ(𝑝, 𝑞)
    𝑝 =
    𝜕ℒ
    𝜕 ሶ
    𝑞
    ℋ = 𝑝𝑞 − ℒ
    ラグランジアンから ハミルトニアンへ

    View Slide

  30. 30
    85
    双対変換の一種
    双対変換とは「何かを入れ替える」変換
    二度行うともとに戻る
    変換で情報は失われない
    双対変換の例
    正多面体の「点」と「面」の入れ替え
    正四面体←→正四面体
    正六面体←→正八面体
    正四面体←→正八面体
    正十二面体←→正二十面体
    ※フーリエ変換なども双対変換

    View Slide

  31. 31
    85
    ルジャンドル変換は(𝑥, 𝑦)空間から(𝑋, 𝑌)への双対変換
    ルジャンドル変換には二つの表式がある
    接線表式
    面線表式
    𝑦 = 𝑋𝑥 + 𝑌
    𝑦 + 𝑌 = 𝑥𝑋

    View Slide

  32. 32
    85
    O
    𝑥
    𝑦
    𝑦 = 𝑋𝑥 + 𝑌
    𝑥
    𝑦 = 𝑦(𝑥)
    関数𝑦 = 𝑦(𝑥)の各点に接線をひく
    接線の傾きを𝑋、y切片を𝑌とする
    𝑋 =
    𝑑𝑦
    𝑑𝑥
    𝑌 = 𝑦 − 𝑥𝑋

    View Slide

  33. 33
    85
    𝑋 =
    𝑑𝑦
    𝑑𝑥
    𝑌 = 𝑦 − 𝑥𝑋
    𝑥 = −
    𝑑𝑌
    𝑑𝑋
    𝑦 = 𝑌 − 𝑥𝑋
    逆変換を考えてやると
    ここに余計な負符号がつく
    最終的に辻褄が合えばどうでも良いが、これを嫌って
    𝑌 = 𝑥𝑋 − 𝑦
    と定義すると、変換と逆変換が対称になる

    View Slide

  34. 34
    85
    O
    𝑋 = 𝑓 𝑥
    𝑋
    𝑔 𝑋
    𝑥
    𝑥 = 𝑔 𝑋
    𝑋
    𝑥
    𝑓 𝑥
    (𝑥, 𝑋)空間上の曲線を考える
    xを自由変数にとった場合
    Xを自由変数にとった場合

    View Slide

  35. 35
    85
    O
    𝑋
    𝑥
    𝑋
    𝑥
    𝑦 = න
    0
    𝑥
    𝑓 𝑥 𝑑𝑥
    𝑌 = න
    0
    𝑋
    𝑔 𝑋 𝑑𝑋
    𝑦 + 𝑌 = 𝑥𝑋
    長方形の面積
    曲線の下側の面積 曲線の左側の面積
    𝑋 = 𝑓 𝑥
    𝑥 = 𝑔 𝑋

    View Slide

  36. 36
    85
    𝑦 + 𝑌 = 𝑥𝑋
    𝑌 = න
    0
    𝑋
    𝑔 𝑋 𝑑𝑋
    𝑋 =
    𝑑𝑦
    𝑑𝑥
    𝑦 = න
    0
    𝑥
    𝑓 𝑥 𝑑𝑥
    𝑥 =
    𝑑𝑌
    𝑑𝑋
    𝑋 = 𝑓 𝑥
    𝑥 = 𝑔 𝑋
    変換公式が自然に対称になる

    View Slide

  37. 37
    85
    ルジャンドル変換は双対変換
    双対変換は情報を保存する
    →ハミルトニアンへの変換で情報は増えない
    なぜハミルトニアンを考えるか?
    見通しが良くなるから

    𝑑
    𝑑𝑡
    𝜕ℒ

    𝜕𝑞
    +
    𝜕ℒ
    𝜕𝑞
    = 0
    𝑑𝑝
    𝑑𝑡
    = −
    𝜕ℋ
    𝜕𝑞
    𝑑𝑞
    𝑑𝑡
    =
    𝜕ℋ
    𝜕𝑝
    オイラー・ラグランジュ方程式 ハミルトンの運動方程式

    View Slide

  38. 38
    85
    𝑑𝑝
    𝑑𝑡
    = −
    𝜕ℋ
    𝜕𝑞
    𝑑𝑞
    𝑑𝑡
    =
    𝜕ℋ
    𝜕𝑝
    ハミルトンの運動方程式
    ℋの時間微分
    𝑑ℋ
    𝑑𝑡
    =
    𝜕ℋ
    𝜕𝑝

    𝑝 +
    𝜕ℋ
    𝜕𝑞

    𝑞
    = −
    𝜕ℋ
    𝜕𝑝
    𝜕ℋ
    𝜕𝑞
    +
    𝜕ℋ
    𝜕𝑞
    𝜕ℋ
    𝜕𝑝
    = 0
    エネルギーの時間微分がゼロ→エネルギーが保存する
    この事実の幾何学的な意味を考える

    View Slide

  39. 39
    85
    pとqをまとめて一つのベクトルで表現する
    Ԧ
    𝑧 =
    𝑝
    𝑞
    ハミルトニアンの勾配(gradient)を求める
    ∇ℋ =
    𝜕𝑝

    𝜕𝑞

    運動方程式がベクトルの式で表現できる
    𝑑 Ԧ
    𝑧
    𝑑𝑡
    =

    𝑝

    𝑞
    =
    0 −1
    1 0
    ∇ℋ

    View Slide

  40. 40
    85
    ∇ℋ =
    𝜕𝑝

    𝜕𝑞

    p
    q
    スカラー場の勾配は、最も変化が大きい方向へのベクトル
    勾配ベクトルは等高線と直交する
    エネルギーの等高線
    この点での勾配

    View Slide

  41. 41
    85
    𝑑 Ԧ
    𝑧
    𝑑𝑡
    =

    𝑝

    𝑞
    =
    −𝜕𝑞

    𝜕𝑝

    =
    0 −1
    1 0
    ∇ℋ
    勾配ベクトルを
    90度回している
    p
    q
    ∇ℋ
    𝑑 Ԧ
    𝑧
    𝑑𝑡
    運動方向は、勾配ベクトルと直交する
    →等高線に沿って運動する
    →エネルギーが保存する

    View Slide

  42. 42
    85
    𝑑𝑝
    𝑑𝑡
    = −
    𝜕ℋ
    𝜕𝑞
    𝑑𝑞
    𝑑𝑡
    =
    𝜕ℋ
    𝜕𝑝
    ハミルトンの運動方程式

    𝑑 Ԧ
    𝑧
    𝑑𝑡
    =
    𝜕 ሶ
    𝑝
    𝜕𝑝
    +
    𝜕 ሶ
    𝑞
    𝜕𝑞
    = 0
    速度場の発散(divergence)がゼロ
    ハミルトンベクトル場の発散がゼロの意味とは?
    𝑑 Ԧ
    𝑧
    𝑑𝑡
    ハミルトニアンが作るベクトル場
    →ハミルトンベクトル場

    View Slide

  43. 43
    85
    𝑥
    𝑥 𝑥 + ℎ
    𝑣(𝑥 + ℎ)
    𝑣(𝑥)
    速度場𝑣(𝑥)の発散
    入ってくる量 出ていく量
    𝑣 𝑥 + ℎ − 𝑣 𝑥 ∼
    𝑑𝑣
    𝑑𝑥

    出ていく量 入ってくる量
    領域中に残る量

    View Slide

  44. 44
    85
    p
    q
    O

    𝑑 Ԧ
    𝑧
    𝑑𝑡
    =
    𝜕 ሶ
    𝑝
    𝜕𝑝
    +
    𝜕 ሶ
    𝑞
    𝜕𝑞
    p方向の収支
    q方向の収支
    単位時間、単位面積あたりの収支
    速度場の発散は、その点における量の非保存を意味する

    View Slide

  45. 45
    85
    𝑑𝑝
    𝑑𝑡
    = −
    𝜕ℋ
    𝜕𝑞
    𝑑𝑞
    𝑑𝑡
    =
    𝜕ℋ
    𝜕𝑝
    ハミルトンの運動方程式

    𝑑 Ԧ
    𝑧
    𝑑𝑡
    =
    𝜕 ሶ
    𝑝
    𝜕𝑝
    +
    𝜕 ሶ
    𝑞
    𝜕𝑞
    = 0
    速度場の発散(divergence)がゼロ
    速度場の発散がゼロ
    →各点で入ってくる量と出ていく量が等しい
    →流れに沿って密度が変化しない
    →ハミルトンベクトル場は非圧縮流

    View Slide

  46. 46
    85
    O
    非圧縮流:流れに沿って位相空間体積が変化しない

    View Slide

  47. 47
    85
    • ハミルトンの運動方程式は、位相空間に流れ場
    を作る(ハミルトンベクトル場)
    • ハミルトンベクトル場は、回転する流れを表現
    している
    • ハミルトンベクトル場は、非圧縮流になってい

    時間発展とは、位相空間における回転

    View Slide

  48. 48
    85
    一次元一自由度系 𝑝, 𝑞 における物理量𝐴を考える
    点(𝑝(𝑡), 𝑞(𝑡))は運動方程式に従って時間発展する
    この世界の物理量は全て(𝑝, 𝑞)の関数として表現できる
    →𝐴は(𝑝, 𝑞)の関数
    𝐴(𝑝 𝑡 , 𝑞 𝑡 ) 𝐴(𝑡) と略記
    𝑝
    𝑞
    𝐴 𝐴(𝑡)
    (𝑝 𝑡 , 𝑞(𝑡))
    系の軌道
    物理量の時間発展

    View Slide

  49. 49
    85
    時刻𝑡 + ℎにおける値𝐴(𝑡 + ℎ)のテイラー展開を考える
    𝐴 𝑡 + ℎ = 𝐴 𝑡 + ℎ
    𝑑𝐴
    𝑑𝑡
    +
    ℎ2
    2
    𝑑2𝐴
    𝑑𝑡2
    + ⋯
    = ෍
    𝑘=0
    ℎ𝑘
    𝑘!
    𝑑𝑘
    𝑑𝑡𝑘
    𝐴(𝑡)
    = exp ℎ
    𝑑
    𝑑𝑡
    𝐴(𝑡)
    ≡ 𝑈 ℎ 𝐴(𝑡)
    𝑈 ℎ

    View Slide

  50. 50
    85
    𝑈 ℎ 𝐴 𝑡 = 𝐴(𝑡 + ℎ)
    𝐴(𝑡)に𝑈(ℎ)をかけると 時刻がℎだけ進む
    𝑈 ℎ = exp ℎ
    𝑑
    𝑑𝑡
    は時間発展演算子

    View Slide

  51. 51
    85
    𝐴の時間微分を考える
    𝐴は(𝑝, 𝑞)の関数であったから
    𝑑𝐴
    𝑑𝑡
    =
    𝜕𝐴
    𝜕𝑝

    𝑝 +
    𝜕𝐴
    𝜕𝑞

    𝑞
    = ሶ
    𝑝
    𝜕
    𝜕𝑝
    + ሶ
    𝑞
    𝜕
    𝜕𝑞
    𝐴
    −𝑖𝐿
    = −𝑖𝐿𝐴

    View Slide

  52. 52
    85
    −𝑖𝐿𝐴 =
    𝑑𝐴
    𝑑𝑡
    𝐴(𝑡)に−𝑖𝐿をかけると 時間微分が得られる
    −𝑖𝐿 = ሶ
    𝑝
    𝜕
    𝜕𝑝
    + ሶ
    𝑞
    𝜕
    𝜕𝑞
    は時間微分演算子
    これをリュービル演算子(Liouville operator)と呼ぶ
    ※ 虚数単位をつけるのは、演算子をエルミートにするため
    ※ 物理量にかける場合は負符号をつける

    View Slide

  53. 53
    85
    −𝑖𝐿𝐴 =
    𝑑𝐴
    𝑑𝑡
    形式的に積分する
    微分して自分自身が出てくるので指数関数
    𝐴(𝑡 + ℎ) = exp −𝑖ℎ𝐿 𝐴(𝑡)
    𝐴(𝑡 + ℎ) = 𝑈(ℎ)𝐴(𝑡) 時間発展演算子の定義
    𝑈(ℎ) = exp −𝑖ℎ𝐿
    時間発展演算子はリュービル演算子を指数関数の肩に載せたもの

    View Slide

  54. 54
    85
    𝑑𝐴
    𝑑𝑡
    = −𝑖𝐿𝐴
    任意の物理量𝐴について なので
    𝑑
    𝑑𝑡
    = −𝑖𝐿
    𝑈 ℎ = exp ℎ
    𝑑
    𝑑𝑡
    微分演算子をリュービル演算子
    で置き換えた

    View Slide

  55. 55
    85
    一次元調和振動子の運動方程式

    𝑝 = −𝑞

    𝑞 = 𝑝
    Ԧ
    𝑧(𝑡) =
    𝑝(𝑡)
    𝑞(𝑡)
    とベクトル表示すると
    𝑑
    𝑑𝑡
    Ԧ
    𝑧 = 𝐿 Ԧ
    𝑧 𝐿 =
    0 −1
    1 0
    リュービル演算子の行列表示
    ※計算の煩雑さをさけるため、負符号と虚数単位は省いてある

    View Slide

  56. 56
    85
    𝑑
    𝑑𝑡
    Ԧ
    𝑧 = 𝐿 Ԧ
    𝑧 を形式的に解くと
    Ԧ
    𝑧 𝑡 + ℎ = exp(ℎ𝐿) Ԧ
    𝑧 𝑡
    exp ℎ𝐿 = 𝐼 + ℎ𝐿 +
    ℎ2
    2
    𝐿2 + ⋯ +
    ℎ𝑘
    𝑘!
    𝐿𝑘 + ⋯
    𝐿 =
    0 −1
    1 0
    𝐿2 =
    1 0
    0 −1 などから厳密に計算できる
    𝑈 ℎ

    View Slide

  57. 57
    85
    𝐿 =
    0 −1
    1 0
    exp ℎ𝐿 =
    cos ℎ −sin ℎ
    sin ℎ cos ℎ
    回転行列
    回転行列の生成子
    (微小回転)
    リュービル演算子は回転の「方向」を表現している
    時間発展演算子は有限の回転を表現している

    View Slide

  58. 58
    85
    時間を進める演算子が時間発展演算子
    𝑈 ℎ 𝐴 𝑡 = 𝐴(𝑡 + ℎ)
    時間で微分する演算子がリュービル演算子
    −𝑖𝐿𝐴 𝑡 =
    𝑑𝐴
    𝑑𝑡
    リュービル演算子を指数関数の肩に載せたものが
    時間発展演算子
    𝑈(ℎ) = exp −𝑖ℎ𝐿
    時間発展とは回転であり、有限の回転を表現するのが
    時間発展演算子で、無限小の回転を表現するのが
    リュービル演算子

    View Slide

  59. 59
    85
    微分方程式が与えられている時、現時点での値から将来を予測したい
    𝑑𝑥
    𝑑𝑡
    = 𝑓(𝑥)
    𝑥(𝑡)
    𝑥(𝑡 + ℎ)
    方程式
    初期値
    未来の値
    𝑥(𝑡 + ℎ) を 𝑥(𝑡) と 𝑓(𝑥) で近似する
    この積分が厳密評価できない
    𝑥 𝑡 + ℎ = 𝑥 𝑡 + න
    𝑡
    𝑡+ℎ
    𝑓 𝑥 𝑑𝑡
    厳密な表式

    View Slide

  60. 60
    85
    𝑥 𝑡 + ℎ = 𝑥 𝑡 + න
    𝑡
    𝑡+ℎ
    𝑓 𝑥 𝑑𝑡
    ∼ 𝑥 𝑡 + 𝑓 𝑥 ℎ
    積分区間が短いと思って𝑓(𝑥)を定数で近似する
    𝑥 𝑡 + ℎ = 𝑥 𝑡 + 𝑓 𝑥 𝑡 ℎ
    知りたい量 知ってる量
    未知の量𝑥(𝑡 + ℎ)を、既知の量𝑥 𝑡 と𝑓(𝑥(𝑡))だけで表現できた
    この数値積分法をオイラー法と呼ぶ

    View Slide

  61. 61
    85
    𝑑𝑝
    𝑑𝑡
    = −
    𝜕ℋ
    𝜕𝑞
    𝑑𝑞
    𝑑𝑡
    =
    𝜕ℋ
    𝜕𝑝
    運動方程式
    現在の値 𝑝 𝑡 , 𝑞(𝑡)
    未来の値
    𝑝 𝑡 + ℎ , 𝑞(𝑡 + ℎ)
    運動方程式を使って、未来の値を現在の値の関数として表現したい

    View Slide

  62. 62
    85
    𝑝 𝑡 + ℎ = 𝑝 𝑡 + ሶ
    𝑝 𝑡 ℎ = 𝑝 𝑡 −
    𝜕ℋ
    𝑑𝑞

    𝑞 𝑡 + ℎ = 𝑞 𝑡 + ሶ
    𝑞 𝑡 ℎ = 𝑞 𝑡 +
    𝜕ℋ
    𝑑𝑝

    運動方程式にオイラー法を適用する
    現在の値 現在の微分係数
    (傾き)
    t
    現在の値
    真の未来の値
    近似した未来の値
    すごくずれていきそうな気がする

    View Slide

  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
    調和振動子の運動方程式をオイラー法で解くコード
    𝑝 𝑡 + ℎ = 𝑝 𝑡 − 𝑞 𝑡 ℎ
    𝑞 𝑡 + ℎ = 𝑞 𝑡 + 𝑝 𝑡 ℎ

    View Slide

  64. 64
    85
    オイラー法で計算した調和振動子の軌道
    単位円を描くべき軌道が
    徐々に広がっている
    →エネルギーが増えている

    View Slide

  65. 65
    85
    𝑝(𝑡 + ℎ) = 𝑝(𝑡) − 𝑞(𝑡)ℎ
    𝑞 𝑡 + ℎ = 𝑞 𝑡 + 𝑝(𝑡)ℎ
    調和振動子にオイラー法を適用
    行列表示
    𝑝(𝑡 + ℎ)
    𝑞(𝑡 + ℎ)
    =
    1 −ℎ
    ℎ 1
    𝑝(𝑡)
    𝑞(𝑡)
    𝑝(𝑡 + ℎ)
    𝑞(𝑡 + ℎ)
    =
    cos ℎ −sin ℎ
    sin ℎ cos ℎ
    𝑝(𝑡)
    𝑞(𝑡)
    厳密解
    行列式= 1 + ℎ2 > 1
    行列式= 1
    位相空間の面積非保存がエネルギー非保存の原因

    View Slide

  66. 66
    85
    オイラー法は、テイラー展開の一次まで正しいコード
    →高次の積分法を構築する
    𝑞 𝑡 + ℎ = 𝑞 𝑡 + ሶ
    𝑞 𝑡 ℎ +

    𝑞 𝑡 ℎ2
    2

    𝑝 = 𝑓(𝑞), ሶ
    𝑞 = 𝑝
    以下の運動方程式を考える
    位置を二次までテイラー展開
    = 𝑞 𝑡 + 𝑝(𝑡)ℎ +
    𝑓(𝑡)ℎ2
    2
    ポテンシャル力

    View Slide

  67. 67
    85
    𝑝 𝑡 + ℎ = 𝑝 𝑡 + ሶ
    𝑝 𝑡 ℎ +

    𝑝 𝑡 ℎ2
    2
    運動量も二次までテイラー展開
    = 𝑝 𝑡 + 𝑓 𝑡 ℎ +

    𝑓 𝑡 ℎ2
    2
    𝑑𝑓
    𝑑𝑡

    𝑓 𝑡 + ℎ − 𝑓(𝑡)

    力の時間微分を差分近似して代入
    𝑝 𝑡 + ℎ = 𝑝 𝑡 +
    𝑓 𝑡 + ℎ + 𝑓(𝑡)
    2

    先に𝑞(𝑡 + ℎ)が求まっているので値がわかる

    View Slide

  68. 68
    85

    𝑝 = −𝑓(𝑞), ሶ
    𝑞 = 𝑝
    以下の運動方程式について
    𝑝 𝑡 + ℎ = 𝑝 𝑡 +
    𝑓 𝑡 + ℎ + 𝑓(𝑡)
    2

    𝑞(𝑡 + ℎ) = 𝑞 𝑡 + 𝑝(𝑡)ℎ +
    𝑓(𝑡)ℎ2
    2
    以下の二次の数値積分法を得た
    これをVelocity Verlet法と呼ぶ
    ※ 速度ベルレ法や、ベレの方法などとも呼ばれる

    View Slide

  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

    View Slide

  70. 70
    85
    Velocity Verlet法で計算した調和振動子の軌道
    軌道が完全な円を描いている

    View Slide

  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

    View Slide

  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
    エネルギーが真の値のまわりを揺らぎながらも保存する

    View Slide

  73. 73
    85
    𝑑
    𝑑𝑡
    = −𝑖𝐿 = −
    𝜕ℋ
    𝜕𝑞
    𝜕
    𝜕𝑝
    +
    𝜕ℋ
    𝜕𝑝
    𝜕
    𝜕𝑞
    リュービル演算子
    時間発展演算子
    𝑈(ℎ) = exp −𝑖ℎ𝐿
    一般に厳密に計算できない
    リュービル演算子から近似した時間発展演算子を作りたい

    View Slide

  74. 74
    85
    −𝑖𝐿 = −
    𝜕ℋ
    𝜕𝑞
    𝜕
    𝜕𝑝
    +
    𝜕ℋ
    𝜕𝑝
    𝜕
    𝜕𝑞
    = −𝑖𝐿𝑉
    + −𝑖𝐿𝐾
    −𝑖𝐿𝑉
    −𝑖𝐿𝐾
    リュービル演算子も二つに分解する
    ハミルトニアンが運動項とポテンシャル項に分離しているとする
    ℋ 𝑝, 𝑞 = 𝐾 𝑝 + 𝑉(𝑞)
    = −
    𝜕𝑉
    𝜕𝑞
    𝜕
    𝜕𝑝
    +
    𝜕𝐾
    𝜕𝑝
    𝜕
    𝜕𝑞

    View Slide

  75. 75
    85
    −𝑖𝐿𝑉
    を𝑝に二回演算すると消える
    −𝑖𝐿𝑉
    𝑝 = −
    𝜕ℋ
    𝜕𝑞
    𝜕
    𝜕𝑝
    𝑝 = −
    𝜕ℋ
    𝜕𝑞
    −𝑖𝐿𝑉
    2𝑝 = −
    𝜕ℋ
    𝜕𝑞
    𝜕
    𝜕𝑝

    𝜕ℋ
    𝜕𝑞
    = 0
    ℋ 𝑝, 𝑞 = 𝐾 𝑝 + 𝑉(𝑞)
    ハミルトニアンが運動項とポテンシャル項に分かれている時
    qだけの関数
    pによる偏微分
    −𝑖𝐿𝑉
    を𝑞に演算すると消える
    −𝑖𝐿𝑉
    𝑞 = −
    𝜕ℋ
    𝜕𝑞
    𝜕
    𝜕𝑝
    𝑞 = 0
    ※ − 𝑖𝐿𝐾
    も同様

    View Slide

  76. 76
    85
    𝑈𝐾
    ℎ = exp −𝑖ℎ𝐿𝐾
    = 1 + −𝑖ℎ𝐿𝐾
    +
    −𝑖ℎ𝐿𝐾
    2
    2
    + ⋯
    ここがゼロなので これ以降全部ゼロ
    = 1 + −𝑖ℎ𝐿𝐾
    𝑈𝑉
    ℎ = exp −𝑖ℎ𝐿𝑉
    = 1 + (−𝑖ℎ𝐿𝑉
    )
    同様に
    𝑈(ℎ) = exp −𝑖ℎ𝐿 は厳密に求められないが
    exp(−𝑖𝐿𝐾
    )やexp(−𝑖𝐿𝑉
    )は厳密に求めることができる

    View Slide

  77. 77
    85
    一般に演算子𝐴と𝐵について
    exp 𝐴 + 𝐵 ≠ exp 𝐴 exp 𝐵
    しかし、以下が成り立つ(Lie-Trotter公式)
    exp 𝐴 + 𝐵 = lim
    𝑛→∞
    exp 𝐴/𝑛 exp(𝐵/𝑛) 𝑛
    有限のnで止めることで、以下の近似式を得る
    exp ℎ(𝐴 + 𝐵) = exp ℎ𝐴 exp ℎ𝐵 + 𝑂 ℎ2
    exp ℎ(𝐴 + 𝐵) = exp ℎ𝐴/2 exp ℎ𝐵 exp ℎ𝐴/2 + 𝑂 ℎ3
    一次の公式
    二次の公式

    View Slide

  78. 78
    85
    −𝑖𝐿 = −𝑖𝐿𝑉
    + −𝑖𝐿𝐾
    リュービル演算子を二つに分解する
    求めたい時間発展演算子をLie-Trotter公式で分解する
    𝑈 ℎ = exp −𝑖ℎ𝐿
    = exp −𝑖ℎ𝐿𝑉
    exp −𝑖ℎ𝐿𝐾
    + 𝑂(ℎ2)
    = exp −𝑖ℎ𝐿𝑉
    /2 exp −𝑖ℎ𝐿𝐾
    exp −𝑖ℎ𝐿𝑉
    /2 + 𝑂(ℎ3)
    = exp −𝑖ℎ(𝐿𝑉
    + 𝐿𝐾
    )

    𝑈1
    (ℎ)

    𝑈2
    (ℎ)
    一次近似された時間発展演算子
    二次近似された時間発展演算子

    View Slide

  79. 79
    85
    厳密な時間発展
    𝑃
    𝑄
    = 𝑈 ℎ
    𝑝
    𝑞 = exp(−𝑖ℎ𝐿)
    𝑝
    𝑞
    一次近似された時間発展
    𝑃
    𝑄
    = ෩
    𝑈1

    𝑝
    𝑞 = exp(−𝑖ℎ𝐿𝐾
    )exp(−𝑖ℎ𝐿𝑉
    )
    𝑝
    𝑞
    最初にこれを演算して
    次にこれを演算する
    このように、指数分解公式で作られた時間発展演算子を使う
    数値積分法をシンプレクティック積分法と呼ぶ

    View Slide

  80. 80
    85
    𝑈𝑉
    (ℎ)を𝑝に演算してみる
    𝑈𝑉
    ℎ 𝑝 = exp −𝑖ℎ𝐿𝑉
    𝑝
    = (1 − 𝑖ℎ𝐿𝑉
    )𝑝
    = 1 − ℎ
    𝜕ℋ
    𝜕𝑞
    𝜕
    𝜕𝑝
    𝑝
    = 𝑝 − ℎ
    𝜕𝑉
    𝜕𝑞
    それぞれ、お互いを無視して時間ℎだけ運動
    𝑈𝐾
    (ℎ)を𝑞に演算してみる
    𝑈𝐾
    ℎ 𝑞 = exp −𝑖ℎ𝐿𝐾
    𝑞
    = 1 − 𝑖ℎ𝐿𝐾
    𝑞
    = 1 + ℎ
    𝜕ℋ
    𝜕𝑝
    𝜕
    𝜕𝑞
    𝑞
    = 𝑞 + ℎ
    𝜕𝐾
    𝜕𝑞
    𝑝の時間をℎだけ進める 𝑞の時間をℎだけ進める

    View Slide

  81. 81
    85
    調和振動子の場合
    𝑈𝑉

    𝑝
    𝑞 =
    𝑝 − ℎ𝑞
    𝑞
    =
    1 −ℎ
    0 1
    𝑝
    𝑞
    𝑈𝐾

    𝑝
    𝑞 =
    𝑝
    𝑞 + ℎ𝑝 =
    1 0
    ℎ 1
    𝑝
    𝑞
    𝑝の時間をℎだけ進める演算子
    𝑞の時間をℎだけ進める演算子
    二つともかけると
    𝑈𝐾
    ℎ 𝑈𝑉

    𝑝
    𝑞 =
    1 −ℎ
    ℎ 1 − ℎ2
    𝑝
    𝑞
    行列式が1

    View Slide

  82. 82
    85
    𝑝, 𝑞 から(𝑃, 𝑄)への変換
    𝑃
    𝑄
    = 𝑈
    𝑝
    𝑞
    において
    𝑑𝑃𝑑𝑄 = 𝑑𝑝𝑑𝑞
    が満たされる時、この変換をシンプレクティック変換と呼ぶ
    𝑑𝑃𝑑𝑄 =
    𝜕𝑝
    𝑃 𝜕𝑞
    𝑃
    𝜕𝑝
    𝑄 𝜕𝑞
    𝑄
    𝑑𝑝𝑑𝑞 であるから
    シンプレクティック変換とは、変換のヤコビアンが1であること
    ※ 時間発展演算子を行列表示したときの行列式が1

    View Slide

  83. 83
    85

    𝑈2
    (ℎ) = exp −𝑖ℎ𝐿𝑉
    /2 exp −𝑖ℎ𝐿𝐾
    exp −𝑖ℎ𝐿𝑉
    /2
    二次の指数分解公式
    (1)
    (2)
    (3)
    (1)
    𝑝 𝑡 + ℎ/2 = 𝑝(𝑡) +
    𝑓 𝑡
    2
    (2)
    𝑄 = 𝑞(𝑡) + 𝑝 𝑡 + ℎ/2 ℎ
    (3)
    𝑃 = 𝑝(𝑡 + ℎ/2) + 𝑓 𝑡 + ℎ ℎ
    現在の力がそのまま続いたとしてℎ/2だけ時間を進める
    (1)で得た速度がそのまま続いたとしてℎだけ時間をすすめる
    (2)の力がそのまま続いたとしてℎ/2だけ時間を進める

    View Slide

  84. 84
    85
    𝑝 𝑡 + ℎ/2 = 𝑝(𝑡) +
    𝑓 𝑡
    2

    𝑄 = 𝑞(𝑡) + 𝑝 𝑡 + ℎ/2 ℎ
    𝑃 = 𝑝(𝑡 + ℎ/2) + 𝑓 𝑡 + ℎ ℎ
    𝑝(𝑡 + ℎ/2)を消去して整理するとVelocity Verlet公式を得る
    𝑃 = 𝑝 𝑡 +
    𝑓 𝑡 + ℎ + 𝑓(𝑡)
    2

    𝑄 = 𝑞 𝑡 + 𝑝(𝑡)ℎ +
    𝑓(𝑡)ℎ2
    2
    Velocity Verlet法は、シンプレクティック積分

    View Slide

  85. 85
    85
    • 数値積分とは、現在時刻𝑡の点(𝑝, 𝑞)から、時刻𝑡 + ℎの
    点(𝑃, 𝑄)を得る写像である
    • 写像を表す変換のヤコビアンが1であるとき、位相空
    間の体積が保存される。このような写像をシンプレク
    ティック写像と呼ぶ
    • 変換がシンプレクティック写像であるような数値積分
    法をシンプレクティック積分と呼ぶ
    • シンプレクティック積分は、Lie-Trotterの指数分解公
    式から構築できる
    (𝑝, 𝑞)
    (𝑃, 𝑄)
    𝑑𝑃𝑑𝑄 = 𝑑𝑝𝑑𝑞

    View Slide