限制性三體問題中顯式辛格式的構(gòu)造

韋炳威，李銀山

（河北工業(yè)大學(xué) 機(jī)械工程學(xué)院，天津 300130）

限制性三體問題中顯式辛格式的構(gòu)造

韋炳威，李銀山

（河北工業(yè)大學(xué) 機(jī)械工程學(xué)院，天津 300130）

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

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

0 引言

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

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

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

1 辛可分哈密頓系統(tǒng)

1.1 基本理論[6-7]

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

式中φi（u）是n個(gè)變量的標(biāo)量函數(shù)，且上式滿足條件

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

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

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

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

其中有效勢能U（q）為

在計(jì)算共線拉格朗日點(diǎn)附近的軌道時(shí)往往將坐標(biāo)系的原點(diǎn)平移到相應(yīng)的拉格朗日點(diǎn).本文以地-月系統(tǒng)的拉格朗日2（L2）點(diǎn)為例進(jìn)行計(jì)算研究，故將上述坐標(biāo)系的原點(diǎn)平移到L2點(diǎn)，再取L2點(diǎn)到月球的距離γ為1個(gè)單位，哈密頓函數(shù)改寫為

現(xiàn)將哈密頓函數(shù)（10）作如下的分解：

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

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

2 顯式辛格式的構(gòu)造

2.1 共軛積分算子

由上述結(jié)論，CR3BP的哈密頓函數(shù)可以分解成6個(gè)二階冪零哈密頓系統(tǒng)

根據(jù)（5）式，可以通過復(fù)合6個(gè)二階冪零系統(tǒng)的歐拉格式得到1階顯式格式

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

其中因?yàn)棣読是2次的n元函數(shù)，Dφi是一個(gè)線性算子，可以用一個(gè)矩陣Mi表示.考慮到條件（2），有

2.2 高階組合格式

3 實(shí)例軌道計(jì)算

3.1 一般軌道計(jì)算

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

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

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

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

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

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

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

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

圖5 地-月系統(tǒng)一般軌道Fig.5 Regular orbit of Earth-Moon system

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

3.2 Halo軌道計(jì)算

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

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

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

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

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

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

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

4 結(jié)論

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

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

本文提出的顯式辛格式的構(gòu)造方法不僅僅適用于CR3BP，也同樣適用于雙圓型限制性四體問題，原因是本質(zhì)上該問題也是辛可分哈密頓系統(tǒng).因?yàn)樵撓到y(tǒng)不再是自治系統(tǒng)，故不存在能量守恒，但可以通過擴(kuò)充原有的哈密頓系統(tǒng)維數(shù)使其變成自治系統(tǒng)，從而構(gòu)造形式上的守恒“能量”，余下的分析方法與本文討論的一致.

[1]馮康，秦孟兆.哈密頓系統(tǒng)的辛幾何算法[M].杭州：浙江科學(xué)技術(shù)出版社，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]陳云龍，伍歆.力梯度辛方法在圓型限制性三體問題中的應(yīng)用[J].物理學(xué)報(bào)，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]李俊峰，寶音賀西，蔣方華.深空探測動力學(xué)與控制[M].北京：清華大學(xué)出版社，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.

[責(zé)任編輯 田 豐 夏紅梅]

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

國家自然科學(xué)基金（10632040）

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

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