孫 浪，劉福窯，王 穎，孫 威,3

(1.上海工程技術(shù)大學(xué) 數(shù)理與統(tǒng)計(jì)學(xué)院，上海201620；2.上海工程技術(shù)大學(xué) 計(jì)算物理與應(yīng)用研究中心，上海201620；3.安慶師范大學(xué) 資源環(huán)境學(xué)院，安慶246133)

1 引言

哈密頓系統(tǒng)從物理本質(zhì)上具有辛結(jié)構(gòu)的不變性。20世紀(jì)80年代初期，中國學(xué)者馮康先生[1]與Ruth[2]幾乎同時(shí)提出能夠保持哈密頓相流辛結(jié)構(gòu)的數(shù)值積分算法，解決了傳統(tǒng)算法因長期對時(shí)間積分不能保持系統(tǒng)的能量守恒的問題[4]。前者主要對不可分的哈密頓系統(tǒng)采用以隱式中點(diǎn)法為基礎(chǔ)來構(gòu)建高階隱式辛算法，后者主要是針對可分解為動(dòng)能T和勢能V的哈密頓系統(tǒng)建立顯式辛算法。顯式辛算法和隱式辛算法的區(qū)別在于積分過程中是否需要迭代計(jì)算。

在Ruth建立的T+V形式的顯式辛算法的基礎(chǔ)上，高階辛算法的研究和應(yīng)用得到了快速發(fā)展[5]。當(dāng)哈密頓分解成主要的未受攝部分H0和次要的攝動(dòng)部分εH1且兩部分均可積時(shí)[6]，可以通過升高攝動(dòng)項(xiàng)εH1的小參數(shù)ε的階提高精度，即類高階辛算法(pseudo-highorder-symplectic integrator,PSI)[7]和辛校正(wisdom-holman-touma correction,WHT)[8]。力梯度辛算法[6]是在算法的Lie算子中嵌入了力梯度算子，可以有效避免負(fù)時(shí)間系數(shù)的出現(xiàn)。Ruth[2]和Chin[9,10]分別構(gòu)造了三階力梯度辛算法和四階力梯度辛算法。李榮和伍歆[12–14]結(jié)合McLachlan[11]的算法優(yōu)化思想構(gòu)造了三階優(yōu)化力梯度辛算法和四階優(yōu)化力梯度辛算法。陳云龍和伍歆[15]論證了力梯度辛算法在求解旋轉(zhuǎn)坐標(biāo)系下的圓型限制性三體問題的有效性。

哈密頓系統(tǒng)變量不可分離時(shí)，顯式辛算法一般不能直接應(yīng)用，隱式辛算法應(yīng)該適合應(yīng)用。常用隱式辛算法有半隱式辛算法[16,17]和隱式中點(diǎn)法[16–23]。隱式辛算法適合任意哈密頓系統(tǒng)，但計(jì)算效率低，不宜作為首選數(shù)值方法。為使顯式辛算法在不可分哈密頓系統(tǒng)中得到應(yīng)用，Pihajoki[24]在不可分哈密頓系統(tǒng)中增加了一組與原相空間動(dòng)量坐標(biāo)相同的量，提出了相空間擴(kuò)充對稱算法，但系統(tǒng)的辛結(jié)構(gòu)被破壞。Tao[25]在Pihajoki的基礎(chǔ)上將哈密頓系統(tǒng)分成三部分，第三部分是兩組坐標(biāo)動(dòng)量差的平方組成的約束哈密頓，保證了數(shù)值解在擴(kuò)展相空間中的辛結(jié)構(gòu)。在此基礎(chǔ)上，Liu和Wu[26,27]構(gòu)造了連續(xù)坐標(biāo)動(dòng)量置換相空間擴(kuò)大法，算法適合求解后牛頓哈密頓問題和旋轉(zhuǎn)哈密頓問題等哈密頓系統(tǒng)不可分的問題，但在旋轉(zhuǎn)致密雙星后牛頓哈密頓系統(tǒng)中不能保持系統(tǒng)的能量穩(wěn)定。Luo等人[28,29]提出了中點(diǎn)置換相空間擴(kuò)大方法，這種方法優(yōu)越于Pihajoki方法和劉磊等人的方法，它對算法的算子個(gè)數(shù)沒有要求，適用范圍更廣。Wu等人[30]提出的擴(kuò)大相空間優(yōu)化Forest-Ruth算法在平面圓型限制性三體問題和致密雙星問題中與未優(yōu)化算法相比精度得到了明顯改善，且表明采用中點(diǎn)置換能使優(yōu)化的Forest-Ruth算法具有更好的數(shù)值表現(xiàn)。Li和Wu[31]在Mikkola等人[32,33]的基礎(chǔ)上結(jié)合擴(kuò)大相空間中的幾種置換方法，提出了擴(kuò)大相空間的對數(shù)哈密頓顯式對稱法。這些方法在圓型限制性三體問題、三階后牛頓自旋致密雙星系統(tǒng)和Ernst-Schwarzschild黑洞等模型中表現(xiàn)出了高精度和高效率。辛算法在解決各類天文學(xué)物理模型中得到了較為廣泛的應(yīng)用。本文的主要目的是對上述辛算法及類辛算法的發(fā)展、構(gòu)建以及適用物理模型進(jìn)行闡述。

2 辛算法的建立

按照T+V分解形式，假設(shè)哈密頓可以表示為：

