限制性三體問題中顯式辛格式的構造

韋炳威，李銀山

（河北工業大學 機械工程學院，天津 300130）

限制性三體問題中顯式辛格式的構造

韋炳威，李銀山

（河北工業大學 機械工程學院，天津 300130）

針對圓型限制性三體問題（CR3BP）研究了顯式辛積分格式的構造問題.首先通過把CR3BP對應的哈密頓函數拆分成若干個二階冪零哈密頓系統，得到每個二階冪零哈密頓系統對應的顯式歐拉格式.然后證明了每個顯式歐拉格式都是自共軛算子，針對這樣的特點提出了一種由這些歐拉格式復合得到的顯式組合辛格式.最后利用本文提出的顯式辛格式與其他積分器求解了CR3BP下的一般軌道和Halo軌道，驗證了顯式辛格式有效性和優越性.研究發現顯式辛格式能長時間保持系統能量誤差在一定范圍內波動不會出現發散，且計算精度高于同階的非辛算法.

辛幾何算法；限制性三體問題；哈密頓系統；組合格式

0 引言

哈密頓系統具有2類守恒律[1]：一類是Liouville-Poincaré守恒律，即相空間內偶數維體積不變；另一類是能量、動量、角動量的運動不變量的守恒律.對于一個給定的數值計算方法，人們自然希望該算法能夠保持以上2個守恒律.但是研究表明，現有的數值計算方法大多不能保持這2個守恒律.以龍格-庫塔（Runge-Kutta，RK）方法為例，對于不可分的哈密頓系統，顯式的RK方法不可能保持Liouville-Poincaré守恒律.辛幾何算法是一種能夠自動保持Liouville-Poincaré守恒律，對于運動不變量的逼近與算法本身的逼近階相當的幾何積分算法.于是辛幾何算法成為了一個求解哈密頓系統長期演化的最佳方法.

圓型限制性三體問題（Circular Restricted Three-body Problem，CR3BP）是一個特殊的自治哈密頓系統.對于該系統的理論研究已經趨于成熟，許多動力學特性已經給出[11-20].現在一般用于研究該系統的積分器都以RK方法為主.例如美國噴氣推進實驗室所使用的主要積分器之一是變步長的RKF7（8）[2]，它實質是7階13級的RK方法.RKF7（8）擁有較高精度被廣泛使用于天體動力學和任務軌道設計，但是RK方法是非辛算法，系統的能量不能在該算法下被長期跟蹤.通過計算發現，無論使用4階精度的RK方法（RK4）還是使用7階的RKF7（8）方法，系統的能量都會隨著時間的推移而不斷地積累著非辛算法帶來的人為耗散而產生的誤差.若希望跟蹤一條初值敏感度很高的混沌軌道，那么使用傳統RK方法勢必會造成不準確的結果.

要使辛幾何算法在限制性多體問題中得到廣泛應用，就應該找到一種在效率上不亞于傳統RK方法的顯式辛格式.就目前而言，在這方面的研究還甚少.陳云龍[3]從Lie算子出發，將力梯度辛方法應用于CR3BP，構造了一種顯式的辛方法.Su X N[4]運用對數哈密頓算法克服了定步長辛算法在大引力梯度區域的失真問題.然而對于一類比較特殊的辛可分哈密頓系統而言，可以通過組合多個二階冪零哈密頓系統的歐拉顯式格式得到任意偶數階的顯式辛格式[5-7].CR3BP的哈密頓函數經過一定變形之后，可看出其恰巧是辛可分系統，完全可以用上述方法去構造顯式的辛格式.

1 辛可分哈密頓系統

1.1 基本理論[6-7]

哈密頓函數H（p，q）是辛可分的，如果

式中φi（u）是n個變量的標量函數，且上式滿足條件

其中Hi（p，q）稱為二階冪零哈密頓系統.對于每一個二階冪零系統Hi（p，q）都可以構造1個如下的歐拉顯式格式Eτi：式中：In是n階單位矩陣；是哈密頓系統Hi（p，q）的精確相流，故是辛的.以下的組合格式是原辛可分系統H（p，q）的一階顯式辛格式

