はじめに
数値解析学を学ぶと微分方程式の解法として一番簡単なオイラー法が紹介されると思います。そして実際に実装してみると確かに解くことはできるけど、誤差がどんどん大きくなるところに目を瞑っていられなくなります。そこで調和振動子などの物理系に対し、安定して長期間のシミュレーションが行えるシンプレクティック数値積分法を紹介したいと思います。
問題提起
よく古典力学の教科書とかで見る次の状況を考えましょう。
上図のバネをびよーんって伸ばしてから離して振動する質点Pのt秒後の位置x=x(t)と運動量p=p(t)を知りたいです。これは調和振動子と呼ばれる物理系になります。このとき質点の進行方向に対して逆向きにバネの力が働くので運動方程式は次のようになります。
m\frac{d^{2}x}{dt^{2}}=-kx
ここでm,kは質点の質量とバネ定数になりますが簡単のため両方ともm=1,k=1としてしまいましょう。そして運動量pは質量mと速度\frac{dx}{dt}の積であり、
と書けます。そうすると運動方程式は
となります。この\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}は質点の速度です。この系において摩擦などの外的要因がないとすれば、全エネルギー
は保存します。これを少し一般化しましょう。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}
と見事にエネルギーが保存されます。とりあえずこの改良を通じてオイラー法と比較しつつシミュレーションした結果を示します。まずは位置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も位置と運動量になります。このハミルトニアンに調和振動子のものを代入し計算してみてください。記事のはじめに導いた\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*}
と同様にハミルトンの正準方程式が成り立つことを要求します。この正準方程式を成り立たせる変換は正準変換と呼ばれます。正準変換が成り立つためには
\{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と変換後の面積を変化させないことを言います。正準変換が成り立ったり相空間上の面積が保存されるような変換をシンプレクティック変換と呼んだりします。
まとめ
最後にこれまでの話をまとめましょう。
- シンプレクティック数値積分法はシンプレクティック変換を原理とする数値解析法である。
- シンプレクティック数値積分法は厳密にエネルギーを保存しない。
ということです。最後に紹介したシンプレクティックという概念は微分形式を使うと正準変換にしろ面積保存にしろ次のようにたった一つの式で記述することができます。
この\landという記号はウェッジ積や外積と呼び、形式的にヤコビアンを導くことができる演算です。先ほどdx,dpからなるベクトルを作り、いきなりdxdpとかけ合わせましたが、それを演算という形で面積という意味も含ませてヤコビアンを導出することができます。
Appendix
参考文献
十河 清（著）「解析力学」株式会社日本評論社
ソースコード
オイラー法とシンプレクティック数値積分法を比較できるシミュレーションコードです。実行する場合は自己責任でお願いします。
https://github.com/torataro-tiger/symplectic_integrator/blob/main/main.py
Discussion