動点Pはわからなくても質点Pくらいはシミュレーションしたい（シンプレクティック数値積分法）

はじめに

数値解析学を学ぶと微分方程式の解法として一番簡単なオイラー法が紹介されると思います。そして実際に実装してみると確かに解くことはできるけど、誤差がどんどん大きくなるところに目を瞑っていられなくなります。そこで調和振動子などの物理系に対し、安定して長期間のシミュレーションが行えるシンプレクティック数値積分法を紹介したいと思います。

必要とする前提知識
・高校レベルの古典力学の知識。
・全微分あたりまでの微積分の知識と行列の計算
その他
・$f$の定義を$f:=f(x)$と書く。

問題提起

よく古典力学の教科書とかで見る次の状況を考えましょう。

上図のバネをびよーんって伸ばしてから離して振動する質点Pの$t$秒後の位置$x=x(t)$と運動量$p=p(t)$を知りたいです。これは調和振動子と呼ばれる物理系^[1]になります。このとき質点の進行方向に対して逆向きにバネの力が働くので運動方程式は次のようになります。

$$m\frac{d^{2}x}{dt^{2}}=-kx$$

ここで$m,k$は質点の質量とバネ定数になりますが簡単のため両方とも$m=1,k=1$としてしまいましょう。そして運動量$p$は質量$m$と速度$\frac{dx}{dt}$の積であり、

$$\frac{dx}{dt}=p$$

と書けます。そうすると運動方程式は

$$\frac{dp}{dt}=-x$$

となります。この$\frac{dx}{dt}$と$\frac{dp}{dt}$に関する上二つの微分方程式から次の解析解が得られます。

$$x(t)=x_{0}\cos t,\;p(t)=-x_{0}\sin t$$

はじめ質点Pは$x(0)=x_{0}$に位置し、完全に止まっていること$p(0)=0$を想定しました。この解から$x(t)$も$p(t)$も振動することが分かりました。

次に微分方程式の解き方は知らないけど微分は知っているAさんが現れました。Aさんは非常に近い二点間の変化率が微分だと理解していたので、先ほどの式を$\epsilon\approx0$を使って近似してみました。

$$\frac{dx}{dt}\approx\frac{x_{n+1}-x_{n}}{\epsilon}=p_{n},\;\frac{dp}{dt}\approx\frac{p_{n+1}-p_{n}}{\epsilon}=-x_{n}$$

この式から次の漸化式が導けます。

$$\begin{align*} x_{n+1}&=x_{n}+\epsilon p_{n}\\ p_{n+1}&=-\epsilon x_{n}+p_{n} \end{align*}$$

$t=0$に対応する$n$を$n=0$すれば、上式で再帰的に$n$ステップ目の位置と運動量$x_{n},p_{n}$が求まります。微分の原理に基づき時間を$\epsilon$で離散化して再帰的に解く手法をオイラー法と言います。せっかくなのでオイラー法から導かれる位置と運動量の時間変化を見てましょう。ついでに解析解とも比較しておきます。まずは位置$x$です。

位置$x$の時間変化：ピンク実線: 解析解、青破線: オイラー法

横軸が時間、縦軸が位置です。また$\epsilon=0.05,x_{0}=1$で設定しました。時間が経つにつれ解析解から遠のいていくのが分かります。次に運動量です。

運動量$p$の時間変化：ピンク実線: 解析解、青破線: オイラー法

運動量も同様に解析解から遠のいていきます。

この誤差をどう解決するか？

先ほどオイラー法による誤差として、解析解から遠のいていくところを確認しました。そしてより時間が経てばこの誤差は大きなものになっていきます。そこでエネルギー保存の観点からどうすれば解決できるか考えたいと思います。

エネルギーとハミルトニアン

調和振動子におけるエネルギーは、質点の運動エネルギー$K$とバネの弾性エネルギー$V$です。それぞれ次のように書けます。

$$K=\frac{1}{2}mv^{2}=\frac{p^{2}}{2m},\;V=\frac{1}{2}kx^{2}$$

$v=\frac{dx}{dt}$は質点の速度です。この系において摩擦などの外的要因がないとすれば、全エネルギー

$$E=K+V$$

は保存します。これを少し一般化しましょう。$K$はよく見ると運動量$p$の関数です。また$V$は$x$の関数です。そうするとその和は$x$と$p$の関数になります。なので新たに

$$H=H(x,p):=K+V=\frac{p^{2}}{2m}+\frac{1}{2}kx^{2}$$

