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

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

・高校レベルの古典力学の知識。

・全微分あたりまでの微積分の知識と行列の計算

その他

・fffの定義をf:=f(x)f:=f(x)f:=f(x)と書く。

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



上図のバネをびよーんって伸ばしてから離して振動する質点Pのttt秒後の位置x=x(t)x=x(t)x=x(t)と運動量p=p(t)p=p(t)p=p(t)を知りたいです。これは調和振動子と呼ばれる物理系[1]になります。このとき質点の進行方向に対して逆向きにバネの力が働くので運動方程式は次のようになります。
md2xdt2=−kx
m\frac{d^{2}x}{dt^{2}}=-kx
mdt2d2x​=−kxここでm,km,km,kは質点の質量とバネ定数になりますが簡単のため両方ともm=1,k=1m=1,k=1m=1,k=1としてしまいましょう。そして運動量pppは質量mmmと速度dxdt\frac{dx}{dt}dtdx​の積であり、
dxdt=p
\frac{dx}{dt}=p
dtdx​=pと書けます。そうすると運動方程式は
dpdt=−x
\frac{dp}{dt}=-x
dtdp​=−xとなります。このdxdt\frac{dx}{dt}dtdx​とdpdt\frac{dp}{dt}dtdp​に関する上二つの微分方程式から次の解析解が得られます。
x(t)=x0cos⁡t,  p(t)=−x0sin⁡t
x(t)=x_{0}\cos t,\;p(t)=-x_{0}\sin t
x(t)=x0​cost,p(t)=−x0​sintはじめ質点Pはx(0)=x0x(0)=x_{0}x(0)=x0​に位置し、完全に止まっていることp(0)=0p(0)=0p(0)=0を想定しました。この解からx(t)x(t)x(t)もp(t)p(t)p(t)も振動することが分かりました。
次に微分方程式の解き方は知らないけど微分は知っているAさんが現れました。Aさんは非常に近い二点間の変化率が微分だと理解していたので、先ほどの式をϵ≈0\epsilon\approx0ϵ≈0を使って近似してみました。
dxdt≈xn+1−xnϵ=pn,  dpdt≈pn+1−pnϵ=−xn
\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}
dtdx​≈ϵxn+1​−xn​​=pn​,dtdp​≈ϵpn+1​−pn​​=−xn​この式から次の漸化式が導けます。
xn+1=xn+ϵpnpn+1=−ϵxn+pn
\begin{align*}
x_{n+1}&=x_{n}+\epsilon p_{n}\\
p_{n+1}&=-\epsilon x_{n}+p_{n}
\end{align*}
xn+1​pn+1​​=xn​+ϵpn​=−ϵxn​+pn​​t=0t=0t=0に対応するnnnをn=0n=0n=0すれば、上式で再帰的にnnnステップ目の位置と運動量xn,pnx_{n},p_{n}xn​,pn​が求まります。微分の原理に基づき時間をϵ\epsilonϵで離散化して再帰的に解く手法をオイラー法と言います。せっかくなのでオイラー法から導かれる位置と運動量の時間変化を見てましょう。ついでに解析解とも比較しておきます。まずは位置xxxです。



位置xxxの時間変化：ピンク実線: 解析解、青破線: オイラー法
横軸が時間、縦軸が位置です。またϵ=0.05,x0=1\epsilon=0.05,x_{0}=1ϵ=0.05,x0​=1で設定しました。時間が経つにつれ解析解から遠のいていくのが分かります。次に運動量です。



運動量pppの時間変化：ピンク実線: 解析解、青破線: オイラー法
運動量も同様に解析解から遠のいていきます。

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

 エネルギーとハミルトニアン調和振動子におけるエネルギーは、質点の運動エネルギーKKKとバネの弾性エネルギーVVVです。それぞれ次のように書けます。
