積分法SI01aがなぜ正準変換になっているのか、その由来を調べてみよう。
8.4.1 正準方程式の形式的厳密解
以下では相空間中の位置(つまり系の状態)r r=
(q p )
(8.14) の関数に作用する演算子をAˆのようにハットを付けて書くことにする。
正準変数の関数f とgのポアッソン括弧{f, g}を、fにD(g)ˆ という演算子が 作用したものとみることにしよう。つまり
D(g)(ˆ ·) ={·, g} (8.15) である。また、
Dˆ2(g)(·) ={{·, g}, g} (8.16) のように演算子のべきを定義する。ハミルトニアンH(q, p)の系の正準方程式は
dr
dt = ˆD(H)r (8.17)
と書ける。つまりD(Hˆ )は系の状態変化(時間微分)の演算子である。
時刻tの状態r=r(t)からτだけ時間が経った状態をr(t+τ)とする。ここ ではτは有限で、無限小とは限らない。テーラー展開より
r(t+τ) =r+ τ 1!
d dtr+τ2
2!
d2 dt2r+τ3
3!
d3
dt3r+· · · (8.18)
=r+ τ 1!
D(Hˆ )r+τ2 2!
D(Hˆ )2r+τ3 3!
D(Hˆ )3r+· · · (8.19) ここで演算子Aˆの指数関数を
Exp( ˆA) =
∑∞ n=0
1 n!
Aˆn (8.20)
数値計算のための解析力学 148
8.4シンプレクティック性の由来 と定義すると
r(t+τ) = Exp(τD(Hˆ ))r, (8.21) Exp(τD(Hˆ )) =
∑∞ n=0
τn n!
D(Hˆ ) (8.22)
演算子Exp(τD(Hˆ ))はハミルトニアンHの系の状態を時間τだけ進める(時間
推進の)演算子である。式(8.21)が示すように、正準方程式を積分するというこ とは、時間推進演算子Exp(τD(H))ˆ を求めることに他ならない。
前章で示したように、ハミルトニアンH(r)によって定まるr(t)⇒r(t+τ) という写像(つまり運動)はシンプレクティック変換なので、時間推進の演算子 Exp(τD(Hˆ ))はシンプレクティック変換である。
ある座標変換r(0)⇒r′=r(τ)について、
r′ = Exp(τD(h))ˆ r (8.23)
という形で書けるような関数h(q, p)が存在すれば、この変換は正準変換であり、
その変換はhをハミルトニアンとする系の(初期条件rから出発してr′に至る)
運動である。
8.4.2 時間推進演算子が解ける例
Hがpだけの関数、つまり
H =K(p) (8.24)
という形で与えられる場合には時間推進演算子を具体的に書くことが可能である。
これはポテンシャルの存在しない系の運動(自由運動)に相当する。このとき、
式(8.15)より
D(Hˆ ) = ˆD(K) = dK dp
∂
∂q (8.25)
である。一方、正準方程式に戻って考えると dr
dt = (dq
dt dp dt
)
= (dK
dp(p) 0
)
(8.26) である。この微分方程式は簡単に積分できて
r(t+τ) =
(q(t) +τdKdp(p(t)) p(t)
)
(8.27) である。従って時間推進演算子は
Exp(τD(K))ˆ r= (
q+τdK(p)dp p
)
(8.28) である。
数値計算のための解析力学 149
第8 章 シンプレクティック積分法 同様にハミルトニアンがqのみの関数
H=U(q) (8.29)
つまり
D(Hˆ ) = ˆD(U) =−dU dq
∂
∂p (8.30)
の場合についても時間推進演算子は厳密に計算できる。このとき正準方程式 dr
dt = (dq
dt dp dt
)
= ( 0
−dU(q)dq )
(8.31) の厳密解は
r(t+τ) =
( q(t) p(t)−τdUdq(q(t))
)
(8.32) 従って時間推進演算子は厳密に
Exp(τD(Uˆ ))r=
( q p−τdUdq(q)
)
(8.33) である。
8.4.3 合成変換による厳密解の近似
時間推進演算子Exp(τD(H))ˆ は、H がpだけの関数K(p)の時には式(8.28) として、Hがqだけの関数U(q)の時には式(8.33)として明示的に得られた。し かし、Hが
H(q, p) =K(p) +U(q) (8.34)
という形の時の時間推進演算子 Exp(τD(K) +ˆ τK(U)) =ˆ
∑∞ n=0
1 n!
(
τD(K) +ˆ τK(U)ˆ )n
= 1 +τD(K) +ˆ τK(Uˆ ) +τ2
2!
(D(K)ˆ 2+ ˆD(K) ˆD(U) + ˆD(U) ˆD(K) + ˆD(U)2 ) +τ3
3!
(D(K)ˆ 3+ ˆD(K)2D(Uˆ ) + ˆD(K) ˆD(U) ˆD(K) +· · ·)
+· · · (8.35)
を簡単に書くことはできない。
シンプレクティック積分法 SI01a は厳密な時間推進演算子Exp(τD(K) +ˆ τK(Uˆ ))を、(演算の形が既知である)二つの演算子Exp(τD(K))ˆ とExp(τD(Uˆ )) を組み合わせて
Exp(τD(Hˆ )) = Exp(τD(K) +ˆ τK(Uˆ )) (8.36)
∼Exp(τD(K)) Exp(τˆ D(Uˆ )) (8.37) 数値計算のための解析力学 150
8.4シンプレクティック性の由来 と近似する方法である。
r= (q
p )
= (qn
pn )
, r′= (q′
p′ )
= (qn+1
pn+1 )
, τ = ∆t (8.38)
とすると、シンプレクティック積分法SI01a[式(8.9)と(8.10)]は q′=q+τdK
dp (p) (8.39)
p′=p−τdU
dq(q′) (8.40)
である。これが ( q′
p′ )
= Exp(τD(K)) Exp(τˆ D(Uˆ )) ( q
p )
(8.41) であることを確認しよう。まず、
Exp(τD(K)) Exp(τˆ D(Uˆ ))q=Exp(τD(K))ˆ (
1−τdU dq
∂
∂p+· · · )
q
=Exp(τD(K))ˆ q
= (
1 +τdK dp
∂
∂q +τ2 (dK
dp )2
∂2
∂q2 +· · · )
q
=q+τdK dp
=q′ (8.42)
つぎに
Exp(τD(K)) Exp(τˆ D(Uˆ ))p=Exp(τD(K))ˆ (
1−τdU dq
∂
∂p+τ2 (dU
dq )2
∂2
∂p2 +· · · )
p
=Exp(τD(K))ˆ (
p−τdU dq
)
= (
1 +τdK dp
∂
∂q +τ2 (dK
dp )2
∂2
∂q2 +· · · ) (
p−τdU dq
)
=p−τ (
1 +τdK dp
∂
∂q +τ2 (dK
dp )2
∂2
∂q2 +· · · )
dU dq(q)
=p−τdU
dq(q+τdK
dp) [テーラー展開]
=p−τdU dq(q′)
=p′ (8.43)
数値計算のための解析力学 151
第8 章 シンプレクティック積分法 こうして式(8.41)が確認できた。
変換Exp(τD(K))ˆ と変換Exp(τD(Uˆ ))はどちらも正準変換であり、第6.4章で 示したように正準変換の合成変換は正準変換なので、Exp(τD(K)) Exp(τˆ D(Uˆ )) も正準変換である。こうして積分公式SI01aがシンプレクティック性を持つこと は確認するまでもなく自明となった。
近似式(8.37)の誤差を見積もってみよう。
Exp(τD(K)) Exp(τˆ D(Uˆ )) =
∑∞ n=0
τn n!
D(K)ˆ n
∑∞ n=0
τn n!
D(Uˆ )n
= 1 +τD(K) +ˆ τD(Uˆ ) +τ2
2!
(D(K)ˆ 2+ 2 ˆD(K) ˆD(U) + ˆD(U)2 )
+· · · (8.44)
なので、
Exp(τD(K) +ˆ τK(Uˆ )) = Exp(τD(K)) Exp(τˆ D(Uˆ )) + ˆErr(τ) (8.45) とすると、式(8.35)より
Err(τ) =ˆ τ2 2
(D(K) ˆˆ D(H)−D(Uˆ ) ˆD(K) )
+· · · (8.46) つまり1ステップの誤差はO(τ2)である。これは1次オイラー法と同じである。
1ステップあたりの精度は同じなのにSI01aにするとなぜ長時間積分しても全エ ネルギーがずれていかないのかその理由を次に探る。