ハローラスク

時間から自分を無理に引き離すこのような安定性あるいは確実性に到達するためには、時間が必要である。時間なしですませるためには、時間が必要である。時間のなかで時間を従わせ、みずから時間に従わなければならない。

シンプレクティック法をシミュレーションで理解する

なんでこんな記事書いたの?

  • 解析力学をやっていたときシンプレクティック数値積分という考え方を知り興味を持った
  • 簡単な力学系をシミュレーションしたいみたいな気持ちが前からあった
  • Runge-Kutta法のような数値解析をpythonで実装させたいみたいな気持ちが前からあった
  • とにかくAセメスターの試験勉強をやりたくない

やること

正準運動方程式
f:id:HelloRusk:20171221215906p:plain:w120
を様々な方法で数値解析してみます.

調和振動子

f:id:HelloRusk:20171222085540p:plain:w400
一番簡単な例から見てみます. 原点の位置で自然長となるようにばね定数  k のバネを質量  m のおもりにつなげます. 変位を  q とします.  q=1 でそっとおもりから手を離すのを初期条件としましょう. Lagrangianは
f:id:HelloRusk:20171222090603p:plain:w150
また  q の共役運動量を  p とすると  p = \displaystyle\frac{\partial L}{\partial \dot{q}} = m\dot{q} より  \dot{q} = \displaystyle\frac{p}{m} なのでHamiltonianは
f:id:HelloRusk:20171222091825p:plain:w150
ここで, このHamiltonian  H は時間  t に陽に依っていないので時間保存します(これから出てくる例では全部時間保存します). これはいわゆるエネルギー保存則です. エネルギー保存則がこの記事で重要なポイントです.
簡単のため, 今回は  m=k=1 とおいて無次元化します. そうすると
f:id:HelloRusk:20171222092648p:plain:w150
となり, 正準運動方程式
f:id:HelloRusk:20171222092949p:plain:w100
となります. 言うまでもなくこの微分方程式 q=\sin t, p=\cos t という解析解があるのですが, ここでは数値解を求めてみましょう.
常微分方程式の数値解法としてEuler法というものが古典的に知られています. これは
f:id:HelloRusk:20171222095107p:plain:w120

f:id:HelloRusk:20171222095409p:plain:w250
と置き換えるものです(dtは時刻の刻み幅). Taylor展開の2次以降を無視しているのと同じなので1次の精度だと言われたりします. 今回の正準運動方程式にあてはめてみると
f:id:HelloRusk:20171222100822p:plain:w200
では実際にこの計算をやってみて, エネルギーHが保存されているかどうか確認してみましょう. 時刻の刻み幅を0.1, 回数を100回にしてみます(すなわち10秒間やるということです).

import matplotlib.pyplot as plt
import numpy as np

dt = 0.1
times = 100
T = np.array([i * dt for i in range(times + 1)])
q_0 = 1
p_0 = 0


def f(p):
    return p


def g(q):
    return -q


def energy(p, q):
    return (p**2 + q**2)/2


def euler(x, y, h):
    x_next = x + f(y) * h
    y_next = y + g(x) * h
    return x_next, y_next


def answer(func):
    q = np.array([q_0])
    p = np.array([p_0])
    for i in range(times):
        q_temp, p_temp = func(q[i], p[i], dt)
        q = np.append(q, q_temp)
        p = np.append(p, p_temp)
    return q, p


q1, p1 = answer(euler)

plt.figure(figsize=(8, 8))
plt.plot(T, energy(q1, p1), label="Euler", color="Red")
plt.legend()
plt.xlabel("Time")
plt.ylabel("Energy")
plt.show()

出力はこうなります.
f:id:HelloRusk:20171222120012p:plain:w600
全然エネルギーが保存していないことがわかります.  q p を座標軸にとった位相平面を見てみると何が起こっているのかよくわかります. 本来の解析解ならば  H=\displaystyle\frac{1}{2}(p^{2} + q^{2}) が保存することから円を描くはずなのですが
f:id:HelloRusk:20171222121044p:plain:w600
と渦を巻いてしまっています. ここで気分で(?)
f:id:HelloRusk:20171222121858p:plain:w240
とさっきの式を書き換えてみます. 違うのは赤い部分だけです. これは半陰的 (semi-implicit) Euler法と呼ばれるそうです. (一般的に漸化式の右辺がn+1の項を含まない式なら陽的、含む式なら陰的と呼ぶらしいです). これを使って同じことをやってみます.

