シンプレクティック数値積分法 (シンプレクティックすうちせきぶんほう, 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}
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]
次式で定義される変換
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})}
[ 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次のシンプレクティック法は次式で与えられる[ 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次のシンプレクティック積分子は、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] [ 5]
なお、より高次のシンプレクティック積分子の系統的な構成方法は Suzuki (1990)[ 10] [ 11] ベーカー・キャンベル・ハウスドルフの公式 （英語版 ）
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] [ 13] [ 14] [ 1]
^ a b 伊藤孝士. “シンプレクティク数値積分法の天体力学的応用 ”. 2020年7月7日 閲覧。
^ a b c d 吉田春夫「可変時間ステップによるシンプレクティック数値解法(非線形可積分系による応用解析) 」『数理解析研究所講究録』第889巻、京都大学数理解析研究所、1994年11月、70-76頁、CRID 1050001202298760960 、hdl :2433/84362 ISSN 1880-2818 、2024年1月12日 閲覧
^ 牧野淳一郎, 福重俊幸, 小久保英一郎, 川井敦, 台坂博, 杉本大一郎 (2007年3月13日). “N体シミュレーション啓蟄の学校教科書 ”. 国立天文台. p. 51. 2020年5月24日 閲覧。 ^ a b Binney, James; Tremaine, Scott (2008). Galactic Dynamics (Second ed.). Princeton University Press. pp. 196-201. ISBN 978-0-691-13027-9
^ a b Yoshida, H. (1992). “Symplectic Integrators for Hamiltonian Systems: Basic Theory”. IAU Symposium 152 : 407. Bibcode : 1992IAUS..152..407Y .
^ a b 陰山聡. “数値積分法 ”. 2020年7月6日 閲覧。
^ “2008年度・数理解析・計算機数学2・第12回 ”. 2020年7月6日 閲覧。 ^ 宮武勇登. “保存則に即した数値計算手法 ”. 2020年7月7日 閲覧。 ^ Forest, Ronald; Ruth, D. (1990). “Fourth-order symplectic integration”. Physica D 43 (1): 105-117. doi :10.1016/0167-2789(90)90019-L . ^ 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 . ^ 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 . ^ Wisdom, Jack; Holman, Mathew (1991). “Symplectic maps for the N-body problem”. Astronomical Journal 102 : 1528-1538. doi :10.1086/115978 . ^ Wisdom, Jack; Holman, Mathew (1992). “Symplectic Maps for the n-Body Problem: Stability Analysis”. Astronomical Journal 104 : 2022. doi :10.1086/116378 . ^ 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 . 
![{\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)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/596b2d0242ddd33221ab035c2913b288e8fd6795)
![{\displaystyle x_{1}={\frac {1}{2-{\sqrt[{3}]{2}}}},\ \ x_{0}=1-2x_{1}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/0345c9f89971a9022edecd57f7e34bdfef5a15e4)