文獻 [6]中介紹了通過Lie級數得到用低階格式組合成高階格式的方法，并證明了以下結論：若S（τ）是自共軛積分算子，其階數為2n，則當C1、C2滿足

1.2 圓型限制性三體問題[14]

考慮1個質量可忽略的質點P在2個有限質量天體P1、P2引力作用下的運動情況，其中P1的質量為m1，P2的質量為m2，且有m1＞m2.假設P1、P2都相對于它們的質心做平面圓周運動，質點P能夠自由地在P1、P2附近的空間內運動，且不會影響P1、P2的運動.為了方便地描述該系統，取計算單位如下：1個質量單位為m1＋m2，1個長度單位為P1、P2之間的距離；選取1個時間單位使得P1、P2繞它們質心旋轉的周期為2π.在以上無量綱化后，引力常數G=1.P2的質量為μ=m2/（m1＋m2），P1的質量為1-μ.在以P1、是2n＋2階的.P2質心為原點的旋轉坐標系下，系統的哈密頓函數為

其中有效勢能U（q）為

在計算共線拉格朗日點附近的軌道時往往將坐標系的原點平移到相應的拉格朗日點.本文以地-月系統的拉格朗日2（L2）點為例進行計算研究，故將上述坐標系的原點平移到L2點，再取L2點到月球的距離γ為1個單位，哈密頓函數改寫為

現將哈密頓函數（10）作如下的分解：

顯然T（p）和U（q）是二階冪零系統，而H2（p，q）可以作如下處理

這樣可以清楚看出，實際上H2（p，q）可分解為4個二階冪零系統.從而CR3BP實際上是一個辛可分的哈密頓系統，于是理論上可以找到顯式的辛格式.

2 顯式辛格式的構造

2.1 共軛積分算子

由上述結論，CR3BP的哈密頓函數可以分解成6個二階冪零哈密頓系統

根據（5）式，可以通過復合6個二階冪零系統的歐拉格式得到1階顯式格式

若算子SH（τ）可逆，則有，即.于是SH（τ）的共軛算子可以通過求逆得到

其中因為φi是2次的n元函數，Dφi是一個線性算子，可以用一個矩陣Mi表示.考慮到條件（2），有

2.2 高階組合格式

3 實例軌道計算

3.1 一般軌道計算

以地-月系統為例，計算地、月附近的一般軌道，這些軌道對初值敏感程度不高，方便比較各個算法對系統能量的長期跟蹤能力.對于地-月系統μ=0.012 15，γ=0.167 829 913 1，基于哈密頓函數（10），選擇初值

本文采用上述提出的4階顯式組合辛格式（Comp4），與RKF7（8）、RK4和4階Gauss-Legendre隱式辛格式[1]（GL4）分別計算以上初值的軌道，取時間步長為π/1 000，時間總步數11 000步，并計算每個時間步對應的系統能量Ht與初始系統能量H0的差值.由于系統能量誤差很小，故對其取對數ln（|Ht-H0|）進行分析.

圖1 一般軌道的系統能量誤差曲線（RKF7（8））Fig.1 Curve of system energy error of regular orbit（RKF7（8））

圖2 一般軌道的系統能量誤差曲線（RK4）Fig.2 Curve of system energy error of regular orbit（RK4）