def semi_implicit_euler(x, y, h):
    x_next = x + f(y) * h
    y_next = y + g(x_next) * h
    return x_next, y_next

を定義してグラフを書くと
f:id:HelloRusk:20171222123456p:plain:w600
エネルギーが多少の振動はあるものの大きなズレはありません. 位相平面の場合でも
f:id:HelloRusk:20171222123724p:plain:w600
と綺麗な円にはなっていないものの点が元の位置に戻って周期運動をしています. 回数を100回から1000回に増やすともうその差は歴然とします.
f:id:HelloRusk:20171222123945p:plain:w600
半陰的Euler法すっごーい!という感じですが, これは半陰的Euler法がシンプレクティック性を持っているからです. これについてもう少し述べてみたいと思います.

シンプレクティック性とは

あるHamiltonian  H(q, p, t) と正準変数  (q, p) を用いた正準運動方程式
f:id:HelloRusk:20171221215906p:plain:w120
によって記述される系が, 別のHamiltonian  K(Q, P, t) と正準変数  (Q, P) を用いた正準運動方程式
f:id:HelloRusk:20171223200124p:plain:w120
によって同等に記述されるとき,  (q, p) (Q, P) に移すような変換を正準変換と呼びます. 正準変換の必要十分条件の表現形式はシンプレクティック行列, ポアソン括弧, 正準変換母関数など色々バリエーションがありますが, ここではポアソン括弧を紹介してみます. ポアソン括弧とは  (q, p) の関数  f, g に対して
f:id:HelloRusk:20171223205945p:plain:w350
と定義されるものです. このポアソン括弧を用いると, ある変数変換  (q, p) \rightarrow (Q, P) が正準変換である必要十分条件
f:id:HelloRusk:20171223211801p:plain:w400
と表せます. 特に一変数の場合だと
f:id:HelloRusk:20171223212023p:plain:w120
となります. これを使って, 先程の調和振動子の例において  (q_n, p_n) \rightarrow (q_{n+1}, p_{n+1}) が正準変換かどうかを確認してみましょう. まず普通のEuler法
f:id:HelloRusk:20171222100822p:plain:w200
の場合
f:id:HelloRusk:20171223214215p:plain:w400
となり, 正準変換ではないことが分かります. 一方, 半陰的Euler法
f:id:HelloRusk:20171222121858p:plain:w240
の場合
f:id:HelloRusk:20171223215230p:plain:w500
と綺麗に1になって, 正準変換であることが分かります. このように, 正準変換を保った数値解析の方法をシンプレクティック法(シンプレクティック数値積分法)と呼びます.
上記の例はHamiltonianが
f:id:HelloRusk:20171222092648p:plain:w150
のときを考えていますが, 一般にHamiltonianが  q p で完全に分離した形で書ける場合, 半陰的Euler法はシンプレクティック性を持っています.

ここまでシンプレクティック性についてグダグダと述べてきましたが, では一体シンプレクティック性の何が良いんでしょうか? それはシンプレクティック法ではHamiltonian(エネルギー)の計算値と厳密な値の誤差が蓄積することなく, あるオーダーの中に留まるという点にあります(この証明は私の数学力がなくてあまり理解できませんでした. 興味のある人は「吉田春夫 シンプレクティック」でググってみてください).
逆に言うと, シンプレクティック性を持っていない数値解析は, 本来ならばエネルギーが一定のはずの系でエネルギーがどんどん変わってしまう, みたいなことが起こってマズいわけです. 『解析力学量子論』(須藤 靖)という本の61〜62ページにこんなことが書いてあります.

