比较单步形式和两步格式的速度双因子方法

稳定且高精度的数值算法是天体力学和动力天文研究的基本工具． 传统的常微分方程数值算法如经典的四阶龙格库塔法局部截断误差大，长期积分会导致误差线性积累，不利于运动守恒量的保持． Nacozy提出的流形改正原理是在低阶传统数值算法基础上加入稳定项，稳定化处理后误差得到有效控制，不变积分( 守恒量) 或拟积分的数值精度和稳定性都得到保证，从而为天体长期定性研究提供基础． 流形改正方法目前已得到长足发展． 　　 　　1 速度因子方法 　　 　　日心坐标系下受摄二体问题运动微分方程一般形式为: dv/dt=－r/r3+a．其开普勒能量、角动量和拉普拉斯积分表达式为: K=v2/2－ / r，L = xv，P = vL－r / r．由于数值算法存在误差，使得数值解计算得到的积分量偏离了积分原始超曲面，故不再是守恒量． 流形改正方法采用最小二乘法原理将偏离部分逐步拉回到原始位置．速度因子方法是流形改正原理的直接拓展，即稳定项仅施加在速度分量部分． 与其他方法相比，在计算精度相同的前提下，该方法需要更少的 CPU 时间耗费，且形式更为简单． 速度因子方法有严格和近似之分，严格的速度因子方法要求改正因子必须严格满足积分关系，而近似方法要求并不严格． 　　 　　1． 1 单步方法。

从积分定义上看，严格的开普勒能量关系应满足: K*= v* 2/2－ / r．为挽回积分超曲面，可定义速度改正解和数值解之间满足因子关系: v*= sv*．容易求得因子的具体表达式为:【1】 　　 　　该方法即为单步改正的速度单因子方法．1． 2 两步方法。

两步方法根据线性转换关系得到，Fukushima 认为坐标和速度之间满足线性关系: v*= ( v－x) ．由于引入了两个不同因子，所以因子表达式不能在一步内获得，必须借助其他守恒量．第 1 步: 拉普拉斯积分的前半部分为 F=v*L，且严格满足 Fv*= 0． 由此可得到因子 : = ( Fv) /( Fx) ．第 2 步: 因子 可从开普勒能量关系中求得:【2】 　　 　　相比式( 1) ，( 2) 稍显复杂，但正是由于线性关系的存在，使得速度部分在改正过程中同时保持了两个运动积分常数． 　　 　　2 数值模拟 　　 　　为比较两种方法的数值稳定性和精度，本文以日－木－土三体问题为基本研究模型，五阶龙格库塔方法( RK5) 为基本数值工具，固定时间步长为 36． 525 天，相关初始值从 JPL 星历表 DE405 查得． 由于三体问题没有分析解，以高阶变步长方法计算结果作为参考值． 为书写简单，M1 代表速度双因子方法( 两步法) ，M2代表速度单因子方法( 单步法) ． 数值结果如图 1 所示，两者在提高轨道半长轴精度方面，木星大约提高了 3个量级，土星稍小，为 2 个量级． 对偏心率，M1 在积分达到 1000 年后基本能提高 2 个量级，但 M2 基本没什么改变． 由于偏心率与拉普拉斯积分密切相关，足以说明 M1 改正了除能量积分之外的第二个积分，该现象在近点角距中也得到体现． 需要强调的是，两种速度因子改正方法对木星的改正均比较明显，土星稍微差一些，原因在于数值积分步长与木星的轨道周期更为接近． 在相对位置误差的改进方面，如图 2 所示．两者差别不大，在积分达到 10000 年时，基本都提高了两个精度值． 说明只对速度项改正亦能起到对坐标部分改正的作用． 这正是速度因子方法的一个优势，形式上更为简单，最终效果却一样． 单从构造上看，速度因子方法对坐标部分不起作用，但只限于积分第一步，从第二步起，一旦从速度部分将积分拉回到初始轨道，相应的坐标数值解也得到改正．【图1—2.略】 　　 　　3 结论。

基于流形改正原理，速度因子改正方法有效的保持了个别天体的孤立积分． 本文比较了两种严格的速度因子改正方法，一种是单步改正方法，另一种是基于线性转换格式的两步法． 数值模拟发现两者均能有效保持与孤立积分相对应的轨道根数精度，尤其是两步法，显着提高了木星轨道偏心率和近点角距的数值精度．但是由于积分步长恒定，且积分参考值精度不高，总体改进效果仍然不明显，以后应关注在多体问题的应用．。