圖1～圖4所示系統能量誤差隨時間步變化曲線，CR3BP是自治保守系統，理論上系統能量守恒，但因數值積分存在截斷誤差，所以計算結果顯示系統能量隨著時間步的增加都有所變化.其中圖1由于RKF7（8）具有7階精度，系統的能量變化相比較于其他3種積分器都很小.但傳統的RK法存在人工耗散，系統能量誤差在不斷地累積，并且出現了明顯的發散.圖2是傳統4階RK4積分器下的計算結果，由于只有4階精度，誤差相比于RKF7（8）要大，而且同樣存在能量誤差的發散.圖3是4階Gauss-Legendre隱式辛格式得到的曲線，圖4是本文提出的4階顯式組合辛格式，可以發現它們的能量誤差曲線在一定范圍內波動，并沒有發生發散的情況.雖然辛算法因截斷誤差在長時間的積分中也可引起能量誤差的累積，但能量保持性已經遠遠優于非辛算法.另外隱式辛格式GL4與本文提出的顯式組合辛格式Comp4比較，在精度上稍微優于Comp4，這是隱式格式自身的優點所決定的.但不能忽視Comp4的計算成本遠遠小于隱式辛格式GL4這一特點，Comp4在這次計算中總共進行了66 000次函數的計算，而GL4在計算中總共進行了325 208次函數的計算，這顯示出顯式組合辛格式的優越性.

圖3 一般軌道的系統能量誤差曲線（GL4）Fig.3 Curve of system energy error of regular orbit（GL4）

圖4 一般軌道的系統能量誤差曲線（Comp4）Fig.4 Curve of system energy error of regular orbit（Comp4）

圖5所示地-月系統一般軌道的軌跡，坐標原點平移回地-日質心，距離單位為千米.由于軌道對初值不敏感，4種積分器得到的軌道基本一致.因為RKF7（8）的計算結果在整體精度上都要高于其他3種積分器，所以不妨以RKF7（8）得到的軌道視為1條參考軌道，而其他積分器得到的軌道與該軌道相比得到1個參考誤差|Δq|=|q（t）-qf（t）|，其中qf（t）是參考軌道的廣義坐標向量.如圖6所示RK4、Comp4、GL4與參考軌道比較的軌道參考誤差隨時間的變化.RK4得到的軌道在該時間內最大誤差達到140 km，而Comp4得到的軌道最大誤差不到80 km，精度最好的GL4得到的軌道最大誤差只有18 km.這提供了另一個角度說明了顯式組合辛格式和Gauss-Legendre格式這樣的辛格式跟傳統RK方法相比的優越性.

圖5 地-月系統一般軌道Fig.5 Regular orbit of Earth-Moon system

圖6 軌道參考誤差曲線Fig.6 Curves of reference error of orbits

3.2 Halo軌道計算

采用RKF7（8）、RK4、GL4和Comp4分別計算地-月系統中圍繞L2點，旋轉坐標系下x方向振幅為15 700 km的Halo軌道.由于Halo軌道是圍繞拉格朗日點一簇不穩定的周期軌道，所以對初值的敏感程度很高.極小的偏差也將會導致Halo周期軌道的破壞，所以長期保持Halo周期軌道的能力也可反映出積分器的精度好壞.

圖7 Halo軌道的系統能量誤差曲線（RKF7（8））Fig.7 Curve of system energy error of Halo orbit（RKF7（8））

圖8 Halo軌道的系統能量誤差曲線（RK4）Fig.8 Curve of system energy error of Halo orbit（RK4）