たとえば「我々の太陽系の運動が安定かどうか」という本質的な難問がある. 太陽系が誕生してから46億年とされているので, その時間スケールでの安定性は経験的に証明されている. しかしこれから数十億年後, 数百億年後はどうなのであろう. 重力多体問題は3体以上になると解析的には解けないから数値積分するしかないが, その間地球は数十億回, 数百億回も公転する. わずかな誤差でもただちに蓄積し真の解からずれてしまうであろう. とすれば, 太陽系が力学系としてもつ真の不安定性なのか, あるいは数値積分の精度に起因する数値的な不安定性なのかを判別することはきわめて難しい.
高エネルギー実験においてもっとも重要な研究手段である加速器において, その中を運動する荷電粒子の軌道計算も同様である. 我が国の高エネルギー加速器研究機構の電子・陽電子衝突型加速器では, 電子や陽電子は光速の0.999999998倍という速度で1周3kmのリングを10億回まわる. しかもわずか2ミクロンというサイズの標的に衝突するように, 軌道を正確に設計する必要がある. そのための数値計算にはきわめて高い精度が要求されることが理解できよう.
これらの問題においてはシンプレクティック積分法を用いることなしには, 必要な精度の数値計算は行えない.

バネ振り子

御託を並べるのはこれくらいにして, 次の例を見ていきましょう.
f:id:HelloRusk:20171224154452p:plain:w300
単振り子の糸がバネになったバージョンです. 図の通り  q_1,  q_2 を定めます. おもりの質量は  m, バネ定数は  k, バネの自然長は  a とします. 図の通り一旦  (x, y) 座標を設定して考えるとLagrangianは
f:id:HelloRusk:20171224163201p:plain:w320
で, これに
f:id:HelloRusk:20171224161650p:plain:w250
を代入して,
f:id:HelloRusk:20171224163307p:plain:w400
また  q_1, q_2 の共役運動量をそれぞれ  p_1, p_2 とすると,
f:id:HelloRusk:20171225104114p:plain:w250
なのでHamiltonianは
f:id:HelloRusk:20171225105412p:plain:w400
で, 正準運動方程式
f:id:HelloRusk:20171225110848p:plain:w300
となります. なかなか見づらいので例によって無次元化します.  m = g = k = a = 1 とすると,
f:id:HelloRusk:20171225111733p:plain:w360
f:id:HelloRusk:20171225111826p:plain:w270
この方程式は思いっきり非線形で, 見るからに解析解がなさそうです. というわけで今回はEuler法だけでなくRunge-Kutta法も試してみることにします. Runge-Kutta法には様々なバリエーションがありますが, ここで使うのは最も有名な4段4次陽的Runge-Kutta法(RK4)です. これは
f:id:HelloRusk:20171222095107p:plain:w120

f:id:HelloRusk:20171225114631p:plain:w350
と置き換えるもので, 係数が綺麗な形をしている割にTaylor展開の4次項まで再現できるので(=精度が4次なので)重宝されます. ただし, シンプレクティック性は持っていません. ということはシミュレーションしたらエネルギーが発散してしまうのが確認できるはずです.
では, 前章で紹介した半陰的Euler法は今回はどうなるでしょうか. 今回は二変数なので, シンプレクティック性を持つ条件は
f:id:HelloRusk:20171223211801p:plain:w400
となりますが, これを真面目に考えるのは面倒です. そこで, とりあえず q を求める時はnの項を使って,  p を求める時はn+1の項を使えば良いんじゃね?」と雑に考えました. つまり, もとの正準運動方程式
f:id:HelloRusk:20171225131531p:plain:w180
という形をしていますが, これを
f:id:HelloRusk:20171225132126p:plain:w300
と置き換えるというアイディアです(見やすさのために添え字の位置を変えました). これは本来の半陰的Euler法とは別物のような気がするので, 半陰的Euler法もどきと呼びましょう. 半陰的Euler法もどきにおいては, 4式目を計算してからでないと3式目が計算できないので, プログラムを書く時は順番を変えることにします.
以上のEuler法, 半陰的Euler法もどき, Runge-Kutta法の3つの方法で, エネルギーが保存されているかどうかを確かめてみましょう. 初期条件は「自然長で45°角度をつけた状態でそっと手を離す」, すなわち  q_1 =1, q_2 = \displaystyle\frac{\pi}{4}, p_1 = p_2 = 0 とします. また, 時刻の刻み幅を0.01とします. とりあえず50回(0.5秒間)やってみます.