を定義しましょう。この$H$はハミルトニアンと言います。ここで$m=1,k=1$の仮定を思い出すとこの調和振動子のハミルトニアンは

$$H=\frac{1}{2}\left(p^{2}+x^{2}\right)$$

と書けます。このハミルトニアンに$x_{n},p_{n}$と$x_{n+1},p_{n+1}$を代入し$n$と$n+1$ステップの関係を調べると、

$$\frac{1}{2}\left(p_{n+1}^{2}+x_{n+1}^{2}\right)=\frac{1}{2}\left(p_{n}^{2}+x_{n}^{2}\right)(1+\epsilon^{2})$$

ステップを一つ進めるごとにエネルギーが$1+\epsilon^{2}$倍に蓄積することが確認できます。つまりオイラー法では本来満たしてほしいエネルギーの保存が原理的に満たされないわけです。

オイラー法を工夫する

これらの議論からステップを進めてもエネルギーが蓄積しない手法を考えましょう。結論から言うと次のように改良されます。

$$\begin{align*} x_{n+1}&=x_{n}+\epsilon p_{n}\\ p_{n+1}&=-\epsilon x_{n+1}+p_{n}\left(=-\epsilon x_{n} + (1-\epsilon^{2}) p_{n}\right) \end{align*}$$

オイラー法と何が違うかというと$p_{n+1}$への更新の際、$x_{n}$ではなく$x_{n+1}$を用いている点です。このとき確かにエネルギーは保存されるのですが、実際に保存するハミルトニアンは次のような形になります。

$$H=\frac{1}{2}\left(p^{2}+x^{2}\right)+\frac{\epsilon}{2}px$$

$\frac{\epsilon}{2}px$という余計な項が含まれていますが、試しに$n$と$n+1$ステップ目の関係を調べると

$$\frac{1}{2}\left(p_{n+1}^{2}+x_{n+1}^{2}\right)+\frac{\epsilon}{2}p_{n+1}x_{n+1}= \frac{1}{2}\left(p_{n}^{2}+x_{n}^{2}\right)+\frac{\epsilon}{2}p_{n}x_{n}$$

と見事にエネルギーが保存されます^[2]。とりあえずこの改良を通じてオイラー法と比較しつつシミュレーションした結果を示します。まずは位置$x$です。

位置$x$の時間変化：ピンク実線: 解析解、青破線: オイラー法、緑破線：シンプレクティック数値積分法
次に運動量$p$です。

運動量$p$の時間変化：ピンク実線: 解析解、青破線: オイラー法、緑破線：シンプレクティック数値積分法

正しく解析解に沿ってシミュレーションできていることが確認できました！この数値解析法をシンプレクティック数値積分法と言います。

ここで注意しなければいけないのは「じゃあシンプレクティック数値積分法はエネルギーを保存するシミュレーション方法なのね！」と理解してしまうことです。あくまでも調和振動子本来のハミルトニアンは$\frac{1}{2}(p^{2}+x^{2})$なので、これで時間ごとのエネルギーを評価すると厳密に保存せず一定の周期で振動していることが確認できます。

時間ごとのエネルギーの大きさ
ピンク実線: 解析解、青破線: オイラー法、緑破線：シンプレクティック数値積分法

上の図を見るとステップを進めるとオイラー法はエネルギーが大きくなり続けるのに対して、シンプレクティック数値積分法では一定の範囲で振動し続けます。これは$\frac{\epsilon}{2}px$を含むハミルトニアンが保存するからです。いずれにしてもエネルギーが大きくなり続けないという点で、うまくできた数値解析法だといえます。

シンプレクティック数値積分法の原理

ではシンプレクティック数値積分法どんな原理で成り立っているのでしょうか？ここでは解析力学という高校や大学１年生で学ぶ古典力学を数学的に体系化したものを使って説明します。

ハミルトンの正準方程式と正準変換

解析力学で学ぶ物体の運動を表す式には、ラグランジュの運動方程式とハミルトンの正準方程式があります。どちらもとても重要なのですがハミルトンの正準方程式を紹介します。ハミルトンの正準方程式は次のように定義されます。

$$\begin{align*} \frac{dx}{dt}&=\frac{\partial H}{\partial p}\\ \frac{dp}{dt}&=-\frac{\partial H}{\partial x} \end{align*}$$