K=12mv2=p22m,  V=12kx2
K=\frac{1}{2}mv^{2}=\frac{p^{2}}{2m},\;V=\frac{1}{2}kx^{2}
K=21​mv2=2mp2​,V=21​kx2v=dxdtv=\frac{dx}{dt}v=dtdx​は質点の速度です。この系において摩擦などの外的要因がないとすれば、全エネルギー
E=K+V
E=K+V
E=K+Vは保存します。これを少し一般化しましょう。KKKはよく見ると運動量pppの関数です。またVVVはxxxの関数です。そうするとその和はxxxとpppの関数になります。なので新たに
H=H(x,p):=K+V=p22m+12kx2
H=H(x,p):=K+V=\frac{p^{2}}{2m}+\frac{1}{2}kx^{2}
H=H(x,p):=K+V=2mp2​+21​kx2を定義しましょう。このHHHはハミルトニアンと言います。ここでm=1,k=1m=1,k=1m=1,k=1の仮定を思い出すとこの調和振動子のハミルトニアンは
H=12(p2+x2)
H=\frac{1}{2}\left(p^{2}+x^{2}\right)
H=21​(p2+x2)と書けます。このハミルトニアンにxn,pnx_{n},p_{n}xn​,pn​とxn+1,pn+1x_{n+1},p_{n+1}xn+1​,pn+1​を代入しnnnとn+1n+1n+1ステップの関係を調べると、
12(pn+12+xn+12)=12(pn2+xn2)(1+ϵ2)
\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})
21​(pn+12​+xn+12​)=21​(pn2​+xn2​)(1+ϵ2)ステップを一つ進めるごとにエネルギーが1+ϵ21+\epsilon^{2}1+ϵ2倍に蓄積することが確認できます。つまりオイラー法では本来満たしてほしいエネルギーの保存が原理的に満たされないわけです。

 オイラー法を工夫するこれらの議論からステップを進めてもエネルギーが蓄積しない手法を考えましょう。結論から言うと次のように改良されます。
xn+1=xn+ϵpnpn+1=−ϵxn+1+pn(=−ϵxn+(1−ϵ2)pn)
\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*}
xn+1​pn+1​​=xn​+ϵpn​=−ϵxn+1​+pn​(=−ϵxn​+(1−ϵ2)pn​)​オイラー法と何が違うかというとpn+1p_{n+1}pn+1​への更新の際、xnx_{n}xn​ではなくxn+1x_{n+1}xn+1​を用いている点です。このとき確かにエネルギーは保存されるのですが、実際に保存するハミルトニアンは次のような形になります。
H=12(p2+x2)+ϵ2px
H=\frac{1}{2}\left(p^{2}+x^{2}\right)+\frac{\epsilon}{2}px
H=21​(p2+x2)+2ϵ​pxϵ2px\frac{\epsilon}{2}px2ϵ​pxという余計な項が含まれていますが、試しにnnnとn+1n+1n+1ステップ目の関係を調べると
12(pn+12+xn+12)+ϵ2pn+1xn+1=12(pn2+xn2)+ϵ2pnxn
\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}
21​(pn+12​+xn+12​)+2ϵ​pn+1​xn+1​=21​(pn2​+xn2​)+2ϵ​pn​xn​と見事にエネルギーが保存されます[2]。とりあえずこの改良を通じてオイラー法と比較しつつシミュレーションした結果を示します。まずは位置xxxです。



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

次に運動量pppです。



運動量pppの時間変化：ピンク実線: 解析解、青破線: オイラー法、緑破線：シンプレクティック数値積分法
正しく解析解に沿ってシミュレーションできていることが確認できました！この数値解析法をシンプレクティック数値積分法と言います。
ここで注意しなければいけないのは「じゃあシンプレクティック数値積分法はエネルギーを保存するシミュレーション方法なのね！」と理解してしまうことです。あくまでも調和振動子本来のハミルトニアンは12(p2+x2)\frac{1}{2}(p^{2}+x^{2})21​(p2+x2)なので、これで時間ごとのエネルギーを評価すると厳密に保存せず一定の周期で振動していることが確認できます。



時間ごとのエネルギーの大きさ

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

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

 ハミルトンの正準方程式と正準変換解析力学で学ぶ物体の運動を表す式には、ラグランジュの運動方程式とハミルトンの正準方程式があります。どちらもとても重要なのですがハミルトンの正準方程式を紹介します。ハミルトンの正準方程式は次のように定義されます。