import matplotlib.pyplot as plt
import numpy as np

dt = 0.01
times = 50
T = np.array([i * dt for i in range(times + 1)])
q1_0 = 1
q2_0 = np.pi / 4
p1_0 = 0
p2_0 = 0


def f1(p1):
    return p1


def f2(q1, p2):
    return p2 / (q1**2)


def g1(q1, q2, p2):
    return (p2**2) / (q1**3) + np.cos(q2) - (q1 - 1)


def g2(q1, q2):
    return - q1 * np.sin(q2)


def energy(q1, q2, p1, p2):
    return (p1**2) / 2 + (p2**2) / (2 * (q1**2)) - q1 * np.cos(q2) + ((q1 - 1)**2) / 2


def euler(x1, x2, y1, y2, h):
    x1_next = x1 + f1(y1) * h
    x2_next = x2 + f2(x1, y2) * h
    y1_next = y1 + g1(x1, x2, y2) * h
    y2_next = y2 + g2(x1, x2) * h
    return x1_next, x2_next, y1_next, y2_next


def semi_implicit_euler_modoki(x1, x2, y1, y2, h):
    x1_next = x1 + f1(y1) * h
    x2_next = x2 + f2(x1, y2) * h
    y2_next = y2 + g2(x1_next, x2_next) * h
    y1_next = y1 + g1(x1_next, x2_next, y2_next) * h
    return x1_next, x2_next, y1_next, y2_next


def runge_kutta(x1, x2, y1, y2, h):
    k1 = f1(y1)
    k2 = f1(y1 + k1 * (h / 2))
    k3 = f1(y1 + k2 * (h / 2))
    k4 = f1(y1 + k3 * h)
    x1_next = x1 + (k1 + 2 * k2 + 2 * k3 + k4) * (h / 6)
    k1 = f2(x1, y2)
    k2 = f2(x1 + k1 * (h / 2), y2 + k1 * (h / 2))
    k3 = f2(x1 + k2 * (h / 2), y2 + k2 * (h / 2))
    k4 = f2(x1 + k3 * h, y2 + k3 * h)
    x2_next = x2 + (k1 + 2 * k2 + 2 * k3 + k4) * (h / 6)
    k1 = g1(x1, x2, y2)
    k2 = g1(x1 + k1 * (h / 2), x2 + k1 * (h / 2), y2 + k1 * (h / 2))
    k3 = g1(x1 + k2 * (h / 2), x2 + k2 * (h / 2), y2 + k2 * (h / 2))
    k4 = g1(x1 + k3 * h, x2 + k3 * h, y2 + k3 * h)
    y1_next = y1 + (k1 + 2 * k2 + 2 * k3 + k4) * (h / 6)
    k1 = g2(x1, x2)
    k2 = g2(x1 + k1 * (h / 2), x2 + k1 * (h / 2))
    k3 = g2(x1 + k2 * (h / 2), x2 + k2 * (h / 2))
    k4 = g2(x1 + k3 * h, x2 + k3 * h)
    y2_next = y2 + (k1 + 2 * k2 + 2 * k3 + k4) * (h / 6)
    return x1_next, x2_next, y1_next, y2_next


def answer(func):
    q1 = np.array([q1_0])
    q2 = np.array([q2_0])
    p1 = np.array([p1_0])
    p2 = np.array([p2_0])
    for i in range(times):
        q1_temp, q2_temp, p1_temp, p2_temp = func(q1[i], q2[i], p1[i], p2[i], dt)
        q1 = np.append(q1, q1_temp)
        q2 = np.append(q2, q2_temp)
        p1 = np.append(p1, p1_temp)
        p2 = np.append(p2, p2_temp)
    return q1, q2, p1, p2