式中，p,q分別表示廣義坐標(biāo)和廣義動(dòng)量。相應(yīng)的正則運(yùn)動(dòng)方程：

或取z=(p,q)

其中，L H是Lie算子，分別以A和B約定T和V的Lie導(dǎo)數(shù)：

式中，τ為時(shí)間步長。設(shè)W={,H}，將W的指數(shù)冪可以用A和B的指數(shù)冪多次組合而成，即

a i和b i是由特定的階確定的系數(shù)，O(h K+1)指哈密頓函數(shù)的K階截?cái)嗾`差。若使Lie算子的個(gè)數(shù)等于階條件的個(gè)數(shù)，時(shí)間系數(shù)可以很容易被算出。Ruth[2]在文獻(xiàn)中提出了比較常用的二階和三階辛算法。在此基礎(chǔ)上，F(xiàn)orest-Ruth四階非力梯度辛算法[3]和Yoshida高階非力梯度辛算法[5]也相繼被建立。

3 非力梯度顯式辛算法

3.1 T+V分解法

3.1.1 未優(yōu)化的非力梯度顯式辛算法

按照Ruth[2]的方法，對于可積分離的哈密頓系統(tǒng)，辛算法構(gòu)造如下，

(1)一階顯式辛算法

(2)二階顯式對稱辛算法

二階顯式對稱辛算法由三個(gè)單指數(shù)算子構(gòu)成，它的構(gòu)造形式為：

(3)三階顯式辛算法

三階顯式辛算法由6個(gè)單指數(shù)算子構(gòu)成，它的構(gòu)造形式為S3：

(4)四階顯式Forest-Ruth辛算法(fourth-order-Forest-Ruth symplectic integrator,FR)

Forest-Ruth算法從式(6)出發(fā)，它由7個(gè)單指數(shù)構(gòu)成S4：

構(gòu)造高階算法的另一個(gè)思路是直接對二階辛算法進(jìn)行多重積，Yoshida[5]在Ruth[2]的理論基礎(chǔ)上用三個(gè)二階辛算法對稱組合得到四階顯式Y(jié)oshida辛算法S YO：

同年Forest和Ruth[3]提出了Forest-Ruth四階辛算法，這兩種辛算法構(gòu)造可以達(dá)到等價(jià)的效果。按照上述算法構(gòu)造方式，高階辛算法S2n+2可以被構(gòu)造[5]，

3.1.2 優(yōu)化的非力梯度顯式辛算法

McLachlan[11]1992年提出了算法的優(yōu)化概念，通過將主要截?cái)嗾`差項(xiàng)的影響降到最小，得到優(yōu)化系數(shù)，從而構(gòu)造優(yōu)化的二階顯式辛算法和三階顯式辛算法。在此基礎(chǔ)上，四階優(yōu)化非力梯度算法[34]等高階優(yōu)化算法也被建立起來。

(1)二階優(yōu)化的顯式辛算法和三階優(yōu)化顯式辛算法

單一化系數(shù)表示的二階顯式辛算法的主要誤差函數(shù)為：

