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

シンプレクティック数値積分法 (シンプレクティックすうちせきぶんほう, symplectic integrator) とは、正準力学系の運動方程式に特化した常微分方程式の数値解法のことをいう。系のシンプレクティック形式およびハミルトニアンを保存するため、ルンゲ＝クッタ法のような汎用の数値積分法に比べて良い性質を示す。このために天体力学などの分野で採用されている[1]。

概要

オイラー法、ルンゲ＝クッタ法とシンプレクティック積分子による調和振動子の数値解のエネルギー誤差の比較。横軸は周期で規格化した時間、縦軸は数値解のエネルギーの真のエネルギーに対する相対誤差。すべての数値解で時間刻み幅は同一である。オイラー法 (Euler) およびルンゲ＝クッタ法 (RK4) では誤差が単調に増加する一方、シンプレクティック積分法 (Symp1-4) では誤差の増大が生じない。

正準力学系において、 ( q i , p i ) {\displaystyle (q_{i},p_{i})} を正準変数、 H = H ( q i , p i ) {\displaystyle H=H(q_{i},p_{i})} をハミルトニアンとするとき、その運動方程式はハミルトンの正準方程式

d q i d t = ∂ H ∂ p i ( q , p ) ,     d p i d t = − ∂ H ∂ q i ( q , p ) {\displaystyle {\frac {dq_{i}}{dt}}={\frac {\partial H}{\partial p_{i}}}(q,p),\ \ {\frac {dp_{i}}{dt}}=-{\frac {\partial H}{\partial q_{i}}}(q,p)}

である。これらの運動方程式の解 ( q i ( t ) , p i ( t ) ) {\displaystyle \left(q_{i}(t),p_{i}(t)\right)} は一般に次の性質を満足する[2]。

  • シンプレクティック形式 ω = ∑ i d p i ∧ d q i {\displaystyle \omega =\sum _{i}dp_{i}\wedge dq_{i}} を不変に保つ。
  • ハミルトニアン H ( q i ( t ) , p i ( t ) ) {\displaystyle H\left(q_{i}(t),p_{i}(t)\right)} を不変に保つ。

ところが、正準方程式を数値的に解くためにルンゲ＝クッタ法のような汎用の数値積分アルゴリズムを適用すると、一般に数値解においてこれらの性質が破れ、長時間の積分によりエネルギーが保存しないなどの非物理的な結果を生じ得る。シンプレクティック積分法は厳密にシンプレクティック写像であるような（すなわちシンプレクティック形式を保存する）数値積分アルゴリズムであり[2]、ハミルトン力学系の数値解析手法としてより優れた性質を示す。例えば調和振動子

H = 1 2 p 2 + 1 2 q 2 {\displaystyle H={\frac {1}{2}}p^{2}+{\frac {1}{2}}q^{2}}

に2次のシンプレクティック積分法を適用すると、時間刻み幅を h {\displaystyle h} として、真のハミルトニアンの代わりに

H ~ = 1 2 p 2 + 1 2 q 2 − h 2 4 q 2 {\displaystyle {\tilde {H}}={\frac {1}{2}}p^{2}+{\frac {1}{2}}q^{2}-{\frac {h^{2}}{4}}q^{2}}

を保存するため、数値解はある楕円軌道の上に留まり、エネルギーの単調な増加または減少を生じない[3]。さらに、偶数次のシンプレクティック積分子は時間反転対称性を持つという利点も存在する[4]。

一方で、時間刻み幅を動的に変更する適応時間刻みを単純にシンプレクティック積分子に適用するとハミルトニアンの保存が破れることが知られている[2]。

アルゴリズム

ハミルトニアン H {\displaystyle H} が2つの可積分ハミルトニアン H A {\displaystyle H_{A}} , H B {\displaystyle H_{B}} の和であると仮定する。

H ( q , p ) = H A ( q , p ) + H B ( q , p ) {\displaystyle H(q,p)=H_{A}(q,p)+H_{B}(q,p)}