dxdt=∂H∂pdpdt=−∂H∂x
\begin{align*}
\frac{dx}{dt}&=\frac{\partial H}{\partial p}\\
\frac{dp}{dt}&=-\frac{\partial H}{\partial x}
\end{align*}
dtdx​dtdp​​=∂p∂H​=−∂x∂H​​HHHは例のごとくハミルトニアンです。x,px,px,pも位置と運動量になります[3]。このハミルトニアンに調和振動子のものを代入し計算してみてください。記事のはじめに導いたdxdt=p,dpdt=−x\frac{dx}{dt}=p,\frac{dp}{dt}=-xdtdx​=p,dtdp​=−xが導出されます。さてここでx,px,px,pから導かれる新たな位置と運動量をX,PX,PX,Pとしましょう。要するに
X=X(x,p)P=P(x,p)
\begin{align*}
X&=X(x,p)\\
P&=P(x,p)
\end{align*}
XP​=X(x,p)=P(x,p)​X,PX,PX,Pがx,px,px,pの関数になる場合です。これはx,px,px,pからの変換になります。例としてこれまで使ったxn+1,pn+1x_{n+1},p_{n+1}xn+1​,pn+1​は
xn+1=xn+1(xn,pn)pn+1=pn+1(xn,pn)
\begin{align*}
x_{n+1}&=x_{n+1}(x_{n},p_{n})\\
p_{n+1}&=p_{n+1}(x_{n},p_{n})
\end{align*}
xn+1​pn+1​​=xn+1​(xn​,pn​)=pn+1​(xn​,pn​)​と書けます。このX,PX,PX,Pについても
dXdt=∂H∂PdPdt=−∂H∂X
\begin{align*}
\frac{dX}{dt}&=\frac{\partial H}{\partial P}\\
\frac{dP}{dt}&=-\frac{\partial H}{\partial X}
\end{align*}
dtdX​dtdP​​=∂P∂H​=−∂X∂H​​と同様にハミルトンの正準方程式が成り立つことを要求します[4][5]。この正準方程式を成り立たせる変換は正準変換と呼ばれます。正準変換が成り立つためには
{X,P}x,p:=∂X∂x∂P∂p−∂P∂x∂X∂p=1
\{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
{X,P}x,p​:=∂x∂X​∂p∂P​−∂x∂P​∂p∂X​=1を満たす必要があります。{⋅,⋅}x,p\{\cdot,\cdot\}_{x,p}{⋅,⋅}x,p​はポアッソン括弧といいます。ここでシンプレクティック数値積分法においてポアッソン括弧を111で保つか確認しましょう。
∂xn+1∂xn=∂∂xn(xn+ϵpn)=1∂xn+1∂pn=∂∂pn(xn+ϵpn)=ϵ∂pn+1∂xn=∂∂xn(−ϵxn+(1−ϵ2)pn)=−ϵ∂pn+1∂pn=∂∂pn(−ϵxn+(1−ϵ2)pn)=1−ϵ2
\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*}
∂xn​∂xn+1​​∂pn​∂xn+1​​∂xn​∂pn+1​​∂pn​∂pn+1​​​=∂xn​∂​(xn​+ϵpn​)=1=∂pn​∂​(xn​+ϵpn​)=ϵ=∂xn​∂​(−ϵxn​+(1−ϵ2)pn​)=−ϵ=∂pn​∂​(−ϵxn​+(1−ϵ2)pn​)=1−ϵ2​より、
{xn+1,pn+1}xn,pn=∂xn+1∂xn∂pn+1∂pn−∂pn+1∂xn∂xn+1∂pn=1⋅(1−ϵ2)−(−ϵ)ϵ=1
\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*}
{xn+1​,pn+1​}xn​,pn​​​=∂xn​∂xn+1​​∂pn​∂pn+1​​−∂xn​∂pn+1​​∂pn​∂xn+1​​=1⋅(1−ϵ2)−(−ϵ)ϵ=1​と111でポアッソン括弧を保つことが確認できました！これがシンプレクティック数値積分法の原理です。因みにオイラー法では111になりません。

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

 相空間これまで位置と運動量を変数として話を展開させていきました。というのも解析力学（統計力学や量子力学もそうですが）位置と運動量で論じることが多いです。従って位置と運動量からなる空間を作って置けば、その空間内で物理状態を示すことができるようになります。この空間を相空間と呼びます。
相空間の次元は対象とする系によって変わるのですが、例えば今回のような調和振動子の場合一つの質点しか考えないので位置xxxと運動量pppがなす二次元の空間(x,p)(x,p)(x,p)を考えることになります。その空間に一つ点(x0,p0)(x_{0},p_{0})(x0​,p0​)を打っておけば、位置や運動量のみならずエネルギー等も決定するのでそれだけで系の状態を記述することができます。


相空間の図：(x0,p0)(x_{0},p_{0})(x0​,p0​)で点を打つとエネルギーも12(p02+x02)\frac{1}{2}(p_{0}^{2}+x_{0}^{2})21​(p02​+x02​)で決定する。
ここで調和振動子が常にEEEというエネルギーを取ることを考えると、調和振動子のハミルトニアンから
E=12(p2+x2)
E=\frac{1}{2}(p^{2} + x^{2})
E=21​(p2+x2)であるので、これを相空間に反映させると、