如圖7至圖10所示RKF7（8）、RK4、GL4和Comp4計算Halo軌道得到的系統能量誤差隨時間變化曲線，步長為1/2 000倍的Halo軌道周期，積分時間為4.5倍的Halo軌道周期；圖11所示與之對應的Halo軌道在x-y平面的投影，坐標原點建立在L2點.結果顯示，由于RKF7（8）具有很高的精度，在整個積分時間內都能保持很小的系統能量誤差，得到的Halo軌道也保持得最好，直到接近4.5倍周期時才逐漸稍微偏離周期軌道.RK4得到的Halo軌道在不到4倍周期時就開始偏離了周期軌道，系統能量的誤差在隨著時間的推移不斷累積，并且在接近4.5倍周期時出現跳躍現象.這是由于軌道偏離周期軌之后靠近了月球，造成引力梯度突增所至.因為采用本文提出的GL4只有4階精度，得到的Halo軌道也在不到4倍周期就偏離了周期軌道，但是系統能量誤差卻相比RK4得到的結果要穩定許多，沒有能量誤差的積累，保持在較低水平波動.GL4得到的系統能量誤差曲線在接近4.5倍周期的時候也出現了明顯的跳躍，從圖11可知同樣是因為軌道過于接近月球，使得引力梯度突增而造成的結果.若要在該區域保持較高的精度，則需要控制步長的方法，可以在引力梯度增高的區域縮小時間步長.Comp4所得到的Halo軌道也在不到4倍周期時出現了偏移，但由于偏移方向恰好朝著背離月球的方向，故系統的能量誤差一直保持在較低值的范圍振蕩，沒有出現跳躍.綜上所述，對于初值敏感的Halo軌道，4種積分器得到的結果完全不同.RKF7（8）因其自身的高精度性在該積分時段內保持住了Halo軌道，RK4、GL4、Comp4由于只有4階精度，都在不到4倍周期的時間內出現了Halo軌道的偏移，并且偏移方向各不相同.從系統能量誤差曲線看出RK4在每個周期都在累積著人工耗散帶來的誤差，GL4、Comp4都沒有人工耗散帶來的誤差問題.

圖9 Halo軌道的系統能量誤差曲線（GL4）Fig.9 Curve of system energy error of Halo orbit（GL4）

圖10 Halo軌道的系統能量誤差曲線（Comp4）Fig.10 Curve of system energy error of Halo orbit（Comp4）

圖11 4種積分器得到的Halo軌道Fig.11 Halo orbits obtained by four integrators

4 結論

CR3BP實質上是辛可分的哈密頓系統，理論上可以建立顯式辛格式.首先通過組合CR3BP的各個二階冪零哈密頓系統對應的精確歐拉相流，得到1階辛格式.為了得到2階格式，需要求出1階格式對應的共軛積分算子.通過證明各個歐拉相流算子都是自共軛算子，得到1階格式共軛積分算子的一般形式，進而得到了2階顯式辛格式.最后根據秦孟兆在文獻 [6]中的結果，可以構造出CR3BP任意偶數階的顯式辛格式.

通過一般軌道和Halo軌道的實例計算結果，表明了組合方法得到的4階顯式辛格式相比傳統RK4法在保持能量不發散方面，有明顯的優越性.同時相比GL4格式，結果精度稍欠缺，但計算成本大大減少，這說明顯示組合格式在精度方面同隱式辛格式比較還有改進的空間，希望后期對這方面進行完善.

本文提出的顯式辛格式的構造方法不僅僅適用于CR3BP，也同樣適用于雙圓型限制性四體問題，原因是本質上該問題也是辛可分哈密頓系統.因為該系統不再是自治系統，故不存在能量守恒，但可以通過擴充原有的哈密頓系統維數使其變成自治系統，從而構造形式上的守恒“能量”，余下的分析方法與本文討論的一致.

[1]馮康，秦孟兆.哈密頓系統的辛幾何算法[M].杭州：浙江科學技術出版社，2003.

[2]Erwin Fehlberg.Classical fifth-，sixth-，seventh-，and eighth-order runge-kutta formulas with stepsize control[R].NASA Technical Report TR R-287，Huntsville，Alabama：Marshall Space Flight Center，1968.

[3]陳云龍，伍歆.力梯度辛方法在圓型限制性三體問題中的應用[J].物理學報，2013，62（14）：140501-140501.

[4]Su X N，Wu X，Liu F Y.Application of the logarithmic Hamiltonian algorithm to the circular restricted three-body problem with some post-Newtonian terms[J].Astrophysics and Space Science，2016，361（1）：1-12.

[5]Yoshida H.Construction of higher order symplectic integrators[J].Physics Letters A，1990，150（s5-7）：262-268.

[6]Qin M Z，Zhu W J.Construction of higher order symplectic schemes by composition[J].Computing，1992，47（47）：309-321.