例えばポテンシャル V {\displaystyle V} 中の質量 m {\displaystyle m} の粒子という系の場合 H A = 1 2 m p 2 {\displaystyle H_{A}={\frac {1}{2m}}p^{2}} , H B = V ( q ) {\displaystyle H_{B}=V(q)} であり、この仮定を満足する。 H {\displaystyle H} , H A {\displaystyle H_{A}} , H B {\displaystyle H_{B}} に対応するハミルトンベクトル場をそれぞれ D {\displaystyle D} , A {\displaystyle A} , B {\displaystyle B} と書くとき

D = A + B {\displaystyle D=A+B}

が成立し、それぞれのハミルトンベクトル場に沿う時間 t {\displaystyle t} の発展すなわち指数写像を S ( t ) := exp ⁡ ( t D ) {\displaystyle S(t):=\exp \left(tD\right)} , exp ⁡ ( t A ) {\displaystyle \exp(tA)} , exp ⁡ ( t B ) {\displaystyle \exp(tB)} と書く（これはシンプレクティック写像である）。 H A {\displaystyle H_{A}} および H B {\displaystyle H_{B}} が可積分であるという仮定により exp ⁡ ( t A ) {\displaystyle \exp(tA)} , exp ⁡ ( t B ) {\displaystyle \exp(tB)} はその具体的な表示が既知である。ここでの問題は、真のハミルトニアンに関する指数写像

S ( t ) = exp ⁡ ( t D ) = exp ⁡ { t ( A + B ) } {\displaystyle S(t)=\exp \left(tD\right)=\exp \left\{t(A+B)\right\}}

を N {\displaystyle N} 区間の時間積分の集積 S ( t ) = ∏ i = 1 N S ( h ) {\displaystyle S(t)=\prod _{i=1}^{N}S(h)} (ここに h = t / N {\displaystyle h=t/N} とおいた) と書くとき、 S ( h ) {\displaystyle S(h)} を既知のシンプレクティック写像 exp ⁡ ( h A ) {\displaystyle \exp(hA)} , exp ⁡ ( h B ) {\displaystyle \exp(hB)} を用いて

exp ⁡ { h ( A + B ) } = ∏ k [ exp ⁡ ( c k h A ) exp ⁡ ( d k h B ) ] + O ( h n + 1 ) {\displaystyle \exp \left\{h(A+B)\right\}=\prod _{k}\left[\exp(c_{k}hA)\,\exp(d_{k}hB)\right]+{\mathcal {O}}\left(h^{n+1}\right)}

という形に近似することである[5]。この右辺が求める n {\displaystyle n} 次のシンプレクティック積分子であり、シンプレクティック性を満足することが保証され、ハミルトニアン

H ~ = H + h n H n + O { h n + 1 } {\displaystyle {\tilde {H}}=H+h^{n}H_{n}+{\mathcal {O}}\left\{h^{n+1}\right\}}

を保存する[2]。

1次のシンプレクティック法

次式で定義される変換 S 1 s t ( h ) {\displaystyle S_{\mathrm {1st} }(h)} は S ( h ) = S 1 s t ( h ) + O ( h 2 ) {\displaystyle S(h)=S_{\mathrm {1st} }(h)+{\mathcal {O}}(h^{2})} を満足する1次のシンプレクティック積分子である[6]。

S 1 s t ( h ) = exp ⁡ ( h A ) exp ⁡ ( h B ) {\displaystyle S_{\mathrm {1st} }(h)=\exp(hA)\,\exp(hB)}

