シンプレクティック数値積分法
シンプレクティック数値積分法 (シンプレクティックすうちせきぶんほう, symplectic integrator) とは、正準力学系の運動方程式に特化した常微分方程式の数値解法のことをいう。系のシンプレクティック形式およびハミルトニアンを保存するため、ルンゲ=クッタ法のような汎用の数値積分法に比べて良い性質を示す。このために天体力学などの分野で採用されている。
正準力学系において、(q_i, p_i) を正準変数、H = H ( q_i, p_i ) をハミルトニアンとするとき、その運動方程式はハミルトンの正準方程式
\frac{ d q_i }{ d t } = \frac{ \partial H }{ \partial p_i } ( q, p ) , \ \ \frac{ d p_i }{ d t } = - \frac{ \partial H }{ \partial q_i } ( q, p )
である。これらの運動方程式の解 \left( q_i ( t ), p_i ( t ) \right) は一般に次の性質を満足する。
シンプレクティック形式 \omega = \sum_i dp_i \wedge dq_i を不変に保つ。
ハミルトニアン H \left( q_i ( t ) , p_i ( t ) \right) を不変に保つ。
ところが、正準方程式を数値的に解くためにルンゲ=クッタ法のような汎用の数値積分アルゴリズムを適用すると、一般に数値解においてこれらの性質が破れ、長時間の積分によりエネルギーが保存しないなどの非物理的な結果を生じ得る。シンプレクティック積分法は厳密にシンプレクティック写像であるような(すなわちシンプレクティック形式を保存する)数値積分アルゴリズムであり、ハミルトン力学系の数値解析手法としてより優れた性質を示す。例えば調和振動子
H = \frac{ 1 }{ 2 } p^2 + \frac{ 1 }{ 2 } q^2
に2次のシンプレクティック積分法を適用すると、時間刻み幅を h として、真のハミルトニアンの代わりに
\tilde{H} = \frac{ 1 }{ 2 } p^2 + \frac{ 1 }{ 2 } q^2 - \frac{ h^2 }{ 4 } q^2
を保存するため、数値解はある楕円軌道の上に留まり、エネルギーの単調な増加または減少を生じない。さらに、偶数次のシンプレクティック積分子は時間反転対称性を持つという利点も存在する。
一方で、時間刻み幅を動的に変更する適応時間刻みを単純にシンプレクティック積分子に適用するとハミルトニアンの保存が破れることが知られている。
ハミルトニアン H が2つの可積分ハミルトニアン H_A, H_B の和であると仮定する。
H ( q, p ) = H_A ( q, p ) + H_B ( q, p )
例えばポテンシャル V 中の質量 m の粒子という系の場合 H_A = \frac{ 1 }{ 2 m } p^2, H_B = V ( q ) であり、この仮定を満足する。H, H_A, H_B に対応するハミルトンベクトル場をそれぞれ D, A, Bと書くとき
D = A + B
が成立し、それぞれのハミルトンベクトル場に沿う時間 t の発展すなわち指数写像を S ( t ) := \exp \left( t D \right), \exp ( t A ), \exp ( t B ) と書く(これはシンプレクティック写像である)。H_A および H_B が可積分であるという仮定により \exp ( t A ), \exp ( t B ) はその具体的な表示が既知である。ここでの問題は、真のハミルトニアンに関する指数写像
S ( t ) = \exp \left( t D \right) = \exp \left\{ t ( A + B ) \right\}
を N 区間の時間積分の集積 S ( t ) = \prod_{i = 1}^N S ( h ) (ここに h = t / N とおいた) と書くとき、S ( h ) を既知のシンプレクティック写像 \exp ( h A ), \exp ( h B ) を用いて
\exp \left\{ h ( A + B ) \right\} = \prod_k \left[ \exp ( c_k h A ) \, \exp ( d_k h B ) \right] + \mathcal{O} \left( h^{n+1} \right)
という形に近似することである。この右辺が求める n 次のシンプレクティック積分子であり、シンプレクティック性を満足することが保証され、ハミルトニアン
\tilde{H} = H + h^n H_n + \mathcal{O} \left\{ h^{n+1} \right\}
を保存する。
次式で定義される変換 S_\mathrm{1st} ( h ) は S ( h ) = S_\mathrm{1st} ( h ) + \mathcal{O} ( h^2 ) を満足する1次のシンプレクティック積分子である。
S_\mathrm{1st} ( h ) = \exp ( h A ) \, \exp ( h B )
特に H_A = \frac{ 1 }{ 2 } p^2, H_B = V ( q ) の場合、変換 \exp ( h B ) は ( q, p ) \mapsto ( q, p - h V'(q) )、変換 \exp ( h A ) は ( q, p )\mapsto ( q + h p, p ) と表示できるため、このスキーム S_\mathrm{1st} = \exp ( h A ) \exp ( h B ) 全体としては
p_{n+1} = p_n - h \frac{ d V }{ d q } ( q_n )
q_{n+1} = q_n + h p_{n+1}
と表示できる。このスキームはオイラー法を修正したものとみなせるため、シンプレクティックオイラー法と呼ばれることもある。
2次のシンプレクティック法は次式で与えられる。なおこの積分スキームはリープ・フロッグ法あるいはベレの方法、Strörmer法など分野毎に異なった名称で知られている。
S_\mathrm{2nd} ( h ) = \exp \left( \frac{ h }{ 2 } A \right) \, \exp ( h B ) \, \exp \left( \frac{ h }{ 2 } A \right)
上述の H_A = \frac{ 1 }{ 2 } p^2, H_B = V ( q ) の場合にはこれは次のスキームに帰着する。
q_{n+1/2} = q_n + h \frac{ 1 }{ 2 } p_n
p_{n+1} = p_n - h \frac{ d V }{ d q } \left( q_{n+1/2} \right)
q_{n+1} = q_{n+1/2} + h \frac{ 1 }{ 2 } p_{n+1}
4次のシンプレクティック積分子は、2次の積分子を異なる時間ステップで三度適用することにより得られる。
S_\mathrm{4th} ( h ) = S_\mathrm{2nd} ( x_1 h ) \, S_\mathrm{2nd} ( x_0 h ) \, S_\mathrm{2nd} ( x_1 h )
x_1 = \frac{ 1 }{ 2 - \sqrt[3]{2} } , \ \ x_0 = 1 - 2 x_1
これは Forest & Ruth (1990) によって導かれた 後、Yoshida (1990) によって2次のシンプレクティック積分を三度適用したものに等しいことが指摘された。
なお、より高次のシンプレクティック積分子の系統的な構成方法は Suzuki (1990) および Yoshida (1990)によって与えられている。Yoshida (1990) はベーカー・キャンベル・ハウスドルフの公式を適用することにより高次の 2 n 次シンプレクティック積分子を解析的に構成しているが, 次数が増大すると必要なステップ数 (S_\mathrm{2nd} の数) が指数関数的に増大し効率が悪化することも指摘し、より効率的なシンプレクティック積分子をいくつか数値的に求めてもいる。
H = \frac{ 1 }{ 2 m } \mathbf{p}^2 - \frac{ G M m }{ \left| \mathbf{q} \right| }
にシンプレクティック積分法を適用すると、オイラー法やルンゲ=クッタ法では時間の経過とともに数値誤差が累積しエネルギーの誤差が増大するが、シンプレクティック積分法ではエネルギー誤差の累積は見られない(右図)。
中心天体からの重力が支配的であるような天体力学の典型的な問題では、系のハミルトニアンは中心天体によるケプラー問題の部分 H_\mathrm{Kep} と、天体間相互作用による摂動部分 H_\mathrm{int} に分割できる。
H = H_\mathrm{Kep} + H_\mathrm{int}
それを踏まえて、ハミルトニアンをこの二項に分割してシンプレクティック積分子を適用する数値解析手法が Wisdom & Holman (1991, 1992) および Kinoshita et al. (1991) によって提案され、従来の手法より誤差が小さくなるなどの良好な性質を持つことが示された。