在常微分方程 的數值計算中,幾何積分 是一種保留微分方程的流的精確幾何特性的數值方法。
可考慮單擺 運動以引出幾何積分的研究。
設擺錘質量為
m
=
1
{\displaystyle m=1}
,擺杆長度為
ℓ
=
1
{\displaystyle \ell =1}
。設重力加速度為
g
=
1
{\displaystyle g=1}
。用
q
(
t
)
{\displaystyle q(t)}
表示杆偏移垂直方向的角位移,並用
p
(
t
)
{\displaystyle p(t)}
表示擺的動量,則系統的哈密頓量 (動能 與勢能 之和)為
H
(
q
,
p
)
=
T
(
p
)
+
U
(
q
)
=
1
2
p
2
−
cos
q
,
{\displaystyle H(q,p)=T(p)+U(q)={\frac {1}{2}}p^{2}-\cos q,}
其給出哈密頓方程
(
q
˙
,
p
˙
)
=
(
∂
H
/
∂
p
,
−
∂
H
/
∂
q
)
=
(
p
,
−
sin
q
)
.
{\displaystyle ({\dot {q}},{\dot {p}})=(\partial H/\partial p,-\partial H/\partial q)=(p,-\sin q).\,}
很自然,可將所有
q
{\displaystyle q}
的位形空間
Q
{\displaystyle Q}
看做單位圓
S
1
{\displaystyle \mathbb {S} ^{1}}
,這樣
(
q
,
p
)
{\displaystyle (q,p)}
就位於圓柱體
S
1
×
R
{\displaystyle \mathbb {S} ^{1}\times \mathbb {R} }
上。取
(
q
,
p
)
∈
R
2
{\displaystyle (q,p)\in \mathbb {R} ^{2}}
只是因為
(
q
,
p
)
{\displaystyle (q,p)}
空間會更方便繪製。定義
z
(
t
)
=
(
q
(
t
)
,
p
(
t
)
)
T
{\displaystyle z(t)=(q(t),p(t))^{\mathrm {T} }}
、
f
(
z
)
=
(
p
,
−
sin
q
)
T
{\displaystyle f(z)=(p,-\sin q)^{\mathrm {T} }}
。讓我們用一些簡單的數值方法對這個系統進行積分。像往常一樣,選擇常數步長
h
{\displaystyle h}
,對任意非負整數
k
{\displaystyle k}
記
z
k
:=
z
(
k
h
)
{\displaystyle z_{k}:=z(kh)}
。
我們用以下方法:
z
k
+
1
=
z
k
+
h
f
(
z
k
)
{\displaystyle z_{k+1}=z_{k}+hf(z_{k})\,}
(顯式歐拉 );
z
k
+
1
=
z
k
+
h
f
(
z
k
+
1
)
{\displaystyle z_{k+1}=z_{k}+hf(z_{k+1})\,}
(隱式歐拉);
z
k
+
1
=
z
k
+
h
f
(
q
k
,
p
k
+
1
)
{\displaystyle z_{k+1}=z_{k}+hf(q_{k},p_{k+1})\,}
(辛歐拉);
z
k
+
1
=
z
k
+
h
f
(
(
z
k
+
1
+
z
k
)
/
2
)
{\displaystyle z_{k+1}=z_{k}+hf((z_{k+1}+z_{k})/2)\,}
(隱式中點法則)。
(注意,辛歐拉法用顯式歐拉法處理q ,用隱式歐拉法處理
p
{\displaystyle p}
。)
觀察到
H
{\displaystyle H}
在哈密頓方程的解曲線上是常數,於是可以描述系統的精確軌跡,是
p
2
/
2
−
cos
q
{\displaystyle p^{2}/2-\cos q}
的水平曲線 。在
R
2
{\displaystyle \mathbb {R} ^{2}}
中繪製了系統的精確軌跡和數值解。對顯式、隱式歐拉法,分別取
h
=
0.2
{\displaystyle h=0.2}
;z 0 = (0.5, 0)及(1.5, 0);對其他兩種方法,分別取
h
=
0.3
{\displaystyle h=0.3}
、z 0 = (0, 0.7);(0, 1.4)及(0, 2.1)。
單擺軌跡
顯式(或隱式)歐拉法是從原點向外(或向內)的螺旋運動。另兩種方法顯示了正確的定性行為,隱式中點法則與精確解的吻合程度高於辛歐拉法。
回顧一下,具有1自由度的哈密頓系統的精確流
ϕ
t
{\displaystyle \phi _{t}}
是保面積的,即
det
∂
ϕ
t
∂
(
q
0
,
p
0
)
=
1
{\displaystyle \det {\frac {\partial \phi _{t}}{\partial (q_{0},p_{0})}}=1}
for all
t
{\displaystyle t}
.
此式很容易手動驗證。對我們的單擺例子,可以發現,顯式歐拉法的數值流
Φ
e
E
,
h
:
z
k
↦
z
k
+
1
{\displaystyle \Phi _{{\mathrm {eE} },h}:z_{k}\mapsto z_{k+1}}
不 保面積;即
det
∂
∂
(
q
0
,
p
0
)
Φ
e
E
,
h
(
z
0
)
=
|
1
h
−
h
cos
q
0
1
|
=
1
+
h
2
cos
q
0
.
{\displaystyle \det {\frac {\partial }{\partial (q_{0},p_{0})}}\Phi _{{\mathrm {eE} },h}(z_{0})={\begin{vmatrix}1&h\\-h\cos q_{0}&1\end{vmatrix}}=1+h^{2}\cos q_{0}.}
隱式歐拉法也可進行類似計算,行列式為
det
∂
∂
(
q
0
,
p
0
)
Φ
i
E
,
h
(
z
0
)
=
(
1
+
h
2
cos
q
1
)
−
1
.
{\displaystyle \det {\frac {\partial }{\partial (q_{0},p_{0})}}\Phi _{{\mathrm {iE} },h}(z_{0})=(1+h^{2}\cos q_{1})^{-1}.}
辛歐拉法是 保面積的:
(
1
−
h
0
1
)
∂
∂
(
q
0
,
p
0
)
Φ
s
E
,
h
(
z
0
)
=
(
1
0
−
h
cos
q
0
1
)
,
{\displaystyle {\begin{pmatrix}1&-h\\0&1\end{pmatrix}}{\frac {\partial }{\partial (q_{0},p_{0})}}\Phi _{{\mathrm {sE} },h}(z_{0})={\begin{pmatrix}1&0\\-h\cos q_{0}&1\end{pmatrix}},}
於是
det
(
∂
Φ
s
E
,
h
/
∂
(
q
0
,
p
0
)
)
=
1
{\displaystyle \det(\partial \Phi _{{\mathrm {sE} },h}/\partial (q_{0},p_{0}))=1}
。隱式中點法則具有類似的幾何特性。
總結:單擺例表明,除顯式、隱式歐拉法不是解決問題的好方法外,辛歐拉法和隱式中點法則與系統的精確流非常吻合,後者更精確。而且後兩種方案與精確流都保面積,是幾何積分(實際上是辛積分 )的兩個例子。
Hairer, Ernst; Lubich, Christian; Wanner, Gerhard. Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations. Springer-Verlag. 2002. ISBN 3-540-43003-2 .
Leimkuhler, Ben; Reich, Sebastian. Simulating Hamiltonian Dynamics. Cambridge University Press. 2005. ISBN 0-521-77290-7 .
Budd, C.J.; Piggott, M.D. Geometric Integration and its Applications. Handbook of Numerical Analysis 11 . Elsevier. 2003: 35–139. ISBN 9780444512475 . doi:10.1016/S1570-8659(02)11002-7 .
Kim, Pilwon. Invariantization of Numerical Schemes Using Moving Frames. BIT Numerical Mathematics 47 (3). Springer. 2007: 525–546. doi:10.1007/s10543-007-0138-8 .