特に H A = 1 2 p 2 {\displaystyle H_{A}={\frac {1}{2}}p^{2}} , H B = V ( q ) {\displaystyle H_{B}=V(q)} の場合、変換 exp ⁡ ( h B ) {\displaystyle \exp(hB)} は ( q , p ) ↦ ( q , p − h V ′ ( q ) ) {\displaystyle (q,p)\mapsto (q,p-hV'(q))} 、変換 exp ⁡ ( h A ) {\displaystyle \exp(hA)} は ( q , p ) ↦ ( q + h p , p ) {\displaystyle (q,p)\mapsto (q+hp,p)} と表示できるため、このスキーム S 1 s t = exp ⁡ ( h A ) exp ⁡ ( h B ) {\displaystyle S_{\mathrm {1st} }=\exp(hA)\exp(hB)} 全体としては

p n + 1 = p n − h d V d q ( q n ) {\displaystyle p_{n+1}=p_{n}-h{\frac {dV}{dq}}(q_{n})}
q n + 1 = q n + h p n + 1 {\displaystyle q_{n+1}=q_{n}+hp_{n+1}}

と表示できる。このスキームはオイラー法を修正したものとみなせるため、シンプレクティックオイラー法と呼ばれることもある[7]。

2次のシンプレクティック法

2次のシンプレクティック法は次式で与えられる[6]。なおこの積分スキームはリープ・フロッグ法あるいはベレの方法、Strörmer法など分野毎に異なった名称で知られている[8]。

S 2 n d ( h ) = exp ⁡ ( h 2 A ) exp ⁡ ( h B ) exp ⁡ ( h 2 A ) {\displaystyle S_{\mathrm {2nd} }(h)=\exp \left({\frac {h}{2}}A\right)\,\exp(hB)\,\exp \left({\frac {h}{2}}A\right)}

上述の H A = 1 2 p 2 {\displaystyle H_{A}={\frac {1}{2}}p^{2}} , H B = V ( q ) {\displaystyle H_{B}=V(q)} の場合にはこれは次のスキームに帰着する。

q n + 1 / 2 = q n + h 1 2 p n {\displaystyle q_{n+1/2}=q_{n}+h{\frac {1}{2}}p_{n}}
p n + 1 = p n − h d V d q ( q n + 1 / 2 ) {\displaystyle p_{n+1}=p_{n}-h{\frac {dV}{dq}}\left(q_{n+1/2}\right)}
q n + 1 = q n + 1 / 2 + h 1 2 p n + 1 {\displaystyle q_{n+1}=q_{n+1/2}+h{\frac {1}{2}}p_{n+1}}

4次のシンプレクティック法

4次のシンプレクティック積分子は、2次の積分子を異なる時間ステップで三度適用することにより得られる。

S 4 t h ( h ) = S 2 n d ( x 1 h ) S 2 n d ( x 0 h ) S 2 n d ( x 1 h ) {\displaystyle S_{\mathrm {4th} }(h)=S_{\mathrm {2nd} }(x_{1}h)\,S_{\mathrm {2nd} }(x_{0}h)\,S_{\mathrm {2nd} }(x_{1}h)}
x 1 = 1 2 − 2 3 ,     x 0 = 1 − 2 x 1 {\displaystyle x_{1}={\frac {1}{2-{\sqrt[{3}]{2}}}},\ \ x_{0}=1-2x_{1}}

これは Forest & Ruth (1990) によって導かれた[9] 後、Yoshida (1990) によって2次のシンプレクティック積分を三度適用したものに等しいことが指摘された[5]。

なお、より高次のシンプレクティック積分子の系統的な構成方法は Suzuki (1990)[10] および Yoshida (1990)[11]によって与えられている。Yoshida (1990) はベーカー・キャンベル・ハウスドルフの公式（英語版）を適用することにより高次の 2 n {\displaystyle 2n} 次シンプレクティック積分子を解析的に構成しているが, 次数が増大すると必要なステップ数 ( S 2 n d {\displaystyle S_{\mathrm {2nd} }} の数) が指数関数的に増大し効率が悪化することも指摘し、より効率的なシンプレクティック積分子をいくつか数値的に求めてもいる。

応用

オイラー法、ルンゲ＝クッタ法とシンプレクティック積分子によるケプラー問題 (離心率 e = 0.5 {\displaystyle e=0.5} ) の数値解のエネルギー誤差の比較。横軸は軌道周期で規格化した時間、縦軸は数値解のエネルギーの真のエネルギーに対する相対誤差。ルンゲ＝クッタ法などでは誤差が増加する一方、シンプレクティック積分法では誤差の増大が生じない。

ケプラー問題

重力場中の粒子の運動は天体力学や軌道力学、宇宙物理学での重要さからシンプレクティック積分法が適用される典型的な例となっている。例えばケプラー問題

H = 1 2 m p 2 − G M m | q | {\displaystyle H={\frac {1}{2m}}\mathbf {p} ^{2}-{\frac {GMm}{\left|\mathbf {q} \right|}}}

にシンプレクティック積分法を適用すると、オイラー法やルンゲ＝クッタ法では時間の経過とともに数値誤差が累積しエネルギーの誤差が増大するが、シンプレクティック積分法ではエネルギー誤差の累積は見られない[4]（右図）。

天体力学

中心天体からの重力が支配的であるような天体力学の典型的な問題では、系のハミルトニアンは中心天体によるケプラー問題の部分 H K e p {\displaystyle H_{\mathrm {Kep} }} と、天体間相互作用による摂動部分 H i n t {\displaystyle H_{\mathrm {int} }} に分割できる。

H = H K e p + H i n t {\displaystyle H=H_{\mathrm {Kep} }+H_{\mathrm {int} }}

それを踏まえて、ハミルトニアンをこの二項に分割してシンプレクティック積分子を適用する数値解析手法が Wisdom & Holman (1991[12], 1992[13]) および Kinoshita et al. (1991)[14] によって提案され、従来の手法より誤差が小さくなるなどの良好な性質を持つことが示された[1]。

脚注

  1. ^ a b 伊藤孝士. “シンプレクティク数値積分法の天体力学的応用”. 2020年7月7日閲覧。
  2. ^ a b c d 吉田春夫「可変時間ステップによるシンプレクティック数値解法(非線形可積分系による応用解析)」『数理解析研究所講究録』第889巻、京都大学数理解析研究所、1994年11月、70-76頁、CRID 1050001202298760960、hdl:2433/84362、ISSN 1880-2818、2024年1月12日閲覧。 
  3. ^ 牧野淳一郎, 福重俊幸, 小久保英一郎, 川井敦, 台坂博, 杉本大一郎 (2007年3月13日). “N体シミュレーション啓蟄の学校教科書”. 国立天文台. p. 51. 2020年5月24日閲覧。
  4. ^ a b Binney, James; Tremaine, Scott (2008). Galactic Dynamics (Second ed.). Princeton University Press. pp. 196-201. ISBN 978-0-691-13027-9 
  5. ^ a b Yoshida, H. (1992). “Symplectic Integrators for Hamiltonian Systems: Basic Theory”. IAU Symposium 152: 407. Bibcode: 1992IAUS..152..407Y. 
  6. ^ a b 陰山聡. “数値積分法”. 2020年7月6日閲覧。
  7. ^ “2008年度・数理解析・計算機数学2・第12回”. 2020年7月6日閲覧。
  8. ^ 宮武勇登. “保存則に即した数値計算手法”. 2020年7月7日閲覧。
  9. ^ Forest, Ronald; Ruth, D. (1990). “Fourth-order symplectic integration”. Physica D 43 (1): 105-117. doi:10.1016/0167-2789(90)90019-L. 
  10. ^ Suzuki, M. (1990). “Fractal decomposition of exponential operators with applications to many-body theories and Monte Carlo simulations”. Physics Letters A 146 (6): 319-323. doi:10.1016/0375-9601(90)90962-N. 
  11. ^ Yoshida, H. (1990). “Construction of higher order symplectic integrators”. Phys. Lett. A 150 (5–7): 262–268. Bibcode: 1990PhLA..150..262Y. doi:10.1016/0375-9601(90)90092-3. 
  12. ^ Wisdom, Jack; Holman, Mathew (1991). “Symplectic maps for the N-body problem”. Astronomical Journal 102: 1528-1538. doi:10.1086/115978. 
  13. ^ Wisdom, Jack; Holman, Mathew (1992). “Symplectic Maps for the n-Body Problem: Stability Analysis”. Astronomical Journal 104: 2022. doi:10.1086/116378. 
  14. ^ Kinoshita, Hiroshi; Yoshida, Haruo; Nakai, Hiroshi (1991). “Symplectic integrators and their application to dynamical astronomy”. Celestial Mechanics and Dynamical Astronomy 50 (1): 59-71. Bibcode: 1991CeMDA..50...59K. 

関連文献

  • Leimkuhler, Ben; Reich, Sebastian (2005). Simulating Hamiltonian Dynamics. Cambridge University Press. ISBN 0-521-77290-7.
  • Hairer, Ernst; Lubich, Christian; Wanner, Gerhard (2006). Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations (2 ed.). Springer. ISBN 978-3-540-30663-4.

関連項目