自20世纪末期, 对称约化理论提出以来, 一直是几何动力学领域的热门课题之一, 并广泛应用于力学、物理学和工程科学等多个领域, 如固体力学、流体力学、场论、量子力学和广义相对论等^[1]。对称约化的主要目的是利用系统的对称性消除动力学系统的部分局部坐标变量, 从而简化实际的动力学系统。对称约化的思想最早源于Routh^[2]于1884年利用循环坐标对系统进行化简, 这种化简理论对应于现代的Lagrange约化理论。现代约化理论始于20世纪六七十年代Arnold^[3-4]和Smale^[5]的重要工作, Arnold开展了Lie群约化方法, Smale引入动量映射概念。Marsden等^[6]在前人工作的基础上, 利用等变动量映射开展辛流形上的约化理论, 使现代约化理论走向成熟, 并得到广泛应用^[7-8]。目前, 完整约束力学系统的几何动力学约化理论已经非常成熟, 大体上包括:辛约化^[9]、Poission约化^[1]和Lagrange约化^{[1, 10]}。对称性约化理论在研究约束系统几何动力学及应用以及约束系统的保结构算法中发挥着非常重要的作用, 为研究约束系统几何数值积分的几何不变性质提供了新的途径。但是, 在现有的少量研究工作中, 对称约化理论在约束力学系统几何数值积分的研究中, 还没有发挥其应有的作用。较复杂的动力学系统经过对称约化得到一个简单的新的动力学系统, 系统方程的性质发生了变化, 变量所代表的物理含义有时也发生改变, 那么约化系统与原系统的内部联系, 特别是一些重要的性质是否能够很好地得以保留, 对原系统和约化后的系统分别做数值积分, 对称约化对数值积分的结果是否会受到影响?本文就这一问题展开探讨, 通过数值实验验证对称约化对数值积分的影响。首先, 简要介绍完整约束系统的辛约化理论; 然后, 简要介绍完整约束系统的数值积分方法, 并应用数值积分方法计算约化前和约化后系统的动力学方程, 比较约化前后动力学变量的数值结果; 最后得出结论。

1 完整约束系统的辛约化

完整约束力学系统总可以表示为Hamilton形式。定义Hamilton系统为$(M, \; \mathit{\Omega}, \; H)$, 这里$(M, \; \mathit{\Omega})$是辛流形, $\mathit{\Omega} $是辛2-次型, $H$是定义在流形M上的Hamilton函数。系统的运动微分方程可以表示成如下形式:

$\frac{{{\rm{d}}x}}{{{\rm{d}}t}} = {}^\# {\rm{d}}H(x),$

(1)

这里${}^\# {\rm{d}}H$表示Hamilton矢量场。辛流形$(M, \; \mathit{\Omega})$在Lie群$G$的作用下, 可以表示为$\phi :G \times M \to M$, 或者, 对任意的$g \in G$和$x \in M$, 映射$\phi $可以局域的表示为

$ \phi (g, \;z) = {\phi _g}(z), $

(2)

如果对于任意一个$g \in G$, Lie群$G$在辛流形$(M, \mathit{\Omega})$上的作用$\phi $满足如下等式:

$\phi _g^ * \mathit{\Omega} = \mathit{\Omega}, $

(3)

则该作用${\phi _g}$是辛的。

定义(辛约化)  辛流形$(M, \mathit{\Omega})$的辛约化是$M$的子流形(或侵入子流形)$N$到辛流形$(P, {\mathit{\Omega} _P})$的侵没满射: $\pi :N \to P$, 满足如下方程:

${\pi ^ * }{\mathit{\Omega} _P} = {\mathit{\Omega} _N}, $

(4)

这里${\mathit{\Omega} _N}=i_N^ * \mathit{\Omega} $是流形$N$上的$\mathit{\Omega} $诱导的2-形式, $(P, \; {\mathit{\Omega} _{\; P}})$称为约化辛流形。如果$N$是余迷向子流形, 可以称$\pi $是辛约化。

