シンプレクティック積分法についてのメモ - 神戸大学

シンプレクティック積分法について∗
陰山聡
神戸大学システム情報学研究科計算科学専攻
ver. 140110a
1
正準方程式の厳密解
1 自由度系の相空間 (q, p) 中の点を
(
r=
q
p
)
(1)
ˆ などと書くことにする。
と書く。関数 f (r) に作用する演算子をハットをつけて A
ˆ
関数 h(r) をパラメータとし、f (r) に作用する演算子 D(h)
をポアッソン括弧を使い
∂f ∂h ∂h ∂f
ˆ
D(h)f
= {f, h} =
−
∂q ∂p
∂q ∂p
(2)
と定義すると、ハミルトニアンが H で与えられる系の正準方程式は
dr
ˆ
= D(H)
r
dt
(3)
である。
時刻 t の状態を r = r(t)、そこから τ だけ時間が経った状態を r(t + τ ) とする。テーラー展開と式 (3) より
)
(∞
)
(∞
∑ τ n dn
∑ τn
ˆ
D(H)
r
r(t + τ ) =
r=
n! dtn
n!
n=0
ˆ の指数関数を Exp(A)
ˆ =
演算子 A
∑∞
n=0
n=0
Aˆn /n! と定義すると、正準方程式の形式的な厳密解は
ˆ
r(t + τ ) = Exp(τ D(H))
r(t)
(4)
ˆ
と書ける。つまり Exp(τ D(H))
は状態を時間 τ だけ進める時間推進の演算子である。
ˆ
一般のハミルトニアン H(r) に対する時間推進演算子 Exp(τ D(H))
を具体的に計算することができないが、
n
ˆ
H が p だけに依存する関数 K(p) の時は D(K)
r = 0 (n ≥ 2) なので
(
)
(
)
(
)
q + τ dK(p)
dK ∂
dp
ˆ
ˆ
Exp(τ D(K)r = 1 + τ D(K) r = 1 + τ
r=
(5)
dp ∂q
p
と厳密に書ける。任意の τ について成り立つこの式は等速直線運動の演算子である。同様に H = U (q) の時に
も時間推進演算子の厳密な形を具体的に求めることができ、
(
)
)
(
(
)
q
∂
dU
ˆ )r = 1 + τ D(U
ˆ ) r = 1−τ
r=
(6)
Exp(τ D(U
dq ∂p
p − τ dU (q)
dq
∗ このノートは神戸大学工学部の講義「解析力学
B」(2013 年度後期）の講義資料から一部を抜き出したものである。
1
である。
ˆ
ˆ )r のような厳密に解ける
以下に見るように、シンプレクティック積分法は上の Exp(τ D(K)r
や Exp(τ D(U
ˆ
時間推進演算子を組み合わせて（合成変換して）本来の時間推進演算子 Exp(τ D(H)r
を近似する手法である。
ˆ )) はシ
なお、任意の運動は正準変換なので、式 (4) から、ある関数 f (r) から構成される演算子 Exp(τ D(f
ンプレクティック性をもつ演算子であることがわかる。また、r(t) から r ′ = r(t + τ ) への変換が式 (4) の形で
書けるのであればこれは（H をハミルトニアンとする）運動と解釈することができる。
2
シンプレクティック積分法
ハミルトニアンが
H(q, p) = K(p) + U (q)
(7)
のときの 1 次精度シンプレクティック積分法の一つとして以下の公式がある。
dK
(p)
dp
dU ′
p′ = p − τ
(q )
dq
q′ = q + τ
(8)
(9)
ここで τ < 1 である。これは
dK
(p)
dp
dK
p∗ = p − τ
(q)
dq
q∗ = q + τ
dU ∗
(p )
dp
dU ∗
p′ = p∗ − τ
(q )
dq
q′ = q∗ + τ
(10)
(11)
(12)
(13)
という二つの変換に分けられる。式 (10) と (11) は、ハミルトニアン K による厳密な時間推進 (5) であり、式
(12) と (13) は、ハミルトニアン U による（新しい正準座標 r ∗ の下で計算した）厳密な時間推進 (6) である。
つまりこの積分法は演算子を使って
r ′ = Sˆ r
(14)
ˆ
ˆ
ˆ
S = Exp(τ D(K)) Exp(τ D(U ))
(15)
と書ける。先に作用させる K の演算子が左にくることに注意する1 。
式 (15) は次のようにしても確認できる。式 (9) を書き換えると
dU
dK
(p + τ
)
dq
dp
(
)
(
)2 3
(
)3 4
dU
dK d2 U
τ 2 dK
d U
τ 3 dK
d U
=p−τ
+τ
+
+
+ ···
dq
dp dq 2
2! dp
dq 3
3! dp
dq 4
}
{∞
∑ τ n ( dK d )n dU
=p−τ
n! dp dq
dq
n=0
{∞
}(
)
∑ τn
dU
n
ˆ
=
D(K)
p−τ
n!
dq
n=0
(
)
ˆ
ˆ ) p
= Exp(τ D(K))
1 + τ D(U
p′ = p − τ
ˆ
ˆ )) p
= Exp(τ D(K))
Exp(τ D(U
(16)
1 たとえば、3 次元カーテシアン座標 (x, y, z) の z 軸周りの回転 R
ˆ と x 軸方向の平行移動 Tˆ は非可換である。局所座標を回転（R
ˆ）
ˆ Tˆ r である。
してから変換後の x 軸方向に平行移動（Tˆ ）するという操作は r ′ = R
2
また、式 (8) も同じ形に書けることがわかるので、式 (15) が成り立つ。
このシンプレクティック積分法はハミルトニアン H による（未知の）時間推進演算子を、二つの（既知の）
時間推進演算子の合成で
ˆ
ˆ
ˆ
ˆ )) = Sˆ
Exp(τ D(H))
= Exp(D(K
+ U )) ∼ Exp(τ D(K))
Exp(τ D(U
(17)
と近似した方法と言える。
ˆ
ˆ )) はどちらも正準変換であり、正準変換の合成変換は正準変換なので、
変換 Exp(τ D(K))
と変換 Exp(τ D(U
ˆ
この演算子 S はシンプレクティック性を持つことがわかる。
3
計算精度
ˆ
ˆ とすると、二つの演算子の展開
近似式 (17) の誤差を見積もってみよう。Exp(τ D(H))
= Sˆ + Err
)2
(
)
2 (
ˆ
ˆ ) + ···
ˆ
ˆ
ˆ ) + τ D(K)
+ τ D(U
Exp(τ D(H))
= 1 + τ D(K)
+ τ D(U
2!
(
)(
)
2
τ ˆ
τ2 ˆ
2
2
ˆ
ˆ
ˆ
S = 1 + τ D(K) + D(K) + · · ·
1 + τ D(U ) + D(U ) + · · ·
2!
2!
を比較して
)
τ2 ( ˆ
ˆ ) − D(U
ˆ )D(K)
ˆ
ˆ
D(K)D(U
+ ···
Err(τ
)=
2
(18)
つまり 1 ステップの誤差は O(τ 2 ) である。これは 1 次オイラー法と同じである。
4
エネルギーの計算精度
1 ステップあたりの精度は同じなのに、1 次オイラー法とは異なり、シンプレクティック法は長時間積分し
ても全エネルギーが本来の値から大きくずれることはない。その理由は以下の通りである。
シンプレクティック積分法の演算子 Sˆ は正準変換なので、それに対応する（仮想的な）運動があると期待で
きる。実際
ˆ
ˆ )) = Exp(τ D(
ˆ H))
˜
Sˆ = Exp(τ D(K))
Exp(τ D(U
(19)
と書けることが以下のように確認できる。
ˆ B
ˆ に対して以下の公式（Baker-Campbell-Hausdorﬀ 公式）が成り立つ。
一般に非可換演算子 A,
ˆ Exp(B)
ˆ = Exp(C)
ˆ
Exp(A)
(20)
ˆ + 1 [A,
ˆ B]
ˆ + 1 [Aˆ − B,
ˆ [A,
ˆ B]]
ˆ − 1 [B,
ˆ [A,
ˆ [A,
ˆ B]]]
ˆ + ···
Cˆ = Aˆ + B
2
12
24
(21)
ˆ B]
ˆ = AˆB
ˆ −B
ˆ Aˆ
[A,
(22)
ˆ ), D(g)]
ˆ
ˆ
[D(f
= D({g,
f })
(23)
{f, {g, h}} + {g, {h, f }} + {h, {f, g}} = 0
(24)
ˆ
ˆ )) = Exp(C)
ˆ
Sˆ = Exp(τ D(K))
Exp(τ D(U
(25)
α = {U, K}
(26)
ここで
等である。定義から、任意の f , g に対して
が成り立つ。これはヤコビ恒等式
を使えば証明できる。
3
とすると
2
τ3 ˆ
τ4 ˆ
ˆ
ˆ ) + τ D(α)
ˆ
ˆ
ˆ
ˆ
Cˆ = τ D(K)
+ τ D(U
+ [D(K
− U ), D(α)]
− [D(U
), [D(K),
D(α)]]
+ ···
2
12
24
2
τ3 ˆ
τ4 ˆ
ˆ
ˆ ) + τ D(α)
ˆ
+ D({α,
K − U }) − D({{α,
K} , U }) + · · ·
= τ D(K)
+ τ D(U
2
12
24
τ
τ2
τ3
ˆ
= τ D(K
+U + α+
{α, K − U } −
{{α, K} , U } + · · · )
2
12
24
(27)
従って
2
3
˜ = K + U + τ α + τ {α, K − U } − τ {{α, K} , U } + · · ·
H
(28)
2
12
24
ˆ H))
˜ による推進を厳密に実装したものなので、どれだけ
を得る。シンプレクティック積分法は演算子 Exp(τ D(
˜ を一定に保つ。上の式は、その H
˜ が本来の系のエネルギー H とは O(τ ) だけしかずれて
積分しても、常に H
いないことを意味する。つまりシンプレクティック積分法ではエネルギーの誤差は蓄積しない。
5
その他
式 (15) において K による推進と H による推進の順番を入れ替えて
ˆ
ˆ )) Exp(τ D(K))
ˆ
Exp(τ D(H))
∼ Exp(τ D(U
(29)
とすれば別の 1 次精度のシンプレクティック法ができる。具体的なアルゴリズムは
dU
(q)
dq
dK ′
(p )
q′ = q + τ
dp
p′ = p − τ
(30)
(31)
である。
これらのシンプレクティック法は 1 次精度である。精度を上げるためには次のようにすればよい。一般に演
ˆとB
ˆ に対して
算子 A
ˆ Exp(τ B)
ˆ Exp( τ A)
ˆ + O(τ 3 )
ˆ = Exp( τ A)
(32)
Exp(τ Aˆ + τ B)
2
2
なので、そこで、ハミルトニアン H = K(p) + U (q) に対する時間推進演算子を
τ ˆ
ˆ
ˆ )) Exp( τ D(K))
ˆ
Exp(τ D(H))
∼ Exp( D(K))
Exp(τ D(U
2
2
(33)
と近似すれば、2 次精度の数値積分法になっている。右辺は正準変換の合成変換なので、これは明らかにシン
プレクティックな積分法である。より高次のシンプレクティック法を構成方法については [2] を参照のこと。
References
[1] 吉田春男. ハミルトニアン力学系のためのシンプレクティック数値積分法. 共同研究「非線形現象の数理科
学」湘南レクチャー論文集, p. 68–83, 1997.
[2] H. Yoshida. Construction of Higher Order Symplectic Integrators. Phys. Lett., Vol. 102, pp. 262–268.,
1990
4

Download Report