應用三階力梯度辛算法求解純Kepler二體問題

李 榮

(安陽師范學院數學與統計學院，河南安陽 455000）

辛幾何最早可以追溯到19世紀，英國天文學家哈密頓為了研究牛頓力學，引進了廣義坐標和廣義動量來表示系統的能量，稱為哈密頓函數。辛幾何應用于數值方法，最早出現于20世紀80年代初的文獻。我國數學家馮康院士［1－2］于1984年首次提出了求解哈密頓系統的隱式辛算法，他指出:“傳統的數值方法，除了極少數的例外，都不可避免的帶有人為耗散性等歪曲體系特征的缺陷”。而辛算法彌補了這一缺陷，是一種具有高保真性算法。隱式辛算法主要運用正則變換的隱式生成函數從事辛算法的研究，優勢在于適用范圍廣，它適用于幾乎任何形式的哈密頓函數，但由于在計算中要使用迭代方法，而使得它的計算效率大大降低。同一時期，Ruth［3－4］建立了求解哈密頓系統需要分解為動能T部分和勢能V部分的顯式辛算法。考慮到顯式辛算法的計算效率比隱式辛算法高，當哈密頓系統能分解時，顯式辛算法是更好的選擇。

力梯度顯式辛算法［3］是一種獨特的辛算法，它與非力梯度顯式辛算法同時出現。通過力實現力梯度的計算，力梯度顯式辛算法可以做到只含有正積分步長，這滿足了在求解時間不可逆的方程時對正積分步長的需求［5－6］。Chin 和合作者［7－9］給出了滿足算子形式對稱的力梯度算法并把它擴展到四階。本文構造了一系列算子形式不對稱的三階力梯度辛算法，并把它應用到求解純Kepler二體問題［10－12］的軌道長半徑中去，表明其具有很高的誤差精度。

1 三階力梯度顯辛算法的構造

下表中，我們列舉出幾種滿足條件的辛算法。

表1.1算子型三階力梯度顯辛算法的時間系數表

表1.1算子型三階力梯度顯辛算法的時間系數表

算法 a1 a2 b1 b2 b3 k1 23 13 14 3 4 0 1 482 59 49 15 27 40 18 7 4803 7 15 8 15 17 75 112 3 16 19 1344

表1.2 算子型三階力梯度顯辛算法得時間系數表

表1.2 算子型三階力梯度顯辛算法得時間系數表

算法 a1 a2 a3 b1 b2 k1 14 7 12 16 47 37 1 482 15 26 45 29 25 52 27 52 7 6243 17 62 105 4 15 49 124 75 124 19 1488

2 數值模擬

Kepler二體問題［5－6］是天體動力學中最基本的模型之一，是研究天體精確運動的理論基礎，具有重要的研究意義。在零階近似下，天體對外的引力場可被看做質量集中于質心的質點所產生的引力場，所以兩天體的運動問題就可以近似視為兩質點在其相互之間的萬有引力作用下的動力學問題。哈密頓函數為

圖2.1中，選取一條坐標和動量的初始條件分別選擇 1，0，0，1的運動軌道，步長為周期的1/20。通過觀察可知，每一種算法都能夠很好的保持了辛算法的優勢，即軌道長半徑隨時間周期型變化。同時，每一種算法的誤差精度都能夠達到10－7左右，表明這些三階力梯度辛算法非常適用于計算天體力學中軌道參數問題。由于軌道長半徑的相對誤差和能量的相對誤差具有正比例關系，這是因為Kepler能量E和軌道長半徑a具有如下關系E－μ/2a，此處μ為引力常量。所以圖2.1也可以反映系統的能量誤差精度。2.1(a）、(b）、(c）圖對應算子型力梯度辛算法，(d）、(e）、(f）圖對應算子型力梯度辛算法。總體來說，算子型力梯度辛算法的誤差精度相較于算子型要高一個數量級左右，其中以算法最優。

3 小結

在研究哈密頓系統的動力學性質時，傳統的數值方法由于存在人為耗散性等歪曲體系特征導致系統能量得不到保持，這將使得數值解法產生失真。本文構造了一系列算子形式不對稱的三階力梯度辛算法，并通過求解純Kepler二體問題的軌道長半徑來驗證了其有效性。

［1］Feng K.Beijing Symposium on Differential Geometry and Differential Equations［M］.Beijing:Science Press，1985.

［2］馮康，秦孟兆.Hamilton動力系統的 Hamilton算法［J］.自然科學進展 －國家重點實驗室通訊［J］，1991，(2）:102 ～112.

［3］Ruth R D.A canonical integration technique［J］.IEEE Transactions on Nuclear Science，1983，NE － 30(4）:2669～2671.

［4］Forest E，Ruth R D.Fourth －order symplectic integration［J］.Phys.D，1990，43:105 ～117.

［5］Rong Li，Xin Wu.A symmetric product of two optimal third－order force gradient symplectic algorithms［J］.ACTA PHYSICA SINICA，2010，59(10），7135 ～7143.

［6］Rong Li，Xin Wu.Optimized third－order force－gradient symplectic algorithms［J］.SCIENCE CHINA，2010，53(9），1600～1609.

［7］Chin S A.Symplectic integrators from composite operator factorizations［J］.Phys.Lett.A，1997，226:344 ～348.

［8］Chin S A，Chen C R.Forward symplectic integrators for solving gravltational fewbody problems［J］.Celestial Mechanics and Dynamical Astronomy，2005，91:301 ～322.

［9］Chin S A.Physics of symplectic integrators:Perihelion advances and symplectic corrector algorithms［J］.Physical Review E，2007，75:036701

［10］劉福窯，伍歆，陸本魁.幾類辛方法的數值穩定性研究［J］.天文學報，2006，47(4）:418－431.

［11］徐佳.含力梯度顯辛算法拓廣于攝動二體問題［J］.江西科學，2009，27(4）:506－509.

［12］Jia xu and Xin Wu.“Several fourth－order force gradient symplectic algorihhms”.Research in Astronomy and Astrophysics 2010，10(2）:173.