能量稳定化方法在平面圆形限制性三体问题中的应用

赵杨华，马大柱

(湖北民族学院 理学院，湖北 恩施 445000)

复杂动力系统的非线性行为研究依赖于数值计算方法的稳定性和精度.传统低阶数值算法在一定程度上能反映动力系统的相空间结构，但是对于某些系统的混沌行为，由于其对初始条件极为敏感，并且对数值结果的精度要求非常高，所以寻找合适的数值计算方法是当前非线性研究工作的基本问题.为避免数值误差产生的假混沌现象，目前文献中主要采取以下三种方式解决：一是采用高阶算法求解微分方程(组)；二是采用辛算法[1]；三是对传统低阶数值算法加以改造，即后稳定化方法[2].研究表明，高阶算法由于截断误差小，计算耗时较长，不利于广泛应用.辛算法虽然能保持系统的辛结构，但是其主要用于处理可分的哈密顿系统，所以应用有限.相对而言，后稳定化方法不仅节省了计算时间，而且能提高数值积分的稳定性和精度.

1 能量稳定化方法

对一含有不变运动积分的微分动力系统：φ(x)=c，c为积分常数，理论情况下应满足：

ε(x)=φ(x)-c=0

(1)

但是数值误差的存在，实际计算导致ε≠0.原因在于数值算法存在截断误差以及计算机将舍入误差代入到各种算法中.为解决该问题，Nacozy[3]在1971年提出了流形改正方案，即后稳定化方法，以不变积分常数作为控制量，将数值算法得到的超曲面用最小二乘法的原理以最短路径拉回到原始的超曲面上.

简单而言，假设x代表数值方法得到的坐标和速度矢量，x′为改正后的坐标和速度矢量，两者之间应该满足x′=x+Δη，Δη即为改正量，Δη必须满足最小二乘法原理，并且是以最短的距离，所以Δη=-ET(EET)-1ε，式中上标T代表转置矩阵，其中E满足：

(2)

多次研究表明，该方案非常适合用于稳定常微分方程模型.将后稳定化方法推广到含有能量积分模型中即为能量稳定化方法[4]，尤其是只存在单个能量积分情况下，能量稳定化可以简写为标度因子的形式，即：

x′=sx

(3)

s的形式由具体运动的哈密顿形式给出.本文基于以上结论，进一步推广该方案，将雅克比积分常数看做不变积分应用能量稳定化方法讨论平面圆形限制性三体问题.

2 平面圆形限制性三体问题模型

图1 雅克比积分随时间的变化关系

平面圆形限制性三体问题模型[5]是三体其中一体的质量远小于其它两体的特殊三体问题模型.可以不考虑质量最小体对其它天体的引力作用，于是两主天体在相互引力的作用下作圆形运动，对于两主天体引力下的第三体则作限制性运动，系统至多为3自由度，故不可积，但是仍然可以找到一个稳定的积分.当主天体作圆运动时，第三体的运动方程可记为：

(4)

(5)

(6)

该常数是旋转坐标系中的能量，将其等同于能量积分，应用能量稳定化方法，可以得到如下形式的改正因子：

(7)

3 结果分析

如前所述，由于模型不可积，只能用数值方法求解.本文以四阶龙格库塔方法(RK4)为基本的数值算法，根据程序设计原则，能量稳定化程序应加在基本算法之后.为考察稳定化程序的作用，首先观察雅克比积分随时间的变化规律，如图1所示.从理论上分析雅克比积分为常数，在积分运动方程时应该保持不变，但是RK4所得到的误差结果却是直线增加，也就是说误差随时间逐渐积累，即得不到有效控制，由此以来积分时间越长，数值稳定性越差.从计算精度上看，RK4方法只能获得不超过10-6的结果，与理想值有一定差距.然而加入稳定化程序之后，结论大不相同，无论是稳定性还是计算精度都获得了令人满意的结果.容易得出，能量稳定化方法在圆形限制性三体问题中也有重要的应用价值.

图2即为庞加莱截面图，此时将RK4方法和稳定化程序作对比，可以很清晰的看出传统低阶数值方法获得的相空间结构并不真实，截面显示该轨道是非周期的，也就是说该条轨道可能是混沌的，该结果与实际情况相差甚大，加入稳定化可以得到真实情况，轨道是有三个岛链组成的拟周期轨道，没有混沌行为，用RK4得到的混沌现象只是假象.

图2 庞加莱截面图

为验证上述结论是否正确，本文还用高阶算法做出功率谱和偏差轨道的结构演化图，相邻轨道的坐标初始值改变量为0.001，如图3和图4所示.两种动力学分析方法均表明该轨道是周期的.

图3 功率谱图 图4 偏差轨道演化图

4 结论

基于平面圆形限制性三体问题模型存在孤立的不变积分常数即雅克比积分，本文将n体问题中广泛应用的能量稳定化方法应用到该模型，数值实验得出能量稳定化方法对雅克比积分保持了较高的数值稳定性和精度，而且该方法能有效的避免数值误差产生的假混沌现象，与高阶精度算法得出结论一致，从而提供了更为简单有效的讨论高维动力系统非线性行为的工具.

[1] Feng K,Sym.Diff.Geometry & Diff. Equation[M].Beijing:Science press,1985.

[2] Wu X, Huang T Y, Wan X S,et al.Comparison among correction methods of individual Kepler energies in n-body simulations[J].A J,2007,133(6):2643-2653.

[3] Ma D Z, Wu X, Zhu J F.Velocity scaling method to correct individual Kepler energies[J].New Astron,2008,13(5):216-223.

[4] Nacozy P E. The use of Integrals in Numerical Integrations of the n-body problem[J].Ap&SS,1971,14(11):40-51.

[5] Simo C,Stuchi T J.Central stable/unstable manifolds and the destruction of KAM tori in the planar hill problem[D].Physica D,2000,140(6):1-32.