q1_e, q2_e, p1_e, p2_e = answer(euler)
q1_s, q2_s, p1_s, p2_s = answer(semi_implicit_euler_modoki)
q1_r, q2_r, p1_r, p2_r = answer(runge_kutta)

plt.figure(figsize=(8, 8))
plt.plot(T, energy(q1_e, q2_e, p1_e, p2_e), label="Euler", color="Red")
plt.plot(T, energy(q1_s, q2_s, p1_s, p2_s), label="SemiImplicitEulerModoki", color="Blue")
plt.plot(T, energy(q1_r, q2_r, p1_r, p2_r), label="RungeKutta", color="Orange")
plt.legend()
plt.xlabel("Time")
plt.ylabel("Energy")
plt.show()

出力はこうなります.
f:id:HelloRusk:20171225142314p:plain:w600
今の所はRunge-Kutta法が一番エネルギーのズレが少なく優勢です. 流石4次の精度を誇るだけあります. では100回に増やしてみましょう.
f:id:HelloRusk:20171225142722p:plain:w600
Runge-Kutta法の優勢に陰りが見え始めてきました. 500回に増やしてみましょう.
f:id:HelloRusk:20171225143036p:plain:w600
ここまで来ると, Euler法, Runge-Kutta法のエネルギーのズレが明白になってしまいました. 一方半陰的Euler法もどきは快調な動きを見せています. 5000回に増やしてみましょう.
f:id:HelloRusk:20171225143431p:plain:w600
一目瞭然です. Euler法, Runge-Kutta法がシンプレクティック性を持っていないこと, そして半陰的Euler法もどきは(おそらく)シンプレクティック性を持っていることがよく分かります. ちなみに半陰的Euler法もどきだけで20000回でグラフを作るとこんな感じです.
f:id:HelloRusk:20171225185753p:plain:w600
エネルギーの誤差は  \pm 0.005 の範囲に収まっています. かなり正確ですね.
(実を言うと半陰的Euler法よりもさらに精度の高いシンプレクティック法として, Verlet法やRuthの公式, Sanz-Sernaの公式など色々存在するらしいのですが, 今回の半陰的Euler法もどきを発展させてそれらの方法を適用しようとしたところ, エネルギーが発散してしまい失敗しました. おそらく「もどき」の方法を無理やり適用しようとしてシンプレクティック性の条件を満たさなくなってしまったようです. これらの公式についてちゃんと知りたいという人は「ハミルトン系に対する離散変数法」でググってみてください.)
さて, ここまでエネルギーの話ばかりしてきましたが, じゃあ実際どんな運動をしているんだ? というのが気になります. というわけでmatplotlibでアニメーションを作って確認してみましょう. 時刻の刻み幅が0.01だと動画を作るのが難しそうだったので今回は0.1にしてあります. また動画があまりに長いとファイルが大きすぎてブログに上げられなくなってしまう恐れがあるので, 時間は30秒間にしました.

import matplotlib.pyplot as plt
import matplotlib.animation as anime
import numpy as np

dt = 0.1
times = 300
T = np.array([i * dt for i in range(times + 1)])
q1_0 = 1
q2_0 = np.pi / 4
p1_0 = 0
p2_0 = 0


def f1(p1):
    return p1


def f2(q1, p2):
    return p2 / (q1**2)


def g1(q1, q2, p2):
    return (p2**2) / (q1**3) + np.cos(q2) - (q1 - 1)


def g2(q1, q2):
    return - q1 * np.sin(q2)


def energy(q1, q2, p1, p2):
    return (p1**2) / 2 + (p2**2) / (2 * (q1**2)) - q1 * np.cos(q2) + ((q1 - 1)**2) / 2


def semi_implicit_euler_modoki(x1, x2, y1, y2, h):
    x1_next = x1 + f1(y1) * h
    x2_next = x2 + f2(x1, y2) * h
    y2_next = y2 + g2(x1_next, x2_next) * h
    y1_next = y1 + g1(x1_next, x2_next, y2_next) * h
    return x1_next, x2_next, y1_next, y2_next


