−pk+1+pk−∂H
∂q (qk+1, pk) ∆t= 0 (7.47) まとめると
qk+1=qk+ ∆t∂H
∂p(qk+1, pk) (7.48)
pk+1=pk−∆t∂H
∂q (qk+1, pk) (7.49)
を得る。
ハミルトニアンが
H(q, p) =K(p) +U(q) (7.50)
の形になっているときには式(7.48)と(7.49)は qk+1=qk+ ∆t∂K
∂p(pk) (7.51)
pk+1=pk−∆t∂U
∂q(qk+1) (7.52)
となり、陽的なアルゴリズムになる。これを（一次陽的）シンプレクティック積 分法とよぶ。
7.5 オイラー = ラグランジュ方程式
xの関数y(x)が、固定された端点
ya=y(a), yb=y(b) (7.53) という条件の下で、汎関数
I(y, y′, x) =
∫ b a
L(y(x), y′(x), x)dx (7.54) がy=y(x)で極値をとる、つまり
δI= 0 (7.55)
が成り立つ条件を求めておこう。
y(x)を
y(x)→y(x) +ϵ η(x) (7.56)
と変えたときの汎関数Iは I(y+ϵη, y′+ϵη′, x) =
∫ b a
L(y(x) +ϵη(x), y′(x) +ϵη′(x), x)dx (7.57) 数値計算のための解析力学 137
第7 章 ハミルトンの原理
である。ここで、被積分関数をϵでテーラー展開すると L(y+ϵη, y′+ϵη′, x) =L(y, y′, x) +∂L
∂y ϵη(x) + ∂L
∂y′ϵη′(x) +O(ϵ2) (7.58) なので
I(y+ϵη, y′+ϵη′, x) =
∫ b a
L(y, y′, x)dx +ϵ
∫ b a
( η(x)∂L
∂y dx+η′(x)∂L
∂y′ )
dx+O(ϵ2) (7.59)
=I(y, y′) +ϵ δI(y, y′) +O(ϵ2) (7.60) つまり、
δI=δ
∫ b a
L(y, y′, x)dx (7.61)
=
∫ b a
( η(x)∂L
∂y +η′(x)∂L
∂y′ )
dx (7.62)
(7.63) である。ここで、式(7.53)より
η(a) =η(b) = 0 (7.64)
なので、
∫ b a
η′(x) (∂L
∂y′ )
dx= [
η(x) (∂L
∂y′ )]b
a
−
∫ b a
η(x) d dx
(∂L
∂y′ )
dx (7.65) の右辺第一項はゼロになる。すると
δI=δ
∫ b a
L(y, y′, x)dx=
∫ b a
{
−d dx
(∂L
∂y′ )
+∂L
∂y }
η(x)dx= 0 (7.66) となり、任意のη(x)に対して変分δI= 0なので
d dx
(∂L
∂y′ )
−∂L
∂y = 0 (7.67)
が成り立つ。これをオイラー=ラグランジュの方程式という。
関数がN個、yi(x) (i= 1, . . . , N)あるときも同様である。yi(a)とyi(b) (i= 1, . . . , N)が固定という端点条件の下、汎関数
I(y1, y2,· · ·, yN, y1′, y2′,· · · , y′N) =
∫ b a
L(y1, y2,· · ·, yN, y1′, y′2,· · · , y′N, x)dx 数値計算のための解析力学 138
7.5オイラー=ラグランジュ方程式 が極値をとる関数yi(x) (i= 1, . . . , N)は、N個のオイラー=ラグランジュの方 程式
d dx
(∂L
∂yi′ )
−∂L
∂yi
= 0 (i= 1,· · · , N) を満たす。
例として、p.129で紹介した極小曲面を求めてみよう。関数y =y(x)をx軸 の周りに回転してできる側面の面積をSとすると、
S(y, y′) =
∫ a 0
y(x)√
1 +y′(x)2dx (7.68) であった。
L(y, y′) =y√
1 +y′2 (7.69)
とすれば、極小曲面は、オイラー=ラグランジュの方程式(7.67)の解である。実 際に計算してみよう。
∂L
∂y =√
1 +y′2 (7.70)
∂L
∂y′ = yy′
√1 +y′2 (7.71)
ここまでは簡単。次が計算ミスしやすいので注意。
d dx
(∂L
∂y′ )
= d dx
( yy′
√1 +y′2 )
(7.72)
=y′2+yy′′
√1 +y′2 − · · · (7.73) 上の結果を、オイラー＝ラグランジュの方程式(7.67)に代入し整理すると、
yy′′= 1 +y′2 (7.74)
となる。なお、この微分方程式の解は双曲線関数y= cosh(x)である。
ハミルトンの原理
δS=δ
∫ tb ta
L(q,q), dt˙ = 0 (7.75) と端点条件
δq(ta) =δq(tb) = 0 (7.76)
から得られるオイラー=ラグランジュの方程式がラグランジュの運動方程式
∂L
∂q − d dt
(∂L
∂q˙ )
= 0 (7.77)
である。
数値計算のための解析力学 139
第7 章 ハミルトンの原理
最後に、相空間中の軌道q(t), p(t)が極値をとるという拡張ハミルトンの原理 δS∗=δ
∫ tb ta
(pq˙−H(q, p))dt= 0 (7.78) と端点条件
δq(ta) =δq(tb) = 0 (7.79)
から得られるオイラー=ラグランジュの方程式が正準方程式 dq
dt = ∂H
∂p (7.80)
dp
dt =−∂H
∂q (7.81)
になるが、この際、汎関数S∗にはp˙が現れないので、pに関する端点条件は不要 であることに注意しよう。
数値計算のための解析力学 140
Chapter 8
シンプレクティック積分法
8.1 シンプレクティック積分法
前の章の最後に、変分原理（修正ハミルトンの原理）からシンプレクティック 積分法という積分法が導出されることを見た。この章ではこのスキームについて 詳しく調べる。なお、この章では簡単のため、1自由度系について述べる。多自 由度系への一般化は容易である。また、ハミルトニアンが運動エネルギーKとポ テンシャルUの和
H(q, p) =K(p) +U(q) (8.1)
で書かれており、KとU がそれぞれpとqだけに依存するような場合について 主に考察する。
なお、この章ではいくつかのサンプルプログラムと共に、実際の数値計算例 を（エネルギーの時間変化のグラフなどと共に）示すが、そこではそれぞれ のスキームの特徴を示すために時間刻み幅∆tをあえて大きくとっているこ とに注意すること。
正準座標
r= (q, p) から、別の正準座標
R= (Q, P)
への変換が正準変換になっているとき、この変換はシンプレクティック条件を満 たすことは以前示した。 ある変換が正準変換になっていることを、「変換がシン プレクティックである」、あるいは「変換がシンプレクティック性をもつ」と表 現する。シンプレクティック性をもつ数値積分法をシンプレクティック積分法と いう。
例題として、以前も考えたリング上を質点が滑るバネ=質点系をもう一度とり あげよう。x軸とy軸に接する半径1の円周上を質量mの質点が滑らかに滑る。
原点と質点はバネ（ばね定数k、自然長0）で結ばれている。図のようにx軸と 141
第8 章 シンプレクティック積分法
1 1
q
なす角度qとする。ハミルトニアンを再掲すると H(q, p) = p2
2m+k(cosq+ sinq) (8.2)
このハミルトニアンは式(8.1)の形になっていて K(p) = p2
2m (8.3)
U(q) =k(cosq+ sinq) (8.4)
である。
8.2 1 次陽的シンプレクティック法
前章で確認したように、1次（陽的）オイラー法はシンプレクティックになっ ていない。そのためにどのような悪影響があるのか、実際にリング上のバネ質点
系(8.2)を、1次オイラー法で解いてみよう。正準方程式を再掲すると、
dq dt = p
m (8.5)
dp
dt =k(sinq−cosq) (8.6)
プログラムの中心部分は以下の通りである：
Listing 8.1: one ball on ring 1stEuler.cpp 1
2 v o i d e u l e r 1 s t ( s t r u c t p a r t i c l e ∗p a r t i c l e , d o u b l e dt ) 3 {
4 // H a m i l t o n i a n
数値計算のための解析力学 142
8.21次陽的シンプレクティック法 5 // H( q , p ) = p ˆ 2 / ( 2m) + k ( c o s q + s i n q )
6 //
7 d o u b l e q = p a r t i c l e−>pos [ 0 ] ; 8 d o u b l e p = p a r t i c l e−>pos [ 1 ] ; 9 d o u b l e dq , dp ;
10 dq = ( p ∗ MASS INV ) ∗ dt ;
11 dp = ( SPRING K ∗ ( s i n ( q ) − c o s ( q ) ) ) ∗ dt ;
12 q += dq ;
13 p += dp ;
14
15 p a r t i c l e−>pos [ 0 ] = q ; 16 p a r t i c l e−>pos [ 1 ] = p ; 17
18 s t d : : c o u t << ” t o t a l e n e r g y : ” 19 << s t d : : s c i e n t i f i c
20 << t o t a l e n e r g y ( p a r t i c l e−>pos ) 21 << s t d : : e n d l ;
22 }
たとえば、(q, p) = (2,0)の初期条件の下で、これを実行すると・・・明らか に変なことが起きている。はじめはq=−3π/4を中心に振動していた質点が、次 第にその振幅を大きくしていき、ついには回転運動をはじめてしまう。これは系 の全エネルギーが増加してしまっていることを意味する。実際に全エネルギーを プロットしてみるとそれが確認できる。
1次オイラー法は、1ステップでO(∆t2)の誤差がある。有限時間T まで積分 するにはN/∆tステップだけ1次オイラー法を繰り返し適用するので、最終的に
はO(∆t)の誤差が生じるのは仕方がない。しかしながらこの例のように全エネル
ギーが変わると解の振る舞いが（振動状態と回転状態の違いのように）定性的に 異なる系に対しては、どれだけ長時間（多ステップ）積分しても全エネルギーは 数値計算のための解析力学 143
第8 章 シンプレクティック積分法 なるべく変わらないようにしたい。
1ステップでO(∆t)の精度をもつ数値積分法（＝1次精度積分法）は、（O(∆t2) だけの自由度があるので）無限のバリエーションが存在する。1次オイラー法は その中の一つにすぎない。この章では無数にある1次精度積分法の中から「シン プレクティック性をもつ」という条件を付けるとエネルギーの保存性がどうなる か調べよう。なお後で見るように1次精度のシンプレクティックな数値積分法も 複数存在する。
ハミルトニアンが式(8.1)で与えられるとき、1次オイラー法は、
qn+1=qn+ ∆tdK
dp(pn) (8.7)
pn+1=pn−∆tdU
dq(qn) (8.8)
と書ける。これをほんの少し（赤字の部分）だけ変えて 1次シンプレクティック積分法: SI01a
qn+1=qn+ ∆tdK
dp(pn) (8.9)
pn+1=pn−∆tdU
dq(qn+1) (8.10)
とする。実はこの積分法はシンプレクティックな積分法になっている。
r= (qn, pn)からR= (qn+1, pn+1)への変換が正準変換であることをポアッ ソン括弧で調べてみよう。1自由度系の場合には{Q, Q}={P, P}= 0は自明な ので、{Q, P}だけを調べれば良い。
Q=q+ ∆tdK(p) dp P =p−∆tdU(Q)
dq =p−∆tdU
dq (q+ ∆tdK(p) dp ) に注意すると、
{Q, P}= ∂Q
∂q
∂P
∂p −∂P
∂q
∂Q
∂p
= 1× [
1− (
∆td2U dq2
) (
∆td2K dp2
)]
− (
−∆td2U dq2
) (
∆td2K dp2
)
= 1−∆t2d2U dq2
d2K
dp2 + ∆t2d2U dq2
d2K dp2
= 1 (8.11)
従って式(8.9)と(8.10)は正準変換である。この積分法は1次精度の陽的なシン
プレクティック法である。ここではこのスキームをSI01aと呼ぶことにする。
数値計算のための解析力学 144
8.21次陽的シンプレクティック法
式(8.9)と(8.10)に基づいてコードを作ってみよう。見ての通り、このコード
は1次オイラー法とほとんど同じである。
Listing 8.2: one ball on ring 1stSymplectic.cpp 1
2 v o i d s y m p l e c t i c 1 s t ( s t r u c t p a r t i c l e ∗p a r t i c l e , d o u b l e dt ) 3 {
4 // H a m i l t o n i a n
5 // H( q , p ) = p ˆ 2 / ( 2m) + k ( c o s q + s i n q )
6 //
7 d o u b l e p = p a r t i c l e−>pos [ 1 ] ; 8 d o u b l e q = p a r t i c l e−>pos [ 0 ] ; 9 d o u b l e dq , dp ;
10
11 dq = ( p ∗ MASS INV ) ∗ dt ;
12 q += dq ;
13
14 dp = ( SPRING K ∗ ( s i n ( q ) − c o s ( q ) ) ) ∗ dt ;
15 p += dp ;
16
17 p a r t i c l e−>pos [ 0 ] = q ; 18 p a r t i c l e−>pos [ 1 ] = p ; 19
20 s t d : : c o u t << ” t o t a l e n e r g y : ” 21 << s t d : : s c i e n t i f i c
22 << t o t a l e n e r g y ( p a r t i c l e−>pos ) 23 << s t d : : e n d l ;
24 }
このコードを実行すると、図のようにエネルギーは（振動こそすれ）オイラー 法の時のように増大し続けることはなく、質点の運動が振動運動が回転運動に変 わってしまうような明らかに定性的におかしな解になることもない。
1次シンプレクティック法SI01a［式(8.9)と(8.10)］と1次オイラー法［(8.7)
と(8.8)］との違いはほんのわずかである。式(8.10)で前の時間ステップの値qn
を使うと1次オイラー法となる。このわずかな違いによって積分がシンプレク ティック（正準変換）になること、そしてその効果が大きいことは実に印象的で ある。
1次シンプレクティック法SI01a［式(8.9)と(8.10)］では、式(8.9)でまずqn の値をqn+1に更新してからその新しいqn+1を使って式(8.10)をつかってpnを 更新する。この順番を逆にした以下のスキーム（SI01b）もシンプレクティック 性を持つことを確認することができる。
数値計算のための解析力学 145
第8 章 シンプレクティック積分法
1次シンプレクティック積分法: SI01b
pn+1=pn−∆tdU
dq(qn) (8.12)
qn+1=qn+ ∆tdK
dp(pn+1) (8.13)