相空間：ピンク実線: 解析解
半径2E\sqrt{2E}2E​の円が出来上がります。まぎれもなくこの円は解として導かれる
(x(t),p(t))=(x0cos⁡t,−x0sin⁡t)
(x(t),p(t))=(x_{0}\cos{t},-x_{0}\sin{t})
(x(t),p(t))=(x0​cost,−x0​sint)の軌跡です。質点は必ずこの円に沿って運動するので、時間が経過してもエネルギーが増え続ける（減り続ける）ことはないことが明らかになります。次にシンプレクティック数値積分法で導かれる結果を相空間に反映してみましょう。


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


相空間：ピンク実線: 解析解、青破線: オイラー法、緑破線：シンプレクティック数値積分法
オイラー法の場合は円や楕円で閉じることなく外へ発散していくため、エネルギーが保存しないことが分かります。

 シンプレクティックと正準変換唐突ですがX=X(x,p)X=X(x,p)X=X(x,p)とP=P(x,p)P=P(x,p)P=P(x,p)の全微分を計算しましょう。
dX=∂X∂xdx+∂X∂pdpdP=∂P∂xdx+∂P∂pdp
\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*}
dXdP​=∂x∂X​dx+∂p∂X​dp=∂x∂P​dx+∂p∂P​dp​これを行列の形式に書き直すと、
(dXdP)=(∂X∂x∂X∂p∂P∂x∂P∂p)(dxdp)
\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}
(dXdP​)=​∂x∂X​∂x∂P​​∂p∂X​∂p∂P​​​(dxdp​)となります。ここで次のヤコビ行列が定義できます。
∂(X,P)∂(x,p):=(∂X∂x∂X∂p∂P∂x∂P∂p)
\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}
∂(x,p)∂(X,P)​:=​∂x∂X​∂x∂P​​∂p∂X​∂p∂P​​​そしてヤコビ行列の行列式（ヤコビアン）は面積dxdp,dXdPdxdp,dXdPdxdp,dXdP間の比率を表しており、
dXdP=det(∂(X,P)∂(x,p))dxdp=(∂X∂x∂P∂p−∂P∂x∂X∂p)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
dXdP=det(∂(x,p)∂(X,P)​)dxdp=(∂x∂X​∂p∂P​−∂x∂P​∂p∂X​)dxdpの関係を与えます。ここであれ？と思うのがヤコビアンがまさしくポアッソン括弧の定義になっているということです。つまり相空間上で正準変換が成り立つというのは
det(∂(X,P)∂(x,p))={X,P}x,p=1
\text{det}\left(\frac{\partial (X,P)}{\partial (x,p)}\right)=\{X,P\}_{x,p}=1
det(∂(x,p)∂(X,P)​)={X,P}x,p​=1dxdp=dXdPdxdp=dXdPdxdp=dXdPと変換後の面積を変化させないことを言います。正準変換が成り立ったり相空間上の面積が保存されるような変換をシンプレクティック変換と呼んだりします。

 まとめ最後にこれまでの話をまとめましょう。
シンプレクティック数値積分法はシンプレクティック変換を原理とする数値解析法である。
シンプレクティック数値積分法は厳密にエネルギーを保存しない。
ということです。最後に紹介したシンプレクティックという概念は微分形式を使うと正準変換にしろ面積保存にしろ次のようにたった一つの式で記述することができます。
dx∧dp=dX∧dP
dx\land dp=dX\land dP
dx∧dp=dX∧dPこの∧\land∧という記号はウェッジ積や外積と呼び、形式的にヤコビアンを導くことができる演算です。先ほどdx,dpdx,dpdx,dpからなるベクトルを作り、いきなりdxdpdxdpdxdpとかけ合わせましたが、それを演算という形で面積という意味も含ませてヤコビアンを導出することができます。

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

 ソースコードオイラー法とシンプレクティック数値積分法を比較できるシミュレーションコードです。実行する場合は自己責任でお願いします。
https://github.com/torataro-tiger/symplectic_integrator/blob/main/main.py

脚注
物理的に考慮すべき対象のこと。単に系と呼ぶことにします。 ↩︎
この式変形は少し骨が折れます。めんどくさい人はmathematicaやsympyなどをつかって楽するのもありです。 ↩︎
正確には一般化座標、一般化運動量と言います。 ↩︎
物理には最小作用の原理というものがあります。簡単に説明すると力学系（時間変化する物理系）は無駄な運動をしないという物理のルールで、例えばボールを落とした時、突然曲がったりすることなく必ず一直線で落下します。ハミルトンの正準方程式が最小作用の原理から導かれることを踏まえると、正準方程式が成り立たないならば最小作用の原理が満たされないということになります。 ↩︎
もっと一般化するとこの正準方程式はHHHではなくH′H'H′とか別のハミルトニアンで記述するのが良いでしょう。このハミルトニアンH′H'H′は母関数というもの使って探します。 ↩︎