表1 幾種顯式辛算法的截?cái)嗾`差比較[11]

由此，McLachlan將n階辛算法的主要誤差函數(shù)表示為：

f j和g i表示誤差常數(shù)。令該式的系數(shù)平方和最小，得到三階最優(yōu)化辛算法O3：

(2)四階顯式優(yōu)化Forest-Ruth算法(fourth-order-optimal-Forest-Ruth integrator,OFR)Omelyan等人[34]將四階顯式Forest-Ruth算法拓展分解為：

ξ,χ,λ為3個(gè)步長系數(shù)，它們之間存在兩個(gè)約束條件，A的系數(shù)α(ξ,χ,λ)=0，B的系數(shù)β(ξ,χ,λ)=0。由于仍存在一個(gè)自由度，步長系數(shù)的取值將有多種解。從McLachlan的優(yōu)化思想出發(fā)，令四階顯式Forest-Ruth算法的三階截?cái)嗾`差為零及五階截?cái)嗾`差項(xiàng)系數(shù)平方和最小，得到最優(yōu)解：

該組合系數(shù)解下的OFR算法的截?cái)嗾`差值為γmin=0.000 92，四階FR辛算法的截?cái)嗾`差值為γFR=0.039，優(yōu)化后的精度得到了顯著提高。

(3)優(yōu)化的四階類Suziki算法(fourth-order-pseudo-Suziki integrator,SU)

Omelyan等人[34]在Yoshida構(gòu)建的四階辛算法兩端各加入兩個(gè)相同算子構(gòu)建了SU算法，它由5個(gè)算子組合而成：

考慮滿足階條件和截?cái)嗾`差h5系數(shù)平方和最小，得到的哈密頓截?cái)嗾`差值為γ=0.001 1，與OFR算法相比，誤差精度低一個(gè)量級。

3.2 H=H0+εH1分解法

其中，H0和H1可積，H0是主要項(xiàng)，H1是攝動(dòng)項(xiàng)，小參數(shù)ε表示H1與H0有量級ε的差別，具體分解過程在文獻(xiàn)[35]中有詳細(xì)描述。

規(guī)定A={,H0}和B={,H1}，由辛算法的構(gòu)造公式eτW=eτAeτB[2]并反復(fù)利用Baker-Campbell-Hausdroff(BCH)公式[36,37]得到哈密頓函數(shù)的誤差表達(dá)式：

a1,b1是誤差項(xiàng)系數(shù)，從上式可以得到三種提高辛積分器的精度方法：(1)提高τ的冪次得到高階辛算法；(2)改進(jìn)哈密頓的分解方式和構(gòu)造方式，ε越小，精度越高，比如哈密頓函數(shù)攝動(dòng)分解[6]、時(shí)間變換以及擴(kuò)展相空間；(3)升高ε的冪次，得到辛校正[8]。若誤差項(xiàng)中的O(ε2τn+1)仍存在，升高τ的冪次得到類高階辛方法[7]。這兩種類辛算法不需要減小積分步長，也不用升高積分的階就能優(yōu)化積分精度和效率。一般情況下，哈密頓函數(shù)按照H=H0+εH1分解時(shí)，與T+V分解相比，算法精度有顯著提高。如圖1所示，以純開普勒問題為例，攝動(dòng)分解時(shí)的算法的數(shù)值精度比T+V分解提高了接近5個(gè)量級。

圖1 純開普勒問題在兩種哈密頓分解形式中的能量誤差

3.2.1 辛校正

Wisdom和Holman[6]在雅可比坐標(biāo)系建立了如下的二階辛算法(second-order-wisdomverlet,MV2)：

其中，N=1,b0=b1=1/2,a0=1。MV2的數(shù)值結(jié)果對應(yīng)的哈密頓函數(shù)和真實(shí)哈密頓函數(shù)H的關(guān)系為=H+Herr。Wisdom等人[8]運(yùn)用Lie級數(shù)算子構(gòu)造一個(gè)正則變換的生成函數(shù)W，Herr不含攝動(dòng)小參數(shù)ε的一或者二次冪項(xiàng)，從而構(gòu)造辛校正公式。辛校正具有如下特點(diǎn)：(1)在不需要減少步長也不需要升高階數(shù)的情況下，通過消去誤差項(xiàng)中的ε或ε2項(xiàng)，實(shí)現(xiàn)一階或者二階辛校正來顯著提高數(shù)值精度；(2)辛校正不嚴(yán)格保辛，但沒有引入耗散機(jī)制；(3)計(jì)算效率高。

WHT計(jì)算生成函數(shù)W的過程較為繁瑣，而且校正子C只適用于上述辛算法，不具備一般性。Mikkola和Palmer[39]借助Euler-Mclaurin式推導(dǎo)出辛校正的生成函數(shù)W。同年，Chambers和Murisun[7]建立了類高階辛算法。Duncan等人[40,41]借助日心坐標(biāo)系將太陽系動(dòng)力學(xué)中的哈密頓函數(shù)分為可積的三部分。Malhotra[42]研究伽利略衛(wèi)星的潮汐演化時(shí)，將哈密頓系統(tǒng)分為了14個(gè)可積部分。Wu等人[43,44]令哈密頓函數(shù)的多部分分解更加一般化，將辛算法中的H0分成了N個(gè)子步，并用數(shù)值方法討論了哈密頓在具有N+1個(gè)可積部分時(shí)的算法精度(主要是一階和二階)。

Wu等人[44]利用伯努利多項(xiàng)式[7]推出辛算法的哈密頓的誤差函數(shù)：從哈密頓函數(shù)誤差項(xiàng)出發(fā)，消去ε或ε2項(xiàng)，直接確定校正子C而不計(jì)算生成函數(shù)W，分別得到一階辛校正和二階辛校正。

(1)ε的一階辛校正

由誤差函數(shù)得到一階辛校正算子[44]，從而消去截?cái)嗾`差項(xiàng)中εi的一階項(xiàng)，一階校正子的表達(dá)式為：

其中，A=L H0，B i=L Hi，互易子[A,B]=AB?BA，[A,B,C]=[A,[B,C]]，又記[A,B]=?AB，[A,A,B,C,C]=?2A?B?CC，B(k)(x)是伯努利公式，滿足遞推公式[8]：

式(20)適用于有n個(gè)可積部分的哈密頓系統(tǒng)的任意階辛算法，哈密頓的截?cái)嗾`差可以通過matlab或mathematica等數(shù)學(xué)軟件直接確定[44]，提高了計(jì)算效率。

(2)ε的二階辛校正

Wu等人[44]討論了二階顯式辛算法S2X的二階辛校正過程，該推導(dǎo)過程適用于一般n階辛算法。S2X關(guān)于ε2的校正公式如下：

其中，

C2可用Cˉ2近似代替[44]。在太陽-木星-土星模型中，從能量誤差、計(jì)算效率和平均經(jīng)度誤差方面進(jìn)行數(shù)值比較，辛校正的計(jì)算精度比顯式辛算法提高了一個(gè)量級并具有更好的數(shù)值穩(wěn)定性。

3.2.2 類高階辛算法

類高階辛算法是在攝動(dòng)分解的基礎(chǔ)上形成高效的積分方法。相比較高階辛算法的子步隨時(shí)間步長的增加而迅速增加，類高階辛算法的子步更少且計(jì)算效率更高。Chambers和Murison[7]從四階和六階辛算法出發(fā)，通過構(gòu)造合適系數(shù)使誤差項(xiàng)中不含ετ4或ετ6，構(gòu)造過程如下：

當(dāng)

時(shí)，可消去ετ3項(xiàng)。

文中列舉了類四階辛算法和類六階辛算法[7]：

八階和十階等高階算法也可被構(gòu)造。從算法的構(gòu)成特點(diǎn)看，類高階辛算法能保持系統(tǒng)能量穩(wěn)定[45,46]，計(jì)算效率也明顯優(yōu)于同階的通常辛算法[11,39]。從算法的截?cái)嗾`差項(xiàng)來看，它可以看作是對二階辛算法的部分校正。類高階辛算法到達(dá)一定階時(shí)，精度將不再隨階數(shù)的增加而提高。Laskar和Robutel[47]指出類六階或八階的數(shù)值效果最好。Wu等人[44]在日-木-土三體模型中對類四階辛算法PSI4、它的校正算法及通常四階辛算法S4三種算法進(jìn)行數(shù)值模擬，結(jié)果表明三者的能量誤差精度表現(xiàn)為同一量級。