def answer(func):
    q1 = np.array([q1_0])
    q2 = np.array([q2_0])
    p1 = np.array([p1_0])
    p2 = np.array([p2_0])
    for i in range(times):
        q1_temp, q2_temp, p1_temp, p2_temp = func(q1[i], q2[i], p1[i], p2[i], dt)
        q1 = np.append(q1, q1_temp)
        q2 = np.append(q2, q2_temp)
        p1 = np.append(p1, p1_temp)
        p2 = np.append(p2, p2_temp)
    return q1, q2, p1, p2


q1, q2, p1, p2 = answer(semi_implicit_euler_modoki)

fig = plt.figure(figsize=(8, 8))

def update(i):
    if i != 0:
        plt.cla()
    plt.xlim([-1.5, 1.5])
    plt.ylim([-3, 0])
    plt.title("T = " + str(round(i * dt, 1)) + " , Nobi: " + str(round(q1[i] - 1, 1)))
    plt.plot(q1[i] * np.sin(q2[i]), - q1[i] * np.cos(q2[i]), "o", color="blue")
    plt.plot([0, q1[i] * np.sin(q2[i])], [0, - q1[i] * np.cos(q2[i])], "--", color="blue")

ani = anime.FuncAnimation(fig, update, interval=1, frames=times + 1)
# ani.save("mov.gif", writer="imagemagick") で動画を保存
plt.show()

出力はこんな感じです.

f:id:HelloRusk:20171226130629g:plain:w600
k = 1, q1(0) = 1, q2(0) = π/4, p1(0) = p2(0) = 0

グラフの上に時刻とその時のバネの伸びも表示してみました(少し技術的な話をすると, matplotlib.animationにはArtistAnimationとFuncAnimationの2種類があるのですが, FuncAnimationの場合1フレームごとに関数を実行してくれるので, このようにその時々の時刻や伸びを表示させたいときに都合が良いです).
さて, 出力結果を見ると予想以上に重力の影響が大きい気がします. そこで,  k=5 に変えてみます.
f:id:HelloRusk:20171226202224g:plain:w600
k = 5, q1(0) = 1, q2(0) = π/4, p1(0) = p2(0) = 0

大分伸びにくくなった気がします. このバネ定数のまま, 今度は伸ばした状態から始めたり縮めた状態から始めたりしてみます.
f:id:HelloRusk:20171226203215g:plain:w600
k = 5, q1(0) = 1.5, q2(0) = π/4, p1(0) = p2(0) = 0

f:id:HelloRusk:20171226203640g:plain:w600
k = 5, q1(0) = 0.5, q2(0) = π/4, p1(0) = p2(0) = 0

動きが面白くなってきました.  p の初期値をいじると初速をつけることもできます. 例えばこういうものが作れます.

f:id:HelloRusk:20171227123031g:plain:w600
k = 5, q1(0) = 0.5, q2(0) = π, p1(0) = 1, p2(0) = 2

目の保養になりますね.

本当はこの後二重振り子とか二重バネ振り子とかもやってみようと思っていたのですが, 思っていたよりも面倒そうだったのでやめました. 気が向いたらいつかやるかもしれません.

シンプレクティック法の弱点

ここまで「エネルギーの誤差が蓄積しない」というシンプレクティック法の強みについて述べてきたわけですが, シンプレクティック法は決して万能というわけではありません. 実は, 「エネルギーの誤差が蓄積しない」という特性は, 時刻の刻み幅を局所で変えていくと成り立たなくなってしまいます.
これは例えば星の軌道を精密に計算しようとするときに問題で, というのも実際に軌道をシミュレーションする際, 計算精度を落とさないために可変時間ステップ(時刻の刻み幅を可変にすること)で計算するのが普通だからだそうです.
いかに天体運動を精密に計算するか, という課題に関しては現在も様々な研究が続けられているようです(例えば「高度なN体シミュレーション法」でググってみると色々な数値解析の手法が出てきます). こういう研究においてはもはや物理学と情報科学が不可分の関係にあるということが素人目でも分かります.