辛约化满足Marseden-Weinstein-Meyer辛约化定理^{[2, 6]}。取点$x \in M, \mu \in {\mathit{\mathcal{G}}^ * }$是辛流形$M$和Lie代数上$\mathcal{G}^ *$的点, 其在对称性Lie群G作用$\phi $下的轨迹分别为$G \cdot x$和${\mathcal{O}_\mu }$, 它们对应的迷向群分别为${G_x}$和${G_\mu }$, 对应的Lie代数分别为${\mathcal{G}_x }$和${\mathcal{G}_\mu }$, 则可以得到如下的辛约化定理。

定理   设Lie群$G$作用在辛流形$(M, \mathit{\Omega})$上, 并且该作用有一个等变动量映射

$J:M \to {\mathcal{G}^ * }, $

如果存在正则值

$\mu \in J(M) \subset {\mathcal{G}^ * }, $

且存在伴随余迷向子群

${G_\mu } = \left\{ {g \in G\left| {Ad_g^ * \mu = \mu } \right.} \right\}, $

则${G_\mu }$自由而正常的作用于

${M_\mu }(x) = {J^{ - 1}}(\mu )\pi (x)$

上; 如果存在内映射

${i_\mu }:{J^{ - 1}}(\mu ) \to M$

和正则投影映射

${\pi _\mu }:{J^{ - 1}}(\mu ) \to {\hat M_\mu }, $

则一定存在具有辛结构${\mathit{\hat \Omega_\mu} }$的商流形

${\hat M_\mu } = {J^{ - 1}}(\mu )/{G_\mu }, $

满足

$\pi _\mu ^ * \circ \mathit{\hat \Omega _\nu } = i_\mu ^ * \circ \mathit{\Omega} 。$

并且存在一个侵入映射$\varphi $使得

${\hat M_\mu } = {J^{ - 1}}(\mu )/{G_\mu }$

是$M/G$的一个子流形。这里$M/G$由自然投影映射

$\pi :M \to M/G$

确定。定理证明参见文献[2, 6]。

这个定理保证了约化后的商空间${\hat M_\mu }$是一个辛流形, 同时也给出进行辛约化的方法, 即利用对称性Lie群作用于辛流形, 从而找到有效的对称性与守恒量, 进行对称约化。

2 完整约束系统的辛几何算法

完整约束系统总可以表示为Hamilton正则方程(1)的形式, 如取Hamilton函数$H$为$2n$个变量(p₁, …, p_n; q¹, …, qⁿ)的可微函数, 并令

$ z = {\left( {{p_1}, {\rm{ }} \ldots, {p_n};{q^1}, {\rm{ }} \ldots, {q^n}} \right)^T}, $

则方程(1)可以表示成如下坐标形式:

$\frac{{{\rm{d}}z}}{{{\rm{d}}t}} = {\mathit{\boldsymbol{J}}^{ - 1}}\frac{{\partial H(z)}}{{\partial z}}, $

(5)

这里, $\mathit{\boldsymbol{J}}=\left ({\begin{array}{*{20}{c}} {\bf{0}} & {{\mathit{\boldsymbol{I}}_n}}\\ { -{\mathit{\boldsymbol{I}}_n}} & {\bf{0}} \end{array}} \right)$是标准辛矩阵, ${\mathit{\boldsymbol{I}}_n}$和0是n阶单位矩阵和零矩阵。基于Hamilton力学的基本原理:正则方程(5)的解由一个单参数辛群$\left\{ {g_H^t} \right\}$生成, 冯康等^[11]提出Hamilton方程的辛算法:基于恒等变换的生成函数, 可以得到一个辛变换, 从而构造辛差分格式, 如2阶的Euler中点格式:

$\begin{array}{l} \frac{{{z^{i + 1}} - {z^i}}}{\tau } = {J^{ - 1}}{\left( {\frac{{\partial H}}{{\partial z}}} \right)_{\frac{{{z^{i + 1}} + {z^i}}}{2}}}\\ (i = 0, \;1, \;...\;, \;n - 1), \end{array}$

(6)

这里$\tau $是时间步长, 还可以在此基础上构造具有更高精度的差分格式, 如4阶辛差分格式^[12]。

3 数值算例和结论