4 力梯度辛算法

通常辛算法在積分過程中難免會(huì)出現(xiàn)負(fù)積分步長，但在不可逆的力學(xué)問題中需使得積分步長全為正。Ruth[2]在三階辛算法的算子組合中引入力梯度算子使得每個(gè)積分子步的系數(shù)均為正。Chin等人[9,10,48–50]把力梯度辛算法發(fā)展到四階，與同階的非力梯度辛算法相比精度提高了10～80倍[9]，在圓型限制性三體問題、含時(shí)引力場問題[10]及量子力學(xué)問題[48]中得到了很好的應(yīng)用。Sun等人[51]在Chin提出的四階力梯度辛算法基礎(chǔ)上構(gòu)造了兩種含有三階導(dǎo)數(shù)項(xiàng)的四階力梯度辛算法，在Hnon-Heiles系統(tǒng)、四極矩核-殼模型以及含扁率限制性三體問題中進(jìn)行數(shù)值模擬，精度和計(jì)算效率都得到了提高。Xu和Wu[52,53]應(yīng)用這種力梯度辛算法求解攝動(dòng)二體問題和雅可比坐標(biāo)系下的攝動(dòng)N體問題。Omelyan等人[54,55]將最優(yōu)化思想引入到力梯度辛算法，按照對稱組合的形式構(gòu)造了三階和四階最優(yōu)化力梯度辛算法，并將這兩種方法運(yùn)用到天體力學(xué)、分子動(dòng)力學(xué)和量子力學(xué)等領(lǐng)域，與通常辛算法相比，數(shù)值精度得到了明顯的改善，但是步長系數(shù)還是不可避免地用到了負(fù)數(shù)。李榮等人[12–14,56,57]構(gòu)造了步長系數(shù)全為正步長的不對稱的三階、四階最優(yōu)化力梯度辛算法，這種方法在求解攝動(dòng)開普勒混沌問題[57]的能量精度和定態(tài)薛定諤方程能量本征值問題中[12–14,56]表現(xiàn)出明顯的精度優(yōu)勢。陳云龍和伍歆[15]論證了質(zhì)心旋轉(zhuǎn)坐標(biāo)系中動(dòng)能T不是動(dòng)量P的嚴(yán)格二次型時(shí)，力梯度辛算法仍然適合求解。

定義力梯度算子為：

其中A、B如式(4)所設(shè)，將算子C與A、B算子組合可構(gòu)造出力梯度辛算法：

4.1 非優(yōu)化力梯度辛算法

(1)三階力梯度辛算法F3

這種算法的積分步長都為正，但力梯度算子C的計(jì)算過程較復(fù)雜，早期很少受到重視。

(2)四階力梯度辛算法F4

Chin及其合作者在Ruth的基礎(chǔ)上將力梯度辛算法發(fā)展到了四階[9,10,48–50]，

與通常四階辛算法相比，精度提高了10～80倍。

Xu和Wu[52]按照Yoshida構(gòu)造標(biāo)準(zhǔn)高階辛算法的思路將算子C與A,B按照兩種不同方式對稱組合得到了兩類新的四階力梯度辛算法。

(a)力梯度算子嵌套在正中間時(shí)，表示為：

(b)力梯度算子嵌套在兩端時(shí)，表示為：

根據(jù)BCH公式，從上式組合算子的中間開始實(shí)施循環(huán)推導(dǎo)得到W的具體表達(dá)式，由最優(yōu)化思想得到一系列新的正積分系數(shù)，如表2和表3所示。A1和B1型算法對應(yīng)于Chin的四階C和D算法。Xu和Wu[52]用這兩種新的四階力梯度辛算法求解攝動(dòng)開普勒問題和雅可比坐標(biāo)系下的N體問題，以及深入研究了行星磁氣圈內(nèi)的帶電粒子的混沌運(yùn)動(dòng)。如圖2所示，兩種分解形式中力梯度辛算法的數(shù)值性能沒有明顯差異，攝動(dòng)分解形式中算法精度有顯著提高。攝動(dòng)分解形式中的力梯度算子的計(jì)算比計(jì)算H0更加簡單，與H1相比不會(huì)增加大量額外的計(jì)算成本，所以是非常值得推薦其用于研究雅可比坐標(biāo)系下的N體哈密頓問題[15,53]。

圖2 FR和A,B兩種類型四階力梯度辛算法在兩種哈密頓分解方式中的能量誤差[53]

表2 A型四階力梯度辛算法的幾組正積分系數(shù)[53]

表3 B型四階力梯度辛算法的幾組正積分系數(shù)[53]

4.2 優(yōu)化力梯度辛算法

(1)三階最優(yōu)化力梯度辛算法

Lie算子系數(shù)的求解只適用于對稱辛算法，給不對稱的力梯度辛算法的構(gòu)造增加了困難。李榮等人構(gòu)造了不對稱的三階和四階最優(yōu)力梯度辛算法[12–14,52]。根據(jù)Chin的思路，文中構(gòu)造了如下的不對稱算子組合：

