穩定且高精度的數值算法是天體力學和動力天文研究的基本工具． 傳統的常微分方程數值算法如經典的四階龍格庫塔法局部截斷誤差大，長期積分會導致誤差線性積累，不利于運動守恒量的保持． Nacozy提出的流形改正原理是在低階傳統數值算法基礎上加入穩定項，穩定化處理后誤差得到有效控制，不變積分\\( 守恒量\\) 或擬積分的數值精度和穩定性都得到保證，從而為天體長期定性研究提供基礎． 流形改正方法目前已得到長足發展．

1 速度因子方法

日心坐標系下受攝二體問題運動微分方程一般形式為: dv/dt=－μr/r3+a．其開普勒能量、角動量和拉普拉斯積分表達式為: K=v2/2－μ / r，L = x×v，P = v×L－μr / r．由于數值算法存在誤差，使得數值解計算得到的積分量偏離了積分原始超曲面，故不再是守恒量． 流形改正方法采用最小二乘法原理將偏離部分逐步拉回到原始位置．速度因子方法是流形改正原理的直接拓展，即穩定項僅施加在速度分量部分． 與其他方法相比，在計算精度相同的前提下，該方法需要更少的 CPU 時間耗費，且形式更為簡單． 速度因子方法有嚴格和近似之分，嚴格的速度因子方法要求改正因子必須嚴格滿足積分關系，而近似方法要求并不嚴格．

1． 1 單步方法

從積分定義上看，嚴格的開普勒能量關系應滿足: K*= v* 2/2－μ / r．為挽回積分超曲面，可定義速度改正解和數值解之間滿足因子關系: v*= sv*．容易求得因子的具體表達式為:【1】


該方法即為單步改正的速度單因子方法．1． 2 兩步方法。。

兩步方法根據線性轉換關系得到，Fukushima 認為坐標和速度之間滿足線性關系: v*= α\\( v－βx\\) ．由于引入了兩個不同因子，所以因子表達式不能在一步內獲得，必須借助其他守恒量．第 1 步: 拉普拉斯積分的前半部分為 F=v*×L，且嚴格滿足 F·v*= 0． 由此可得到因子 β:β = \\( F·v\\) /\\( F·x\\) ．第 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 結論

基于流形改正原理，速度因子改正方法有效的保持了個別天體的孤立積分． 本文比較了兩種嚴格的速度因子改正方法，一種是單步改正方法，另一種是基于線性轉換格式的兩步法． 數值模擬發現兩者均能有效保持與孤立積分相對應的軌道根數精度，尤其是兩步法，顯著提高了木星軌道偏心率和近點角距的數值精度．但是由于積分步長恒定，且積分參考值精度不高，總體改進效果仍然不明顯，以后應關注在多體問題的應用．