$H$は例のごとくハミルトニアンです。$x,p$も位置と運動量になります^[3]。このハミルトニアンに調和振動子のものを代入し計算してみてください。記事のはじめに導いた$\frac{dx}{dt}=p,\frac{dp}{dt}=-x$が導出されます。さてここで$x,p$から導かれる新たな位置と運動量を$X,P$としましょう。要するに

$$\begin{align*} X&=X(x,p)\\ P&=P(x,p) \end{align*}$$

$X,P$が$x,p$の関数になる場合です。これは$x,p$からの変換になります。例としてこれまで使った$x_{n+1},p_{n+1}$は

$$\begin{align*} x_{n+1}&=x_{n+1}(x_{n},p_{n})\\ p_{n+1}&=p_{n+1}(x_{n},p_{n}) \end{align*}$$

と書けます。この$X,P$についても

$$\begin{align*} \frac{dX}{dt}&=\frac{\partial H}{\partial P}\\ \frac{dP}{dt}&=-\frac{\partial H}{\partial X} \end{align*}$$

と同様にハミルトンの正準方程式が成り立つことを要求します^[4][5]。この正準方程式を成り立たせる変換は正準変換と呼ばれます。正準変換が成り立つためには

$$\{X,P\}_{x,p}:=\frac{\partial X}{\partial x}\frac{\partial P}{\partial p}-\frac{\partial P}{\partial x}\frac{\partial X}{\partial p}=1$$

を満たす必要があります。$\{\cdot,\cdot\}_{x,p}$はポアッソン括弧といいます。ここでシンプレクティック数値積分法においてポアッソン括弧を$1$で保つか確認しましょう。

$$\begin{align*} \frac{\partial x_{n+1}}{\partial x_{n}}&=\frac{\partial}{\partial x_{n}}\left(x_{n}+\epsilon p_{n}\right)=1\\ \frac{\partial x_{n+1}}{\partial p_{n}}&=\frac{\partial}{\partial p_{n}}\left(x_{n}+\epsilon p_{n}\right)=\epsilon\\ \frac{\partial p_{n+1}}{\partial x_{n}}&=\frac{\partial}{\partial x_{n}}\left(-\epsilon x_{n} + (1-\epsilon^{2}) p_{n}\right)=-\epsilon\\ \frac{\partial p_{n+1}}{\partial p_{n}}&=\frac{\partial}{\partial p_{n}}\left(-\epsilon x_{n} + (1-\epsilon^{2}) p_{n}\right)=1-\epsilon^{2} \end{align*}$$

より、

$$\begin{align*} \{x_{n+1},p_{n+1}\}_{x_{n},p_{n}}&=\frac{\partial x_{n+1}}{\partial x_{n}}\frac{\partial p_{n+1}}{\partial p_{n}}-\frac{\partial p_{n+1}}{\partial x_{n}}\frac{\partial x_{n+1}}{\partial p_{n}}\\ &=1\cdot\left(1-\epsilon^{2}\right)-\left(-\epsilon\right)\epsilon\\ &=1 \end{align*}$$

と$1$でポアッソン括弧を保つことが確認できました！これがシンプレクティック数値積分法の原理です。因みにオイラー法では$1$になりません。

"シンプレクティック"って？

シンプレクティック数値積分法の"シンプレクティック"の由来を説明しようと思います。

相空間

これまで位置と運動量を変数として話を展開させていきました。というのも解析力学（統計力学や量子力学もそうですが）位置と運動量で論じることが多いです。従って位置と運動量からなる空間を作って置けば、その空間内で物理状態を示すことができるようになります。この空間を相空間と呼びます。

相空間の次元は対象とする系によって変わるのですが、例えば今回のような調和振動子の場合一つの質点しか考えないので位置$x$と運動量$p$がなす二次元の空間$(x,p)$を考えることになります。その空間に一つ点$(x_{0},p_{0})$を打っておけば、位置や運動量のみならずエネルギー等も決定するのでそれだけで系の状態を記述することができます。

相空間の図：$(x_{0},p_{0})$で点を打つとエネルギーも$\frac{1}{2}(p_{0}^{2}+x_{0}^{2})$で決定する。

ここで調和振動子が常に$E$というエネルギーを取ることを考えると、調和振動子のハミルトニアンから

$$E=\frac{1}{2}(p^{2} + x^{2})$$

であるので、これを相空間に反映させると、

相空間：ピンク実線: 解析解

半径$\sqrt{2E}$の円が出来上がります。まぎれもなくこの円は解として導かれる