根據(jù)BCH公式得到階條件的關(guān)系式，由最優(yōu)化思想，即加入附加的約束條件[13](令其四階截?cái)嗾`差系數(shù)平方和最小)，得到新的三階最優(yōu)化力梯度辛算法：

和

李榮[12]分別在諧振子模型、單擺模型、開普勒二體模型、Hnon-Heiles系統(tǒng)以及量子力學(xué)模型中對新構(gòu)造的三階最優(yōu)化力梯度辛算法與其他辛算法進(jìn)行數(shù)值比較。在經(jīng)典力學(xué)模型中，該算法的能量誤差計(jì)算占有絕對的精度優(yōu)勢，并且可以準(zhǔn)確地識別Hnon-Heiles系統(tǒng)的混沌軌道，因此可以有效避免數(shù)值誤差引起的虛假混沌現(xiàn)象。在定態(tài)薛定諤方程的能量本征值問題中，該算法同樣有明顯的精度優(yōu)勢，甚至與四階標(biāo)準(zhǔn)辛算法相比具有更高精度。

(2)四階最優(yōu)化力梯度辛算法

Omelyan等人[54,55]構(gòu)造的對稱的四階最優(yōu)化力梯度辛算法如下：

李榮[12]將得到的三階最優(yōu)化力梯度辛算法以對稱組合的方式得到兩個(gè)四階辛算法，

和

在開普勒二體問題和攝動(dòng)二體問題等經(jīng)典力學(xué)模型，以及量子力學(xué)模型中的Morse勢模型中進(jìn)行數(shù)值分析比較，新的四階力梯度算法表現(xiàn)出了比較好的數(shù)值精度，甚至還要明顯優(yōu)越于已有的四階最優(yōu)化力梯度辛算法OF4[55]。

旋轉(zhuǎn)坐標(biāo)系中的平面圓型限制性三體問題的動(dòng)能不是動(dòng)量的嚴(yán)格二次型，而是受到旋轉(zhuǎn)坐標(biāo)系影響存在動(dòng)量與坐標(biāo)的交叉項(xiàng)，這給力梯度算子的加入帶來了困難。陳云龍和伍歆[15]嚴(yán)格論證了力梯度算子在旋轉(zhuǎn)坐標(biāo)系中的形式仍與質(zhì)心坐標(biāo)中相同，均是引力的梯度而不是引力與非慣性力所得合力的梯度，只是將A算子變成了如下算子，式(26)中的C算子仍然適用，

應(yīng)用四階力梯度辛算法、最優(yōu)化四階力梯度辛算法和Forest-Ruth辛算法分別求解該問題，數(shù)值結(jié)果顯示最優(yōu)化四階力梯度辛算法能夠取得最好精度。

5 隱式辛算法

5.1 隱式辛算法的構(gòu)造

當(dāng)系統(tǒng)的哈密頓變量不可分離時(shí)，顯式算法一般無法直接應(yīng)用?？紤]不可分離哈密頓函數(shù)：

與3.2節(jié)中的約定類似，約定Lie算子A={,H0}和B={,H1}。

正則方程表示為：

以時(shí)間τ為步長，對整個(gè)哈密頓系統(tǒng)按照一階Euler部分向前差分公式，得到一階隱式辛算法M1和M?1、二階隱式中點(diǎn)法M2[2]：

(1)一階位置隱式Euler法

(2)一階動(dòng)量隱式Euler法

(3)二階隱式中點(diǎn)法

(4)二階顯隱混合辛算法

A算子是有解析解的算子，而B是隱式中點(diǎn)法迭代求解的算子，兩者結(jié)合組成：

(5)高階隱式辛算法

按照Yoshida[5]和Ruth[2]構(gòu)造高階顯式辛算法的思路對全局哈密頓函數(shù)構(gòu)造高階隱式辛算法：

和

全局隱式辛算法使計(jì)算效率不理想。若通過變量分離使得H0可積，H1不可積，再根據(jù)式(44)―(46)構(gòu)造顯隱式混合辛算法，顯然顯隱式混合辛算法的效率要高于全隱辛算法。

Zhong等人[18–21]將二階隱式中點(diǎn)法及其對稱組合運(yùn)用求解整個(gè)后牛頓哈密頓系統(tǒng)，數(shù)值試驗(yàn)表明，顯隱混合中點(diǎn)嵌入法在能量精度方面始終優(yōu)于全隱中點(diǎn)法，且前者的計(jì)算效率更好。Zhong和Wu[19]另外對中點(diǎn)嵌入法進(jìn)行優(yōu)化，進(jìn)一步提高了計(jì)算精度。對于自旋后牛頓哈密頓系統(tǒng)：

開普勒部分和后牛頓部分分別表示為：

對旋轉(zhuǎn)標(biāo)量引入Wu和Xie[18]的正則旋轉(zhuǎn)坐標(biāo)，得到新的正則哈密頓函數(shù)Γ(r,θ,p,ξ)，將隱式中點(diǎn)法直接應(yīng)用于求解變量Γ(r,θ,p,ξ)，記作S2A。另外，Zhong和Wu[19]分別對開普勒部分和后牛頓部分求解析解和數(shù)值解，數(shù)值解部分應(yīng)用隱式中點(diǎn)法，按照Yoshida的思路構(gòu)建四階正則顯隱式混合辛算法：

Mei等人[22,23]對四階Forest-Ruth型顯隱混合辛算法和四階Yoshida顯隱混合辛算法進(jìn)行理論分析，證明前一種算法在大多數(shù)情況下只有二階精度，后者可以達(dá)到四階精度。在一維不可分系統(tǒng)和致密雙星后牛頓哈密頓系統(tǒng)中得到了相同的結(jié)論。他們指出，在實(shí)際計(jì)算中，當(dāng)哈密頓的部分變量可積且解析解容易獲得時(shí)，應(yīng)采用S4(FR)算法或?qū)⒔馕鼋馀c隱式數(shù)值解結(jié)合的Yoshida算法S4(YO)，以獲得較高的計(jì)算精度和效率。當(dāng)可分哈密頓的變量部分可積，但其解析解難以求解，或者哈密頓變量部分不可積，則采用將解析解與隱式數(shù)值解結(jié)合的Yoshida算法S4(YO)。

Lubich等人[58]提出了四階非正則的顯隱式混合辛算法，主要在自旋致密雙星后牛頓哈密頓系統(tǒng)中應(yīng)用，將哈密頓分成了三大部分，分別是軌道HOrb、自旋?軌道HSO和自旋?自旋HSS：

按照Suzuki[17]用二階積分器的五重積構(gòu)建四階算法的思路構(gòu)造如下算法：

再分別將這三部分分解為2,3,4部分運(yùn)用顯隱式混合辛算法求解。將式(54)發(fā)展到四階：

5.2 幾種隱式辛算法的數(shù)值穩(wěn)定性比較

當(dāng)哈密頓兩部分均可積時(shí)，劉福窯等人[59]分別對上述一階隱式辛算法M?1、隱式中點(diǎn)法M2以及一階顯式算法，二階顯式算法四種算法的數(shù)值穩(wěn)定性進(jìn)行分析比較，將這幾種算法作用于線性哈密頓系統(tǒng)：

與T+V分解相比，H0+εH1分解中各算法的穩(wěn)定區(qū)域擴(kuò)大了。2009年，劉福窯和錢曉明[60]比較了Liao[16]提出的兩種顯隱混合辛算法的數(shù)值穩(wěn)定性，即將一階隱式辛算法M?1和隱式中點(diǎn)法M2分別嵌入到一階和二階顯辛算法中得到顯隱混合辛算法，前者的數(shù)值穩(wěn)定性要優(yōu)于后者。

鐘雙英[61]分別以一維耦合振子、經(jīng)典圓型限制性三體問題和后牛頓近似的致密雙星系統(tǒng)為數(shù)值研究對象，對顯隱混合辛算法MSI1和MSI2的數(shù)值優(yōu)劣性能進(jìn)行比較。一維耦合振子模型和圓型限制性三體問題中MSI2算法的穩(wěn)定性均要明顯優(yōu)于MSI1算法。兩種算法在不同的哈密頓分解形式中表現(xiàn)有很大不同，在T+V分解中，兩種嵌入法都能保持?jǐn)?shù)值穩(wěn)定；在H0+H1分解中，MSI1算法在有序軌道和混沌軌道中均不能保持?jǐn)?shù)值穩(wěn)定。自旋致密雙星后牛頓哈密頓系統(tǒng)中，MSI1算法的計(jì)算效率比MSI2算法稍高，但穩(wěn)定性不如MSI2算法。綜合諸因素可知，中點(diǎn)嵌入法更適合于求解各種相對論后牛頓動(dòng)力學(xué)問題。

隱辛算法適用于任何哈密頓系統(tǒng)，尤其在使用經(jīng)典算法表現(xiàn)出較差的穩(wěn)定性時(shí)，隱辛算法是理想的選擇，但是求解過程中由于迭代使得計(jì)算效率大大降低，因此應(yīng)盡量不要作為首選數(shù)值方法。

6 擴(kuò)展相空間顯辛算法

Pihajoki[24]提出一種相空間擴(kuò)充方法。在不可分離的哈密頓系統(tǒng)中分別增加一組與原空間相同的動(dòng)量和坐標(biāo)得到一個(gè)新的可分離的哈密頓量：

上式右邊的兩部分相互獨(dú)立且可積，可以分別進(jìn)行解析求解，由正則方程得到：

或

其中，M1,M2表示兩組變量之間的置換映射，即：

此外，需借助投影算符W將擴(kuò)展相空間的解返回到原始相空間，

6.1 高階連續(xù)坐標(biāo)動(dòng)量置換相空間擴(kuò)充顯式類辛算法

Liu和Wu[26]用偶數(shù)重積構(gòu)建了不同于Yoshida構(gòu)建高階辛算法思路的顯式高階類辛算法，用六個(gè)偶數(shù)低階算法構(gòu)造高階算法S4A：

這種構(gòu)建模式由四部分組成：置換坐標(biāo)，沒有置換的三個(gè)偶數(shù)低階蛙跳格式的積分，置換動(dòng)量，沒有置換的三個(gè)偶數(shù)低階蛙跳格式的積分[26]。

在平面圓型限制性三體問題、Chin[9]考慮的簡單模型和無旋轉(zhuǎn)致密雙星的后牛頓哈密頓系統(tǒng)中，高階連續(xù)坐標(biāo)動(dòng)量置換相空間擴(kuò)充顯式類辛算法在數(shù)值精度上優(yōu)于Yoshida的四階顯式辛算法、Chin[9]的顯式辛算法以及顯隱式混合辛算法，計(jì)算效率也遠(yuǎn)遠(yuǎn)高于顯隱式混合辛算法。這類顯式類辛算法也適合求解后牛頓哈密頓問題和旋轉(zhuǎn)哈密頓問題等坐標(biāo)和動(dòng)量不可分離的哈密頓系統(tǒng)。

6.2 中點(diǎn)置換相空間擴(kuò)充顯式類辛算法

Liu和Wu[26]提出的高階類辛算法需要進(jìn)行兩次三重積才能達(dá)到高精度效果，且在旋轉(zhuǎn)致密雙星后牛頓哈密頓系統(tǒng)的某些軌道積分中會(huì)失敗[52]。駱俊杰等人[28,29]將連續(xù)坐標(biāo)動(dòng)量置換改為一次中點(diǎn)置換(原變量與它對應(yīng)擴(kuò)展變量之間的中點(diǎn))，每一步積分都調(diào)整原變量及其對應(yīng)擴(kuò)展變量的值為他們的中點(diǎn)值，即：

再結(jié)合Yoshida的三重積構(gòu)造高階類辛算法S4B：

以不考慮引力耗散的二階后牛頓近似的自旋致密雙星為模型，在周期軌道和混沌軌道中，新構(gòu)造的中點(diǎn)置換相空間擴(kuò)充法S4B的數(shù)值精度與隱式中點(diǎn)法IM4、連續(xù)坐標(biāo)動(dòng)量置換相空間擴(kuò)充法S4A相比最高。在混沌軌道中，連續(xù)坐標(biāo)動(dòng)量置換相空間擴(kuò)充法S4A在積分一半時(shí)間時(shí)產(chǎn)生了劇烈的誤差偏移。由于H1,H2的不同，在混沌軌道中，使得S4A在置換過程中兩部分哈密頓函數(shù)的誤差不停積累導(dǎo)致失效。

中點(diǎn)置換法可以推廣到非保守的引力耗散哈密頓系統(tǒng)中，分別以引力耗散的2.5PN后牛頓近似的自旋致密雙星、阻尼諧振子和受到Poynying-Robertson阻力的塵埃粒子為模型，中點(diǎn)置換法仍然保持了最佳的優(yōu)越性，該方法的計(jì)算精度、計(jì)算效率以及普適性都優(yōu)于其他算法，值得被推廣使用。

6.3 擴(kuò)大相空間FR最優(yōu)化算法

吳亞林等人結(jié)合優(yōu)化思想和擴(kuò)充相空間的思想，構(gòu)造了幾種優(yōu)化的相空間擴(kuò)充類辛算法；并分別在一維不可分可積系統(tǒng)、平面圓型限制性三體問題、三階后牛頓致密雙星系統(tǒng)中與其它辛算法進(jìn)行了數(shù)值比較[30,62]。

6.3.1 幾種擴(kuò)充相空間類辛算法

(1)連續(xù)坐標(biāo)動(dòng)量置換擴(kuò)充相空間優(yōu)化算法

將FR算法和其優(yōu)化算法OFR中間算子按照Liu和Wu[26]構(gòu)造的高階類辛算法進(jìn)行拆分，分別將(10)式和(16)式改造為：

(2)中點(diǎn)置換擴(kuò)充相空間優(yōu)化算法

按照駱俊杰[29]提出的中點(diǎn)置換法將式(10),(11),(16),(17)分別改造為：

(3)Tao保辛顯式算法

Tao[25]改造了Pihajoki的方法，將不可分哈密頓系統(tǒng)改造為：

上式ωH3部分是Tao人為增加的約束哈密頓函數(shù)，ω是控制新舊相空間變量之間間隔的參數(shù)，H1和H2有解析解，上文第二節(jié)中給出的算子A,B算子仍然適用，H3部分可積。這里再令算子Bˉ為：

其中D算子是H3部分的算子。得到新的二階顯式辛算法：

吳亞林等人對FR算法以及Suziki算法用Tao的方法進(jìn)行改造得：

由于算法中沒有加入置換因子，算法的辛結(jié)構(gòu)將不會(huì)被破壞，ωH3起到約束和調(diào)節(jié)新舊空間變量的作用，ω的選取對精度將產(chǎn)生很大的影響。

6.3.2 算法的數(shù)值檢驗(yàn)

分別在一維不可分哈密頓系統(tǒng)、平面圓型限制性三體問題和三階后牛頓非自旋致密雙星系統(tǒng)中對以上算法數(shù)值比較，精度由高到低是M OFR≈M SU>M FR≈M YO>IM4>EOFR≈TSU>EFR≈TYO?？偟膩碚f，優(yōu)化后的算法比未優(yōu)化的精度高，各算法與中點(diǎn)置換結(jié)合比與Tao法結(jié)合精度高，能達(dá)到優(yōu)化四階通常顯辛算法的精度，且四階隱式辛算法的計(jì)算效率最低。TSU算法與系數(shù)ω有關(guān)，ω=100時(shí)精度達(dá)到四階精度。

6.4 擴(kuò)大相空間的對數(shù)哈密頓顯式對稱法

對數(shù)哈密頓辛算法在解決高偏心率軌道問題時(shí)可以通過調(diào)節(jié)步長避免數(shù)值失真，構(gòu)造過程較一般的時(shí)間變換辛算法更簡單。擴(kuò)大相空間使得顯式辛算法可以應(yīng)用于不可分哈密頓系統(tǒng)，且精度理想。Li和Wu[31,63]結(jié)合這兩種思想，在Mikkola等人[32,64,65]提出的對數(shù)哈密頓方法的基礎(chǔ)上加上一個(gè)常數(shù)或者函數(shù)，再結(jié)合Pihajoki的相空間擴(kuò)充思想，構(gòu)造了適用范圍更廣的擴(kuò)充相空間的對數(shù)哈密頓類辛算法，它的基本的哈密頓形式為：

其中，p0=?H，q0=t，C為積分常數(shù)，參數(shù)C>0。當(dāng)C=0時(shí)，是Mikkola以及蘇湘寧[66]所構(gòu)建和使用的對數(shù)哈密頓方法。時(shí)間變換函數(shù)g=T+p0+C和g=U+C。哈密頓的動(dòng)能和勢能部分可分時(shí)，一般的顯式辛算法可以使用；兩者不可分時(shí)，擴(kuò)充相空間的思想可以被考慮。

6.4.1 擴(kuò)大相空間的對數(shù)哈密頓辛算法

根據(jù)系統(tǒng)顯含時(shí)間t和不顯含時(shí)間t可將哈密頓函數(shù)分別擴(kuò)充為：

(1)不顯含時(shí)間

6.4.2 數(shù)值檢驗(yàn)

將四階顯式辛算法FR、四階隱式中點(diǎn)法IS4(與6.3.2節(jié)中的IM4相同)及三種擴(kuò)充相空間算法S4A,S4B,S4C分別應(yīng)用求解開普勒攝動(dòng)二體問題、圓型限制性三體問題和無自旋三階后牛頓致密雙星系統(tǒng)。在三種模型中，S4C算法的表現(xiàn)始終是最好的，隱式中點(diǎn)法IS4算法雖然也表現(xiàn)出了很好的計(jì)算精度，但是在計(jì)算成本上最高；S4B算法與IS4算法精度同量級且計(jì)算成本最低。擴(kuò)大相空間的對數(shù)哈密頓顯式類辛算法在數(shù)值穩(wěn)定性、數(shù)值精度以及計(jì)算效率上均占有很大的優(yōu)勢，且適用于高偏心率軌道。本文主要在圓型限制性三體問題中對以上算法進(jìn)行數(shù)值比較，這幾種算法分別是四階顯式辛算法FR、四階隱式中點(diǎn)法IS4、連續(xù)坐標(biāo)動(dòng)量置換相空間擴(kuò)充顯式算法EFR及優(yōu)化算法EOFR、中點(diǎn)置換相空間擴(kuò)充顯式算法MFR及優(yōu)化算法MOFR。如圖3所示，對于低偏心率軌道，四階顯式辛算法FR在長時(shí)間積分過程中能保持能量誤差穩(wěn)定且精度較高，但在高偏心率軌道中，F(xiàn)R算法在長時(shí)間積分后能量誤差開始增長。c=0時(shí)，連續(xù)坐標(biāo)動(dòng)量置換EFR四階算法與其優(yōu)化算法EOFR均不適用于高偏心率軌道，它們在較短時(shí)間誤差就開始增長。c=1.7時(shí)，各相空間擴(kuò)大法的精度均有提高，且EFR算法和EOFR算法在高偏心率軌道中使得能量誤差不再線性增長，說明系數(shù)c的選取很重要。另外，中點(diǎn)置換相空間擴(kuò)大法MFR及其優(yōu)化算法MOFR的計(jì)算效率最高，四階隱式算法IS4的計(jì)算效率最低。

圖3 圓型限制性三體問題中的兩條軌道的能量誤差

7 總結(jié)與展望

辛算法自從被提出以來，特別是可分哈密頓系統(tǒng)的顯式辛算法，在動(dòng)力學(xué)天文中被廣泛應(yīng)用，已充分顯示出傳統(tǒng)算法不可比擬的優(yōu)點(diǎn)。在定量天文學(xué)計(jì)算問題中，辛算法可以有效控制軌道沿跡誤差的持續(xù)快速增長，但是對于同階算法，辛算法的截?cái)嗾`差要比Runge-Kutta法等傳統(tǒng)方法大，在較短時(shí)間內(nèi)優(yōu)點(diǎn)無法突出。辛校正和類高階辛算法則通過誤差校正，在不縮小積分步長的前提下使得算法精度提高，時(shí)間變化辛算法使用變步長，解決了天體緊密交會(huì)和大偏心率軌道中的數(shù)值解失真問題。對于不可分哈密頓系統(tǒng)，隱式辛算法非常適合應(yīng)用，與傳統(tǒng)算法相比，它也能保證系統(tǒng)的辛結(jié)構(gòu)和能量誤差穩(wěn)定，但積分過程由于多次迭代使得計(jì)算效率顯著降低。擴(kuò)大相空間的類辛算法近年來得到了快速的發(fā)展，其中Luo等人[28,29]構(gòu)造的中點(diǎn)置換擴(kuò)充相空間方法是目前最理想的，它具有高精度、高效率和高適用范圍，之后提出的擴(kuò)大相空間FR最優(yōu)化算法[30,62]以及擴(kuò)大相空間對數(shù)哈密頓顯式類辛算法[31,63]都對此進(jìn)行了證明。力梯度辛算法作為一種獨(dú)特的辛算法，在時(shí)間可逆的哈密頓系統(tǒng)求解中具有很大的優(yōu)勢，但在哈密頓變量不可分時(shí)應(yīng)用有困難。