众所周知, 平面Kepler问题是典型的二体问题, 且满足机械能守恒和角动量守恒。令

$q = ({q_1}, \;{q_2})$

为该二体问题的广义坐标, 并取

$p = ({p_1}, \;{p_2}) = ({\dot q_1}, \;{\dot q_2}), $

则可以得到系统的Hamilton函数:

$H({p_1}, \;{p_2}, \;{q_1}, \;{q_2}) = \frac{1}{2}(p_1^2 + p_2^2) - \frac{1}{{\sqrt {q_1^2 + q_2^2} }}, $

从而得到系统的Hamilton方程:

$\left\{ \begin{array}{l} {{\dot p}_1} = - \frac{{{q_1}}}{{{{(q_1^2 + q_2^2)}^{3/2}}}}, \;\;\;\;{{\dot q}_1} = {p_1};\\ {{\dot p}_2} = - \frac{{{q_2}}}{{{{(q_1^2 + q_2^2)}^{3/2}}}}, \;\;\;\;{{\dot q}_2} = {p_2}。\; \end{array} \right.$

(8)

在Lie群变化下, 方程(8)可以约化为一维谐振子的运动方程:

$ {u''_{}} + u = 1/d, $

(9)

这里, $u=1/r, $r和θ是对应的极坐标系下的坐标, 且满足

$\begin{array}{l} {q_1} = r\cos \theta, \;\;{q_2} = r\sin \theta 。\\ d = L_0^2 = {({q_1}{{\dot q}_2} - {q_2}{{\dot q}_1})^2} = {({r^2}\dot \theta )^2} \end{array}$

是系统的角动量,

″=d²/dθ ²。

令${p_u}=u'$, 则与式(7)等价的对应方程(9)的Hamilton函数为

$H({p_u}, \;u) = \frac{1}{2}(p_u^2 + {u^2}) - \frac{u}{d} - \frac{{{H_0}}}{d}, $

(10)

这里,

${H_0} = \frac{1}{2}(p_1^2 + p_2^2) - \frac{1}{{\sqrt {q_1^2 + q_2^2} }} = \frac{1}{2}({\dot r^2} + {r^2}{\dot \theta ^2}) - \frac{1}{r}, $

是原二体系统的机械能。从而得到方程(9)对应的Hamilton方程:

${p'_u} = - u + 1/d, {\rm{ }}u' = {p_u}。$

(11)

取如下初始条件:

$\begin{array}{l} {q_1}(0)1 - e, \;\;{q_2}(0) = 0, \;\;{{\dot q}_1}(0) = 0, \\ {{\dot q}_2}(0) = \sqrt {(1 + e)/(1 - e)}, \;\;\;e = 0.6。 \end{array}$

从而系统的总能量

${H_0} = - 1/2, $

(13)

总角动量

${L_0} = \sqrt {1 - {e^2}} = 0.8, \;\;d = L_0^2。$

我们采用上面的辛差分格式(式(6)), 并取步长$h=0.001$, 分别计算方程(8)和(11), 并比较数值结果。图 1给出系统的总能量和能量误差, 图 2给出系统的总动量和动量误差。

从图 1和2可以看出, 虽然在原体系框架下计算得到的总能量没有在约化后的体系中得到的结果好, 但是总动量的数值结果却优于约化体系中的数值结果。因此, 在约化前体系下和约化后体系下的数值结果并没有本质上的区别。并且, 我们求得的系统的总能量和总动量都在数值计算的误差允许范围内, 数值结果都很好地保持了原有系统的守恒量。区别在于, 我们在约化后的体系下进行数值研究, 由于体系的自由度减少, 从而减少了编程计算的难度, 同时也减少计算机机时。因此, 对于完整系统, 无论是在原体系框架下, 还是在约化后的系统中, 都可以得到满意的数值结果, 但在约化体系下, 可以获得简洁的计算机程序, 并减少计算机工作时间。这对于复杂的完整动力学系统来说, 对动力学系统进行对称约化, 以减少动力学系统的自由度数, 可以有效地进行数值计算, 从而在约化系统下间接地研究原系统的动力学性质。