$$(x(t),p(t))=(x_{0}\cos{t},-x_{0}\sin{t})$$

の軌跡です。質点は必ずこの円に沿って運動するので、時間が経過してもエネルギーが増え続ける（減り続ける）ことはないことが明らかになります。次にシンプレクティック数値積分法で導かれる結果を相空間に反映してみましょう。

相空間：ピンク実線: 解析解、緑破線：シンプレクティック数値積分法

しっかりと円に沿って描かれていることが分かりますね。分かりづらいですがこの時$x$軸に対してマイナス45度の傾きの方向に長軸を持つ楕円になります。シンプレクティック数値積分法においても軌跡が楕円で閉じるのでエネルギーが増え続けないことが分かります。最後にオイラー法を相空間に反映します。

相空間：ピンク実線: 解析解、青破線: オイラー法、緑破線：シンプレクティック数値積分法

オイラー法の場合は円や楕円で閉じることなく外へ発散していくため、エネルギーが保存しないことが分かります。

シンプレクティックと正準変換

唐突ですが$X=X(x,p)$と$P=P(x,p)$の全微分を計算しましょう。

$$\begin{align*} dX&=\frac{\partial X}{\partial x}dx + \frac{\partial X}{\partial p}dp\\ dP&=\frac{\partial P}{\partial x}dx + \frac{\partial P}{\partial p}dp \end{align*}$$

これを行列の形式に書き直すと、

$$\begin{pmatrix} dX \\ dP \end{pmatrix}= \begin{pmatrix}\displaystyle \frac{\partial X}{\partial x}&\displaystyle \frac{\partial X}{\partial p}\displaystyle \\[1.2 em] \displaystyle \frac{\partial P}{\partial x}&\displaystyle\frac{\partial P}{\partial p} \end{pmatrix} \begin{pmatrix} dx \\ dp \end{pmatrix}$$

となります。ここで次のヤコビ行列が定義できます。

$$\frac{\partial (X,P)}{\partial (x,p)}:=\begin{pmatrix}\displaystyle \frac{\partial X}{\partial x}&\displaystyle \frac{\partial X}{\partial p}\displaystyle \\[1.2 em] \displaystyle \frac{\partial P}{\partial x}&\displaystyle\frac{\partial P}{\partial p} \end{pmatrix}$$

そしてヤコビ行列の行列式（ヤコビアン）は面積$dxdp,dXdP$間の比率を表しており、

$$dXdP=\text{det}\left(\frac{\partial (X,P)}{\partial (x,p)}\right)dxdp=\left(\frac{\partial X}{\partial x}\frac{\partial P}{\partial p}-\frac{\partial P}{\partial x}\frac{\partial X}{\partial p}\right)dxdp$$

の関係を与えます。ここであれ？と思うのがヤコビアンがまさしくポアッソン括弧の定義になっているということです。つまり相空間上で正準変換が成り立つというのは

$$\text{det}\left(\frac{\partial (X,P)}{\partial (x,p)}\right)=\{X,P\}_{x,p}=1$$

$dxdp=dXdP$と変換後の面積を変化させないことを言います。正準変換が成り立ったり相空間上の面積が保存されるような変換をシンプレクティック変換と呼んだりします。

まとめ

最後にこれまでの話をまとめましょう。

• シンプレクティック数値積分法はシンプレクティック変換を原理とする数値解析法である。
• シンプレクティック数値積分法は厳密にエネルギーを保存しない。

ということです。最後に紹介したシンプレクティックという概念は微分形式を使うと正準変換にしろ面積保存にしろ次のようにたった一つの式で記述することができます。

$$dx\land dp=dX\land dP$$

この$\land$という記号はウェッジ積や外積と呼び、形式的にヤコビアンを導くことができる演算です。先ほど$dx,dp$からなるベクトルを作り、いきなり$dxdp$とかけ合わせましたが、それを演算という形で面積という意味も含ませてヤコビアンを導出することができます。

Appendix

参考文献

十河 清（著）「解析力学」株式会社日本評論社

ソースコード

オイラー法とシンプレクティック数値積分法を比較できるシミュレーションコードです。Pythonで書いています。使用する場合はnumpyとmatplotlibもインストールして下さい。使用する場合は自己責任でお願いします。

シミュレーションコード

Python

import numpy as np
import matplotlib.pyplot as plt

x_0 = 1 # 初期位置
epsilon = 0.05 # ステップ幅
t_max = 20 # シミュレーション期間