[7]Feng K，Wang D L.Variations on a theme by Euler[J].Journal of Computational Mathematics，1998，16（2）：97-106.

[8]Feng K，Wu H M，Qing M Z，et al.Construction of canonical difference schemes for Hamiltonian formalism via generating functions[J].Journal of Computational Mathematics，1989，7（1）：71-96.

[9]Ni X T，Wu X.New adaptive time step symplectic integrator：an application to the elliptic restricted three-body problem [J].Research in Astronomy& Astrophysics，2014，14（10）：1329-1342.

[10]Lu W T，Zhang H，Wang S J.Application of symplectic algebraic dynamics algorithm to circular restricted three-body problem[J].Chinese Physics Letters，2008，25（25）：2342-2345.

[11]McGehee R P.Some homoclinic orbits for the restricted three-body problem[D].Madison：University of Wisconsin，1969.

[12]Koon W S，Lo M W，Marsden J E，et al.Heteroclinic connections between periodic orbits and resonance transitions in celestial mechanics[J].Chaos，2000，10（2）：427-469.

[13]Koon W S，Lo M W，Marsden J E，et al.Low energy transfer to the moon[J].Celestial Mechanics&Dynamical Astronomy，2001，81（1-2）：63-73.

[14]Koon W S，Lo M W，Marsden J E，et al.Dynamical systems，the three-body problem，and space mission design [M].Berlin：World Scientific，2000：123-140.

[15]Gong S P，Li J F，Baoyin H X，et al.Lunar landing trajectory design based on invariant manifold[J].Applied Mathematics and Mechanics，2007，28（2）：201-207.

[16]Zhang P，Li J F，Baoyin H X，et al.A low-thrust transfer between the Earth-Moon and Sun-Earth systems based on invariant manifolds[J].Acta Astronautica，2013，91（10）：77-88.

[17]李俊峰，寶音賀西，蔣方華.深空探測動力學與控制[M].北京：清華大學出版社，2014：382-429.

[18]Anderson R L，Lo M W.Spatial approaches to moons from resonance relative to invariant manifolds[J].Acta Astronautica，2014，105（1）：355-372.

[19]Ren Y，Shan J.Low-energy lunar transfers using spatial transit orbits[J].Communications in Nonlinear Science&Numerical Simulation，2014，19（3）：554-569.

[20]Asano Y，Yamada K，Jikuya I.Approximating elliptic halo orbits based on the variation of constants[J].Acta Astronautica，2015，113：169-179.

[責任編輯 田 豐 夏紅梅]

Construction of explicit symplectic scheme in restricted three-body problem

WEI Bingwei，LI Yinshan
（School of Mechanical Engineering,Hebei University of Technology,Tianjin 300130,China）

The problem of construction of explicit symplectic scheme is investigated based on circular restricted threebody problem(CR3BP).To begin with,by separating original CR3BP Hamiltonian function into several systems with nilpotent of degree two,explicit symplectic Euler schemes with respect to different systems with nilpotent of degree two are found.Secondly,all the explicit symplectic Euler schemes are proved to be self-conjugate operators.As a result,explicit composite symplectic schemes are proposed in the paper by combining those Euler schemes.Finally,a numerical simulation study is conducted by using fourth-order explicit composite symplectic scheme and other integrators to calculate the regular orbit and halo orbit in CR3BP,and the availability and superiority of explicit symplectic scheme are verified.The results show that explicit symplectic scheme leads to oscillation of the system energy error in a certain range instead of error dissipation over the long term.Moreover,symplectic algorithms possess higher numerical accuracy than traditional Runge-Kutta methods with the same order.

symplectic geometry algorithm;restricted three-body problem;Hamiltonian system;composition scheme

P132.2；P138.2

A

1007-2373（2017）01-0040-08

10.14081/j.cnki.hgdxb.2017.01.007

2016-11-21

國家自然科學基金（10632040）

韋炳威（1989-），男，碩士研究生.

：李銀山（1961-），男，教授，博士.