在常微分方程 的数值计算中,几何积分 是一种保留微分方程的流的精确几何特性的数值方法。
可考虑单摆 运动以引出几何积分的研究。
设摆锤质量为
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 .