def hamiltonian(x, p):
    return (x**2 + p**2) / 2

t = np.arange(0, t_max, epsilon)
N = len(t)

# a:解析解、e:オイラー法、s:シンプレクティック数値積分法
# 初期化処理
x_a = x_0*np.cos(t)
p_a = -x_0*np.sin(t)
x_e = np.zeros(N)
p_e = np.zeros(N)
x_s = np.zeros(N)
p_s = np.zeros(N)

# 初期値設定
x_e[0] = x_0
x_s[0] = x_0

# エネルギー初期化処理
E_a = hamiltonian(x_a, p_a)
E_e = np.full(N, hamiltonian(x_e[0], p_e[0]))
E_s = np.full(N, hamiltonian(x_s[0], p_s[0]))

for n in range(N-1):
    # オイラー法
    x_e[n+1] = x_e[n] + epsilon*p_e[n]
    p_e[n+1] = -epsilon*x_e[n] + p_e[n]
    E_e[n+1] = hamiltonian(x_e[n+1], p_e[n+1])

    # シンプレクティック数値積分法
    x_s[n+1] = x_s[n] + epsilon*p_s[n]
    p_s[n+1] = -epsilon*x_s[n+1] + p_s[n]
    E_s[n+1] = hamiltonian(x_s[n+1], p_s[n+1])

analytical_color = [1,0.3,0.3]
euler_color = [0,0,1]
symplectic_color = [0,1,0]

# 位置xシミュレーション結果
plt.subplots(figsize=(10, 4))
plt.xlabel('time (sec)')
plt.ylabel('x')
plt.plot(t, x_a, color=analytical_color, linewidth=6, label='analytical')
plt.plot(t, x_e, color=euler_color, linewidth=3,linestyle="dashed", label='euler')
plt.plot(t, x_s, color=symplectic_color, linewidth=3, linestyle="dashed" ,label='symplectic')
plt.legend()
plt.grid()

# 運動量pシミュレーション結果
plt.subplots(figsize=(10, 4))
plt.xlabel('time (sec)')
plt.ylabel('p')
plt.plot(t, p_a, color=analytical_color, linewidth=6, label='analytical')
plt.plot(t, p_e, color=euler_color, linewidth=3,linestyle="dashed", label='euler')
plt.plot(t, p_s, color=symplectic_color, linewidth=3, linestyle="dashed" ,label='symplectic')
plt.legend()
plt.grid()

# エネルギー評価
ymax = np.ceil(np.amax(E_e))
plt.subplots(figsize=(8, 5))
plt.xlabel('time (sec)')
plt.ylabel(r'Energy $\frac{1}{2}(x^{2} + p^{2})$')
plt.plot(t, E_a, color=analytical_color, linewidth=6, label='analytical')
plt.plot(t, E_e, color=euler_color, linewidth=3,linestyle="dashed", label='euler')
plt.plot(t, E_s, color=symplectic_color, linewidth=3, linestyle="dashed" ,label='symplectic')
plt.ylim(0, ymax)
plt.legend()
plt.grid()

# 相空間
plt.subplots(figsize=(5, 5))
plt.plot(x_a, p_a, color=analytical_color, linewidth=6, label='analytical')
plt.plot(x_e, p_e, color=euler_color, linewidth=3,linestyle="dashed", label='euler')
plt.plot(x_s, p_s, color=symplectic_color, linewidth=3, linestyle="dashed" ,label='symplectic')
plt.xlabel('x')
plt.ylabel('p')
plt.xlim(-1.6, 1.6)
plt.ylim(-1.6, 1.6)
plt.legend()
plt.grid()

plt.show()

脚注

1. 物理的に考慮すべき対象のこと。単に系と呼ぶことにします。

2. この式変形は少し骨が折れます。めんどくさい人はmathematicaやsympyなどをつかって楽するのもありです。

3. 正確には一般化座標、一般化運動量と言います。

4. 物理には最小作用の原理というものがあります。簡単に説明すると力学系（時間変化する物理系）は無駄な運動をしないという物理のルールで、例えばボールを落とした時、突然曲がったりすることなく必ず一直線で落下します。ハミルトンの正準方程式が最小作用の原理から導かれることを踏まえると、正準方程式が成り立たないならば最小作用の原理が満たされないということになります。

5. もっと一般化するとこの正準方程式は$H$ではなく$H'$とか別のハミルトニアンで記述するのが良いでしょう。このハミルトニアン$H'$は母関数というもの使って探します。