COMPUTATION OF INVARIANT MANIFOLD BASED ON SYMPLECTIC ALGORITHM OF MIXED LIE OPERATOR
-
摘要: 平动点附近周期轨道的不变流形因其在低能轨道转移中起着重要作用而受到广泛关注.在设计低能轨道过程中不变流形要实时进行能量匹配，但利用传统数值积分方法进行积分时能量会耗散.显式辛算法具有比隐式辛算法计算效率高的优势，但其要求Hamilton系统必须分成两个可积的部分，而旋转坐标系下的圆型限制性三体问题是不可分的，因而显式辛算法难以用于求解旋转坐标系下的圆型限制性三体问题.本文通过引入混合Lie算子，成功实现了带三阶导数项的力梯度辛算法对圆型限制性三体问题的求解，并将基于混合Lie算子的带三阶导数项的辛算法与Runge-Kutta78算法和Runge-Kutta45算法进行仿真对比，仿真结果表明基于混合Lie算子的含有三阶导数项的辛算法位置精度高、能量误差小且计算效率高.利用基于混合Lie算子的带三阶导数项的辛算法计算不变流形，可以实现低能轨道转移过程中轨道拼接点的能量精准匹配.Abstract: Invariant manifolds of periodic orbit near the libration points attract a lot of attentions due to their importance in the low-energy orbits transfer problem. In the process of low-energy orbit design, the energy of the invariant manifolds must be matched, but the energy is dissipated when integrating with traditional numerical integration method. The explicit symplectic algorithm with energy conservation is more efficient than the implicit symplectic algorithm, but it requires the Hamiltonian system to be divided into two integral parts, while the circular restricted three-body problem in the rotating coordinate system being inseparable. It is difficult to solve the circular restricted three-body problem in the rotating coordinate system by explicit symplectic algorithm. In this paper, the mixed Lie derivative operator of kinetic energy is used to solve the circular restricted three-body problem in the rotating coordinate system, and the effectiveness of this explicit symplectic algorithm with the third derivation in dealing with this problem has been showed. Compared with the Runge-Kutta45 method and Runge-Kutta78 method, the symplectic algorithm with the third-order derivative term not only has high precision but also the smallest energy error and the highest efficiency. Finally, the invariant manifolds are calculated by the symplectic algorithm with the third derivative term, the patched point can match accurately with this method.
-
引言
平动点[1]是限制性三体问题中的动平衡点，在随两个天体一起旋转的参考系中处于静止状态.它们是几何意义上的点，其本身的应用价值十分有限.真正引人关注的是其周围大量存在的周期轨道[2]，以及与之相联的管状不变流形[3](包括稳定流形和不稳定流形).周期轨道是完成某些特殊任务的理想场所，而不变流形则提供了往返于主天体和周期轨道之间的低能耗转移途径[4-9].太阳系中不同三体系统平动点周期轨道的不变流形还存在相交的情形，从而拓展了系统间的深空转移轨道设计范围，因此寻找航天器从地球停泊轨道到平动点附近周期轨道的转移轨道显得尤为重要.
20世纪90年代，Gómez等[10]引入非线性动力系统理论，利用稳定流形将航天器送到了平动点附近，是转移轨道设计[11]的重大突破.这一发现不但节省了能量也为星际超级公路理论的诞生奠定了基础.但是不论是日-地系统还是地-月系统，不变流形距离地球都比较远，需要在推力或引力辅助[12]作用下将航天器从地球停泊轨道转移到稳定流形上.由于在优化拼接点时需要多次计算不变流形，如何选取不变流形的拼接点是进行低能轨道转移设计的关键.已经有学者对不变流形的计算进行研究，Zhang等[13]利用三次回旋插值法计算圆型限制性三体问题不稳定平动点周期轨道的不变流形，Lei等[14-17]利用Lindstedt-Poincaré方法求解了限制性三体和四体平动点周期轨道不变流形的高阶近似解析解，Dellnitz等[18]提出了方向集数值计算方法，Howell等[19]利用胞映射法[20]描述和存储不变流形的数据.这些计算方法提高了不变流形的计算速度，但是却没有关注不变流形的能量变化.由于圆型限制性三体问题是一个非线性自治哈密顿系统，没有严格的解析解，而且由于强非线性其对初值误差十分敏感，较小的计算误差将导致微分方程的解较大地偏离实际情况.因此需要寻找一种保能量的不变流形计算方法.
普通数值解法具有人为耗散性，长时间计算会导致系统总能量的耗散随时间的平方增长.在研究三体问题的过程中，对辛结构的破坏无疑会对动力学特性造成很大影响，使得系统的长期演化不能真实反映客观事实.辛积分方法[21-22]由于具有保辛结构和能量守恒两个重要特性，近年来得到了快速发展[23-26].显式辛算法具有比隐式辛算法计算效率高的优势，但其要求Hamilton系统必须可以分成两个可积的部分.旋转坐标系下的圆型限制性三体问题由于受Coriolis力的影响，不能像惯性系下的系统一样可分，因此，从形式上看显式差分格式对其不适用.廖新浩等[27]将圆型限制性三体问题的Hamilton函数分成动能与势能的和以及由坐标旋转产生的非惯性力两部分，然后利用算法合成构造差分结构，计算复杂. Sun等[28]分析了带有三阶导数项的力梯度辛算法在保能量上的优势，但是没有将该算法应用在旋转坐标系下的圆型限制性三体问题中，也没有将辛算法用在不变流形的求解上.
鉴于显式辛算法不适合求解旋转坐标系下圆型限制性三体问题，本文通过重新定义类动能的Lie算子，严格证明基于混合Lie算子的含有三阶导数项的显式辛算法可以求解旋转坐标系下的圆型限制性三体问题.首先，给出文中所用的动力学模型和不变流形的计算方法；然后，建立混合Lie算子并研究显式辛算法在求解Hamilton系统类动能不具有标准二次型的圆型限制性三体问题时的可行性；最后，将本文改进算法与Runge-Kutta78(RK78) 和Runge-Kutta45(RK45)[29]算法进行仿真对比.
1. 空间圆型限制性三体问题
1.1 动力学模型
圆型限制性三体问题描述的是质量为$m_3 $的航天器$P_3 $在两个绕着共同质心做匀速圆周运动的主天体$P_1 $和$P_2 $的引力场下运动，其质量分别为$m_1$和$m_2$，且$m_3 \ll m_1 \ll m_2$，以$P_1$和$ P_2 $的质心为原点，$P_1 $指向$P_2 $的方向为$x$轴，$y$轴在两个主引力体旋转平面上，$z$轴垂直于$x y$平面，满足右手法则，如图 1所示.
假设$P_3 $在$oxy$平面运动，为了简化模型，对单位进行无量纲化处理.令引力常量$G$、主天体$P_{1}$和$P_{2}$之间的距离、旋转角速度$\omega $、质量和均为1，因此$P_1 $的质量为$1 -\mu = {m_1 } /{(m_1 + m_2 })$，位置坐标为$(-\mu, 0, 0)$，$P_2 $的质量为$\mu = {m_2 } /{(m_1 + m_2 })$，位置坐标为$(1 -\mu, 0, 0)$.文中所有量如无特殊说明都是无量纲的.在此旋转坐标系下航天器的运动方程和Largrange函数[30]分别为
$$ \left.\begin{array}{l} \ddot {x} - 2\dot {y} = \varOmega _x \\ \ddot {y} + 2\dot {x} = \varOmega _y \\ \ddot {z} = \varOmega _z \end{array} \right \} $$ (1) 和
$$ L(x, y, z, \dot x, \dot y, \dot z)=\dfrac 12 \big [(\dot x-y)^2+(\dot y+x)^2+\dot z^2 \big]+ \\ \qquad \varOmega (x, y, z) $$ $\varOmega (x, y, z)$是系统的势函数，通常表示为
$$ \varOmega (x, y, z) = {(x^2 + y^2)} / 2 + {(1 - \mu )} / {r_1 } + \mu / {r_2 } $$ 其中$r_1, r_2 $分别表示航天器与$P_1, P_2 $之间的距离
$$ r_1 = \sqrt {(x + \mu )^2 + y^2 + z^2}, r_2 = \sqrt {(x - 1 + \mu )^2 + y^2 + z^2} $$ 通过Legendre变换
$$ p_i = {\partial L} /{\partial \dot {q}_i } \ \ \ (i = 1, 2, 3) \\ H(q_i, p_i ) = \sum\limits_{i = 1}^n {p_i } \dot {q}_i - L(q_i, p_i ) $$ 其中$q_i, p_i \ (i = 1, 2, 3)$分别为航天器的坐标${\pmb q}= (x, y, z)$和动量${\pmb p} =(p_x, p_y, p_z)$，可以得到圆型限制性三体问题的Hamilton函数
$$ H(x, y, z, p_x, p_y, p_z ) = \dfrac{1}{2}(p_x^2 + p_y^2 + p_z^2 ) + \\ \qquad yp_x - xp_y + \varOmega (x, y, z) $$ (2) 该系统具有雅可比能量积分
$$ C(x, y, z, \dot {x}, \dot {y}, \dot {z}) = 2\varOmega (x, y, z) - (\dot {x}^2 + \dot {y}^2 + \dot {z}^2) $$ 由Hamilton函数(2)，可推导出圆型限制性三体问题的正则方程为
$$ \begin{array}{ll} \dot {x} = \dfrac{\partial H}{\partial p_x } = p_x + y, & \dot {p}_x = - \dfrac{\partial H}{\partial x} = p_y - x + \varOmega _x \\ \dot {y} = \dfrac{\partial H}{\partial p_y } = p_y - x, & \dot {p}_y = - \dfrac{\partial H}{\partial y} = - p_x - y + \varOmega _y \\ \dot {z} = \dfrac{\partial H}{\partial p_z } = p_z, & \dot {p}_z = - \dfrac{\partial H}{\partial z} = \varOmega _z \end{array} $$ 1.2 不变流形
不变流形是与平动点周期轨道光滑连接的一簇空间轨道，航天器在不变流形上可以无动力飞行.因此，不变流形可以作为低能轨道转移的通道，在深空探测轨道转移设计中有着重大应用价值.
将动力学方程(1) 在周期轨道上任意一点$\bar {x}_{0} $处线性化，则有
$$ \Delta \dot {\bar {\pmb x}} = A(t)\Delta \bar {\pmb x} $$ 其中$A(t)$为系统(1) 的雅可比矩阵，$\Delta \bar {\pmb x}$为相对于不动点状态的偏移量.由$\bar {\pmb x}_{0} $点出发，积分一个周期后的状态所对应的状态转移矩阵${\pmb\varPhi}(0, T)$称为单值矩阵[31].周期轨道上不同的点对应不同的单值矩阵，每个单值矩阵有一个不稳定特征根$\lambda _s (\lambda _s > 1)$及一个稳定特征根$\lambda _u = 1 / {\lambda _s }(\lambda _u > 1)$，它们分别对应着特征向量$\bar {\pmb v}_s $和$\bar {\pmb v}_u $，它们包含了稳定流形$W_{L_j, po}^s $和不稳定流形$W_{L_j, po}^u $的方向信息.在$\bar {\pmb x}_{0} $点处沿特征向量方向增加一个微小状态扰动，则可以得到不变流形的积分初值.对于不变流形的计算有
$$ \bar {\pmb x}_{i, s/ u} = \bar {\pmb x}_i \pm {\varepsilon} \bar {\pmb v}_{s / u} $$ 其中，$\bar {\pmb x}_{i, s/u} $为不变流形上一点，$\bar {\pmb x}_i $为周期轨道上任意点，$\varepsilon $为微小标量，代表扰动的大小，其取值非常关键，取值过小，流形的计算速度太慢，过大则无法保证其精度，一般取$10^{ -6}$[32]，$\bar {\pmb v}_{s / u} $为特征向量，代表着不变流形的方向.
一般都采用数值积分来计算不变流形.不变流形计算的关键是获得积分状态初值，而此初值和Halo轨道上的点有密切关系.由于圆型限制性三体问题的强非线性使得微分方程组(1) 的解对初始值十分敏感，可能较小的计算误差都将导致不变流形的能量误差巨大.辛算法因具有保能量的特点而受到广泛关注，而显式辛算法具有比隐式辛算法效率高的优势，因此本文考虑利用含有三阶导数项的力梯度辛算法来计算不变流形.力梯度辛算法是将Hamilton系统分解成动能$T({\pmb p})$(仅关于动量${\pmb p} = (p_x, p_y, p_z)$的正定二次函数)和势能$V(q)$(关于坐标${\pmb q} = (x, y, z)$的函数)两个可积部分，然后分别求解.圆型限制性三体问题的Hamilton函数有多种分解形式，最简单的是分解成类动能$T({\pmb p})$和势能$V({\pmb q})$两部分
$$ \left.\begin{array}{l} T({\pmb p}) = \dfrac{1}{2}(p_x^2 + p_y^2 + p_z^2 ) + yp_x - xp_y \\ V({\pmb q}) = - \varOmega (x, y, z) \end{array} \right\} $$ (3) 其中，$yp_x -xp_y $是由于Coriolis力的影响而产生的偏移部分.从式(3) 可以看出，类动能$T({\pmb p})$不具有标准的二次型形式，因此从形式上看力梯度辛算法不适合求解圆型限制性三体问题.后文试图通过改进算法来解决这一难题.
2. 基于混合Lie算子的力梯度辛算法
考虑可以分解为动能$T({\pmb p})$和势能$V({\pmb q})$的Hamilton系统
$$ H({\pmb p}, {\pmb q}) = T({\pmb p}) - V({\pmb q}) $$ 其中，动能$T({\pmb p})$仅仅是关于动量$p$的正定二次函数, 即$T({\pmb p}) = {{\pmb p}^2}/2$，势能$V({\pmb q})$是关于坐标$q$的函数.分别定义$T({\pmb p})$和$V({\pmb q})$的Lie算子[33-34]
$$ {\rm A} = \sum\limits_{i = 1}^3 {\dfrac{\partial T}{\partial p_i }\dfrac{\partial }{\partial q_i }}, \ \ \ {\rm B} = \sum\limits_{i = 1}^3 {\dfrac{\partial V}{\partial q_i }\dfrac{\partial }{\partial p_i }} $$ (4) 由于圆型限制性三体问题Hamilton函数分解中类动能$T({\pmb p})$不是动量${\pmb p}$的标准二次型，因此具有形式(4) 的Lie算子不再适用，考虑将类动能$T({\pmb p})$的Lie算子变为位置与动量的混合型算子
$$ \bar{\rm A} = \sum\limits_{i = 1}^3 \Big ( \dfrac{\partial T}{\partial p_i }\dfrac{\partial }{\partial q_i } - \dfrac{\partial T}{\partial q_i }\dfrac{\partial }{\partial p_i } \Big) = (p_x + y)\dfrac{\partial }{\partial x} + p_y \dfrac{\partial }{\partial p_x } +\\ \qquad (p_y - x)\dfrac{\partial }{\partial y} - p_x \dfrac{\partial }{\partial p_y } + p_z \dfrac{\partial }{\partial z} $$ 它作用于方程组(3) 第一个方程的位置坐标和动量得到微分方程组
$$ \begin{array}{ll} \dot {x} = \bar {\rm A}x = p_x + y, & \dot {p}_x = \bar {\rm A}p_x = p_y - x \\ \dot {y} = \bar {\rm A}y = p_y - x, & \dot {p}_y = \bar {\rm A}p_y = - p_x \\ \dot {z} = \bar {\rm A}z = p_z, & \dot {p}_z = \bar {\rm A}p_z = 0 \end{array} $$ 该方程组的分析解为
$$ \begin{array}{l} x = x_0 \cos t + y_0 \sin t + t(p_{x0} \cos t + p_{y0} \sin t) \\ y = y_0 \cos t - x_0 \sin t + t(p_{y0} \cos t - p_{x0} \sin t) \\ z = p_{z0} t \\ p_x = p_{x0} \cos t + p_{y0} \sin t \\ p_y = p_{y0} \cos t - p_{x0} \sin t \\ p_z = p_{z0} \end{array} $$ 其中, $(x_0, y_0, z_0, p_{x0}, p_{y0}, p_{z0})$为初始状态，$t$是积分时间，于是新定义的混合算子$\bar {\rm A}$可积，说明重新定义的混合Lie算子是合理的.容易验证算子${\rm B}$也可积.
利用混合Lie算子$\bar {\rm A}$和${\rm B}$构造复杂算子
$$ {\rm C}=[\rm B, [\bar A, B]] \\ {\rm D} = [{\rm B}, [\bar {\rm A}, [{\rm B}, [\bar {\rm A}, {\rm B}]]]] =\rm 2{\rm B}{\rm B}\bar {\rm A}{\rm B}\bar {\rm A} - 2\bar {\rm A}{\rm B}\bar {\rm A}{\rm B}{\rm B} -\\ \qquad {\rm B}\bar {\rm A}\bar {\rm A}\rm BB + BBB\bar {\rm A}\bar {\rm A} + \bar {\rm A}\bar {\rm A}BBB - BB\bar {\rm A}\bar {\rm A}B =\\ \qquad - 2\sum\limits_{i = 1}^3 \sum\limits_{k = 1}^3 \sum\limits_{j = 1}^3 V_k (2V_{jk} V_{ij} + V_j V_{ijk} ) \dfrac{\partial }{\partial p_i } =\\ \qquad - 4\sum\limits_{i = 1}^3 \sum\limits_{k = 1}^3 \sum\limits_{j = 1}^3 V_k V_{jk} V_{ij} \dfrac{\partial }{\partial p_i } -\\ \qquad 2\sum\limits_{i = 1}^n \sum\limits_{k = 1}^n \sum\limits_{j = 1}^n V_k V_j V_{ijk} \dfrac{\partial }{\partial p_i } $$ 算子$\rm D$是含有一阶、二阶和三阶导数项的算子，可以与混合Lie算子$\bar {\rm A}$和Lie算子$\rm B$一起构造高阶算法
$$ \begin{array}{l} e^{hW} = e^{ah\bar {\rm A}} \otimes e^{bh\rm B} \otimes e^{dh\bar {\rm A}} \otimes e^{\varepsilon h{\rm B} + fh^3{\rm C} + gh^5{\rm D}}\otimes \\ \qquad e^{dh\bar {\rm A}} \otimes e^{bh\rm B} \otimes e^{ah\bar {\rm A}} \end{array} $$ 其中, $h$是积分步长，各积分系数分别为
$$ a = 0.070 310 871 341 794 26 \\ b = 0.225 673 220 373 456 00 \\ d = 0.429 689 128 658 206 00 \\ \varepsilon = 0.548 653 559 253 087 0 \\ f = 8.247 384 758 086 507 00 \\ g = 0.000 006 164 365 178 93 $$ 其差分格式如下
$$ \begin{array}{l} {\pmb q}^\ast = {\pmb q}_{n - 1} + ah{\pmb p}_{n - 1} \\ {\pmb p}^\ast = {\pmb p}_{n - 1} + bh{\pmb f} ({\pmb q}_{n - 1} ) \\ {\pmb q}^\dagger = {\pmb q}\ast + dh{\pmb p}^\ast \\ {\pmb p}^\dagger = {\pmb p}^\ast + \varepsilon hf({\pmb q}^\ast ) + fh^3\nabla \left| {f({\pmb q}^\ast )} \right|^2 + gh^5\nabla \left| {f({\pmb q}^\ast )} \right|^4 \\ {\pmb q}^\diamondsuit = {\pmb q}^\dagger + dh{\pmb p}^\dagger \\ {\pmb p}_n ={\pmb p}^\dagger + bh{\pmb f} ({\pmb q}^\diamondsuit ) \\ {\pmb q}_n = {\pmb q}^\diamondsuit + ah{\pmb p}_n \end{array} $$ 这里，力${\pmb f} = {\partial V} / {\partial {\pmb q}}$.
由
$$ {\rm B}{\rm B}\bar {\rm A}p_i = 0, \ \ \bar {\rm A}{\rm B}{\rm B}q_i = 0, \ \ {\rm B}\bar {\rm A}{\rm B}q_i = 0 \\ {\rm B}\bar {\rm A}{\rm B}p_i = \sum\limits_{j = 1}^3 \sum\limits_{k = 1}^3 V_{ij} V_k T_{jk} $$ 具体证明过程见附录A, 得到
$$ \begin{array}{l} {\rm D} = - 4\sum\limits_{i = 1}^3 \sum\limits_{k = 1}^3 \sum\limits_{j = 1}^3 V_k V_{jk} V_{ij} \dfrac{\partial }{\partial p_i } - \\ \qquad 2\sum\limits_{i = 1}^n \sum\limits_{k = 1}^n \sum\limits_{j = 1}^n V_k V_j V_{ijk} \dfrac{\partial }{\partial p_i } = - 4{\rm E} - 2{\rm F} \end{array} $$ 式中，${\rm E}$和${\rm F}$的具体表达形式见附录B.这表明算子${\rm D}$是两主天体的引力梯度，而不是旋转坐标系所产生的非惯性力与主天体引力的混合梯度.
因此，通过引入混合Lie算子，分解成形式(3) 的圆型限制性三体问题可以用改进的显式辛算法进行求解.
3. 仿真算例
仿真使用Matlab R2008b软件，计算机为Windows 7系统，配置Intel(R) Core (TM)2 Duo CPU E7500处理器，主频2.93 GHz，内存2.00GB，32位操作系统.利用第2节所得到的改进显式辛算法MTGS研究地-月圆型限制性三体问题并与RK78和RK45法进行仿真对比，初始值参数如表 1所示.
表 1 初始值参数Table 1. The initial value of the parameter
积分步长取0.001，积分两个周期得到的轨道如图 2所示，可以看出本文改进的算法和RK78算法得到的轨道仍然能够闭合，而通过RK45算法得到的轨道开始偏离周期轨道.当积分时长为3个周期时，轨道如图 3所示，改进的显式辛算法、RK78算法与RK45算法都不能保证轨道闭合，RK45算法发散速度最快，改进的算法发散速度最慢.虽然Halo轨道产生偏移的快慢程度与其自身稳定性也有一定关系，但是，在同样时间内，以同样步长计算，不同算法发生偏移的快慢反映了算法的精确度.
![]()
图 2 三种算法积分两个周期得到的轨道Figure 2. The orbits of three algorithms integral two orbit periods![]()
图 3 三种算法积分3个周期得到的轨道Figure 3. The orbits of three algorithms integral three orbit periods积分100 000步时绝对能量误差如图 4所示，可以看出当步长取固定值0.001，积分100 000步时改进的显式辛算法的能量误差最大为${10}^{-9}$量级.虽然RK78算法的初始能量误差比较小，但随着时间推移，能量误差积累变得越来越大，出现骤增的情况，而RK45算法的能量误差始终很大.与之相比，显式辛算法的能量误差始终保持在一定范围内，能量误差突然增大是因为航天器与主天体$P_{2}$的距离突然变得非常小.
![]()
图 4 三个算法积分100 000步的绝对能量误差Figure 4. The absolute energy errors of 100 000 steps with three algorithms利用改进的显式辛算法、RK78和RK45算法分别计算由表 1所示初始参数得到的Halo轨道对应的一条不变流形，综合考虑积分时间和积分精度，积分步长取0.000 1，积分100 000步时三种算法得到的稳定流形如图 5所示，不稳定流形与其关于$y$轴对称，3种算法得到的不变流形从形状来看差别无几，但是不重合.从图 6所示能量相对误差分析，可知利用改进的算法得到的不变流形，可以保证在长时间积分过程中能量误差在很小范围内，RK78算法初始能量误差比较小，但呈增大趋势，而RK45的能量误差一直都很大.这3种算法所耗时长和相对能量误差的最大值如表 2所示，改进的显式辛算法耗时最短，计算效率约是RK78的12.308 4倍，且相对能量误差的最大值最小，因此本文改进的显式辛算法在低能轨道转移过程中可以实现轨道拼接点精准能量匹配.
![]()
图 6 三种算法积分100 000步的相对能量误差Figure 6. The relative energy error of 100 000 steps with three algorithms表 2 三种算法的积分效率Table 2. Eefficiency of three algorithms
4. 结论
在基于流形拼接法设计航天器低能转移轨道的过程中，需要实时进行不变流形的能量计算.本文针对传统数值积分能量耗散问题，研究了显式辛算法在旋转坐标系下的圆型限制性三体问题不变流形的保能量问题.显式辛算法长时间积分时没有长期的变化，能够保持Hamilton系统的辛结构，本文通过改变圆型限制性三体问题分解形式中类动能(动能部分不是严格的正定二次型)Lie算子的形式，证明含有一阶导数项、二阶导数项和三阶导数项的改进显式辛算法可以求解旋转坐标系下的圆型限制性三体问题，这也说明改进的Lie算子形式合理.最后，应用改进的显式辛算法、RK78和RK45算法分别计算了地-月圆型限制性三体问题L1平动点附近的周期轨道和不变流形，并对其能量误差和计算效率进行了分析.仿真结果显示，本文改进的显式辛算法在精度和能量误差两方面都具有较好的优势.
附录A
$$ \begin{array}{l} {\rm B}{\rm B}\bar {\rm A}x = {\rm B}{\rm B}\Big [(p_x + y)\dfrac{\partial }{\partial x} + (p_y- x)\dfrac{\partial }{\partial y} +\\ \qquad p_z \dfrac{\partial }{\partial z} + p_y \dfrac{\partial }{\partial p_x }-p_x \dfrac{\partial }{\partial p_y }\Big]x = {\rm B}{\rm B}(p_x + y) =\\ \qquad {\rm B}\Big [-\Big (\dfrac{\partial }{\partial x}\dfrac{\partial }{\partial p_x } + \dfrac{\partial }{\partial y}\dfrac{\partial }{\partial p_y } + \dfrac{\partial }{\partial z}\dfrac{\partial }{\partial p_z }\Big)\Big](p_x + y) = 0 \end{array} $$ $$ \begin{array}{l} {\rm B}{\rm B}\bar {\rm A}p_x = {\rm B}{\rm B}\Big [(p_x + y)\dfrac{\partial }{\partial x} + (p_y- x)\dfrac{\partial }{\partial y}+\\ \qquad p_z \dfrac{\partial }{\partial z} + p_y \dfrac{\partial }{\partial p_x }-p_x \dfrac{\partial }{\partial p_y }\Big]p_x = \\ \qquad {\rm B}{\rm B}(p_y ) = {\rm B}\Big [-\Big (\dfrac{\partial }{\partial x}\dfrac{\partial }{\partial p_x } + \dfrac{\partial }{\partial y}\dfrac{\partial }{\partial p_y } + \dfrac{\partial }{\partial z}\dfrac{\partial }{\partial p_z }\Big )\Big] p_y = 0 \end{array} $$ $$ \bar{\rm A}{\rm B}{\rm B}x = \bar {\rm A}{\rm B}\Big [-\Big (\dfrac{\partial }{\partial x}\dfrac{\partial }{\partial p_x } + \dfrac{\partial }{\partial y}\dfrac{\partial }{\partial p_y } + \dfrac{\partial }{\partial z}\dfrac{\partial }{\partial p_z }\Big)\Big]x = 0 $$ $$ \begin{array}{l} \bar {\rm A}{\rm B}{\rm B}p_x = \bar {\rm A}{\rm B}\Big [-\Big (\dfrac{\partial }{\partial x}\dfrac{\partial }{\partial p_x } + \dfrac{\partial }{\partial y}\dfrac{\partial }{\partial p_y } + \dfrac{\partial }{\partial z}\dfrac{\partial }{\partial p_z }\Big)\Big]p_x =\\ \qquad \bar {\rm A}{\rm B}\Big [-\Big (\dfrac{\partial }{\partial x}\dfrac{\partial }{\partial p_x }p_x \Big) \Big] = 0 \end{array} $$ $$ {\rm B}\bar {\rm A}{\rm B}x = {\rm B}\bar {\rm A}\Big [-\Big (\dfrac{\partial }{\partial x}\dfrac{\partial }{\partial p_x } + \dfrac{\partial }{\partial y}\dfrac{\partial }{\partial p_y } + \dfrac{\partial }{\partial z}\dfrac{\partial }{\partial p_z }\Big)\Big]x = 0 $$ 同理可验证
$$ {\rm B}{\rm B}\bar {\rm A}y = {\rm B}{\rm B}\bar {\rm A}z = \bar {\rm A}{\rm B}{\rm B}p_x = \bar {\rm A}{\rm B}{\rm B}p_y = 0 $$ 附录B
$$ {\rm E} = 2(V_{xx} V_x + V_{xy} V_y + V_{xz} V_z )\dfrac{\partial }{\partial p_x } + \\ \qquad 2(V_{yx} V_x + V_{yy} V_y + V_{yz} V_z )\dfrac{\partial }{\partial p_y } + \\ \qquad 2(V_{zx} V_x + V_{xz} V_y + V_{zz} V_z )\dfrac{\partial }{\partial p_z } $$ $$ {\rm F} = \sum\limits_{i = 1}^3 {\sum\limits_{j = 1}^3 {\sum\limits_{k = 1}^3 {V_k V_{jk} V_{ij} \dfrac{\partial }{\partial p_i }} } } = \\ \qquad \sum\limits_{j = 1}^3 {\sum\limits_{k = 1}^3 {(V_k V_{jk} V_{xj} \dfrac{\partial }{\partial p_x } + V_k V_{jk} V_{yj} \dfrac{\partial }{\partial p_y } + V_k V_{jk} V_{zj} \dfrac{\partial }{\partial p_z })} } = \\ \qquad (V_x V_{xx} V_{xx} + V_x V_{yx} V_{xy} + V_x V_{zx} V_{xz} + V_y V_{xy} V_{xx} + \\ \qquad V_y V_{yy} V_{xy} + V_y V_{zy} V_{xz} + V_z V_{xz} V_{xx} + V_z V_{yz} V_{xy} + \\ \qquad V_z V_{zz} V_{xz} )\dfrac{\partial }{\partial p_x } + \\ \qquad (V_x V_{xx} V_{yx} + V_x V_{yx} V_{yy} + V_x V_{zx} V_{yz} + V_y V_{xy} V_{yx} + \\ \qquad V_y V_{yy} V_{yy} + V_y V_{zy} V_{yz} + V_z V_{xz} V_{yx} + V_z V_{yz} V_{yy} + \\ \qquad V_z V_{zz} V_{yz} )\dfrac{\partial }{\partial p_y } + (V_x V_{xx} V_{zx} + V_x V_{yx} V_{zy} + V_x V_{zx} V_{zz} + \\ \qquad V_y V_{xy} V_{zx} + V_y V_{yy} V_{zy} + V_y V_{zy} V_{zz} + V_z V_{xz} V_{zx} + \\ \qquad V_z V_{yz} V_{zy} + V_z V_{zz} V_{zz} )\dfrac{\partial }{\partial p_z } $$ -
![]()
图 2 三种算法积分两个周期得到的轨道
Figure 2. The orbits of three algorithms integral two orbit periods
![]()
图 3 三种算法积分3个周期得到的轨道
Figure 3. The orbits of three algorithms integral three orbit periods
![]()
图 4 三个算法积分100 000步的绝对能量误差
Figure 4. The absolute energy errors of 100 000 steps with three algorithms
![]()
图 6 三种算法积分100 000步的相对能量误差
Figure 6. The relative energy error of 100 000 steps with three algorithms
-
[1] Sezebehely V. Theory of Orbits:The Restricted Problem of Three Bodies. New York:Academic Press, 1967
[2] 刘林, 侯锡云.深空探测器轨道力学.北京:电子工业出版社, 2012 Liu Lin, Hou Xiyun. Orbital Mechanics of the Deep Space Probe. Beijing:Publishing House of Electronics Industry, 2012 (in Chinese)
[3] 雷汉伦. 平动点、不变流形及低能轨道. [博士论文]. 南京: 南京大学, 2015 Lei Hanlun. Equilibrium point, invariant manifold and low energy trajectory.[PhD Thesis]. Nanjing:Nanjing University, 2015 (in Chinese)
[4] Farquhar RW. The role of Sun-Earth collinear libration points in future space exploration//AAS Annual Meeting, 1999
[5] Dunham DW, Farquhar RW. Libration point missions, 1978-2002//Libration Point Orbits and Applications. Singapore:World Scientific, 2003
[6] Condon GL, Pearson DP. The role of humans in libration point missions with specific application to an Earth-Moon libration point gateway station//AAS 01-307, AAS/AIAA Astrodynamics Specialist Coference, 2001
[7] Lo MW. The inter-planetary superhighway and the origins program//IEEE Aerospace Conference Proceedings, 2002
[8] Baoyin HX, McInnes CR. Solar sail Halo orbits at the Sun-Earth artificial L1 point. Celestial Mechanics and Dynamical Astronomy, 2006, 94:155-17 doi: 10.1007/s10569-005-4626-3
[9] Xu M. Application of Hamiltonian structure-preserving control to formation flying on a J2-perturbed mean circular orbit. Celestial Mechanics and Dynamical Astronomy, 2012, 113(4):403-433 doi: 10.1007/s10569-012-9430-2
[10] Gómez G, Jorba A, Masdenmont J, et al. Study of the transfer from the Earth to a halo orbit around equilibrium point L1. Celestial Mechanics, 1993, 56:541-562 doi: 10.1007/BF00696185
[11] 袁建平, 孙冲, 方群.基于虚拟中心引力场方法的航天器转移轨道设计.力学学报, 2015, 47(1):180-184 doi: 10.6052/0459-1879-14-112 Yuan Jianping, Sun Chong, Fang Qun. Spacecraft's transfer orbit design based on the virtual central gravity?eld method. Chinese Journal of Theoretical and Applied Mechanics, 2015, 47(1):180-184 (in Chinese) doi: 10.6052/0459-1879-14-112
[12] 贾建华, 吕敬, 王琪.带脉冲的三维引力辅助变轨研究.力学学报, 2016, 48(2):437-446 doi: 10.6052/0459-1879-15-218 Jia Jianhua, Lü Jing, Wang Qi. Powered gravity assist in three dimensions. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(2):437-446 (in Chinese) doi: 10.6052/0459-1879-15-218
[13] Zhang RY, Luo JJ. Attainable sets approach for low-energy, lowthrust Interplanetary transfers//International Astronautical Congress. Beijing, China, 2013
[14] Lei HL, Xu B, Hou XY, et al. High-order solutions around triangular libration points in the elliptic circular restricted three-body problem and applications to low energy transfers. Communications in nonlinear science and Numerical Simulation, 2014, 19:3374-3398 doi: 10.1016/j.cnsns.2014.01.019
[15] Lei HL, Xu B. High-order analytical solutions around triangular libration points in the circular restricted three-body problem. Monthly Notices of the Royal Astronomical Society, 2013, 434(2):1376-1386 doi: 10.1093/mnras/stt1099
[16] Lei HL, Xu B, Hou XY, et al. High-order solutions of invariant manifolds associated with libration point in the elliptic circular restricted three-body system. Celestial Mechanics and Dynamical Astronomy, 2013, 117(4):349-384 doi: 10.1007/s10569-013-9515-6
[17] Lei HL, Xu B. Analytiacl study on the motions around equilibrium points of restricted Four-body problem. Compute Nonlinear Science Number Simulat, 2015, 29:441-458 doi: 10.1016/j.cnsns.2015.05.023
[18] Dellnitz M, Junge O, Post M. On Target for venus:Set oriented computation of energy efficient low thrust trajectories. Celestial Mechanics and Dynamical Astronomy, 2006, 95:357-370 doi: 10.1007/s10569-006-9008-y
[19] Howell K, Beckman M, Patterson C, et al. Representations of invariant manifolds for applications in three-body system//14th AAS/AIAA Space Flight Mechanics Conference, Maui, Hawaii, 2004
[20] Hsu CS. A theory of cell-to-cell mapping dynamical systems. Applied Mechanics, 1980, 147:931-939
[21] Ruth RD. A Canonical integration technique. IEEE Transactions on Nuclear Science, 1983, 30:2669-2671 doi: 10.1109/TNS.1983.4332919
[22] 李渊, 邓子辰, 叶学华等.基于辛理论的载流碳纳米管能带分析.力学学报, 2016, 48(1):135-139 doi: 10.6052/0459-1879-15-164 Li Yuan, Deng Zichen, Ye Xuehua, et al. Analysing the wave scattering in single-walled carbon nanotube conveying fluid based on the symplectic theory. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(1):135-139 (in Chinese) doi: 10.6052/0459-1879-15-164
[23] Pauli Pihajoki. Explicit methods in extended phase space for inseparable Hamiltonian problems. Celestial Mechanics and Dynamical Astronomy, 2015, 121:211-231 doi: 10.1007/s10569-014-9597-9
[24] Liu L, Wu X, Huang GQ, et al. Higher order explicit symmetric integrators for inseparable forms of coordinates and momenta MNRAS 459. Monthly Notices of the Royal Astronomical Society, 2016
[25] Hu WP, Deng ZC, Han SM, et al. Generalized multi-symplectic Integrators for a class of Hamiltonian nonlinear wave PDEs. Journal of Computational Physics, 2013, 235:394-406 doi: 10.1016/j.jcp.2012.10.032
[26] Mei LJ, Wu X, Liu FY. On preference of Yoshida construction over Forest-Ruth fourth-order symplectic algorithm. The Europe Physical Journal C, 2013, 73:2413 doi: 10.1140/epjc/s10052-013-2413-y
[27] 廖新浩, 刘林.辛算法在限制性三体问题数值研究中的应用.计算物理, 1995, 12(1):102-108 http://www.cnki.com.cn/Article/CJFDTOTAL-JSWL501.014.htm Liao Xinhao, Liu Lin. Applications of symplectic algorithms to the numerical researches of restricted three-body problem. Physic Computation, 1995, 12(1):102-108 (in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-JSWL501.014.htm
[28] Sun W, Wu X, Hang GQ. Symplectic integrators with potential derivatives to third order. Research in Astronomy and Astrophysics, 2011, 11:353-368 doi: 10.1088/1674-4527/11/3/009
[29] Hairer E, Wanner G. Solving Ordinary Differential Equation.[s.l.]:Springer Verlag, 1991
[30] 李言俊, 张科, 吕梅柏等.利用拉格朗日点的深空探测技术.西安:西北工业大学出版社, 2014 Li Yanjun, Zhang Ke, Lü Meibo. et al. Deep Space Exploration Using Lagrange Equilibrium Point. Xi'an:Northwestern Polytechnical Chivesity Press, 2014 (in Chinese)
[31] Bernelli-Zazzera F, Topputo F, Massari M. Assessment of Mission Design Including Utilization of Libration Points and Weak Stability Boundaries.[s.l.]:Politecnico di Milano, 2004
[32] Gómez G, Jorba A, Masdemont J, et al. Study refinement of semianalytical halo orbit theory. Final Report, ESOC Contract No:8625/89/D/MD(SC), 1991
[33] Omelyan IP, Mryglod IM, Folk R. Optimized Forest-Ruth-and Suzuki-like algorithms for integration of motion in many-body systems. Computer Physics Communications, 2002, 146:188-202 doi: 10.1016/S0010-4655(02)00451-4
[34] Omelyan IP, Mryglod IM, Folk R. Symplectic analytically integrable decomposition algorithms:Classification, derivation, and application to molecular dynamics, quantum and celestial mechanics simulations. Computer Physics Communications, 2003, 151:272-314 doi: 10.1016/S0010-4655(02)00754-3
-
期刊类型引用(3)
1. 邱志平，姜南. 随机和区间非齐次线性哈密顿系统的比较研究及其应用. 力学学报. 2020(01): 60-72 . 
本站查看
2. 尹婷婷，邓子辰，胡伟鹏，李庆军，曹珊珊. 空间刚性杆-弹簧组合结构轨道、姿态耦合动力学分析. 力学学报. 2018(01): 87-98 . 本站查看
3. 袁国强，李颖晖. 二维稳定流形的自适应推进算法. 力学学报. 2018(02): 405-414 . 本站查看
其他类型引用(0)
下载:
计量
- 文章访问数: 988
- HTML全文浏览量: 129
- PDF下载量: 678
- 被引次数: 3
