趙楊華，馬大柱

(湖北民族學(xué)院 理學(xué)院，湖北 恩施 445000)

復(fù)雜動(dòng)力系統(tǒng)的非線(xiàn)性行為研究依賴(lài)于數(shù)值計(jì)算方法的穩(wěn)定性和精度.傳統(tǒng)低階數(shù)值算法在一定程度上能反映動(dòng)力系統(tǒng)的相空間結(jié)構(gòu)，但是對(duì)于某些系統(tǒng)的混沌行為，由于其對(duì)初始條件極為敏感，并且對(duì)數(shù)值結(jié)果的精度要求非常高，所以尋找合適的數(shù)值計(jì)算方法是當(dāng)前非線(xiàn)性研究工作的基本問(wèn)題.為避免數(shù)值誤差產(chǎn)生的假混沌現(xiàn)象，目前文獻(xiàn)中主要采取以下三種方式解決：一是采用高階算法求解微分方程(組)；二是采用辛算法[1]；三是對(duì)傳統(tǒng)低階數(shù)值算法加以改造，即后穩(wěn)定化方法[2].研究表明，高階算法由于截?cái)嗾`差小，計(jì)算耗時(shí)較長(zhǎng)，不利于廣泛應(yīng)用.辛算法雖然能保持系統(tǒng)的辛結(jié)構(gòu)，但是其主要用于處理可分的哈密頓系統(tǒng)，所以應(yīng)用有限.相對(duì)而言，后穩(wěn)定化方法不僅節(jié)省了計(jì)算時(shí)間，而且能提高數(shù)值積分的穩(wěn)定性和精度.

1 能量穩(wěn)定化方法

對(duì)一含有不變運(yùn)動(dòng)積分的微分動(dòng)力系統(tǒng)：φ(x)=c，c為積分常數(shù)，理論情況下應(yīng)滿(mǎn)足：

ε(x)=φ(x)-c=0

(1)

但是數(shù)值誤差的存在，實(shí)際計(jì)算導(dǎo)致ε≠0.原因在于數(shù)值算法存在截?cái)嗾`差以及計(jì)算機(jī)將舍入誤差代入到各種算法中.為解決該問(wèn)題，Nacozy[3]在1971年提出了流形改正方案，即后穩(wěn)定化方法，以不變積分常數(shù)作為控制量，將數(shù)值算法得到的超曲面用最小二乘法的原理以最短路徑拉回到原始的超曲面上.

簡(jiǎn)單而言，假設(shè)x代表數(shù)值方法得到的坐標(biāo)和速度矢量，x′為改正后的坐標(biāo)和速度矢量，兩者之間應(yīng)該滿(mǎn)足x′=x+Δη，Δη即為改正量，Δη必須滿(mǎn)足最小二乘法原理，并且是以最短的距離，所以Δη=-ET(EET)-1ε，式中上標(biāo)T代表轉(zhuǎn)置矩陣，其中E滿(mǎn)足：

(2)

多次研究表明，該方案非常適合用于穩(wěn)定常微分方程模型.將后穩(wěn)定化方法推廣到含有能量積分模型中即為能量穩(wěn)定化方法[4]，尤其是只存在單個(gè)能量積分情況下，能量穩(wěn)定化可以簡(jiǎn)寫(xiě)為標(biāo)度因子的形式，即：

x′=sx

(3)

s的形式由具體運(yùn)動(dòng)的哈密頓形式給出.本文基于以上結(jié)論，進(jìn)一步推廣該方案，將雅克比積分常數(shù)看做不變積分應(yīng)用能量穩(wěn)定化方法討論平面圓形限制性三體問(wèn)題.

2 平面圓形限制性三體問(wèn)題模型

圖1 雅克比積分隨時(shí)間的變化關(guān)系

平面圓形限制性三體問(wèn)題模型[5]是三體其中一體的質(zhì)量遠(yuǎn)小于其它兩體的特殊三體問(wèn)題模型.可以不考慮質(zhì)量最小體對(duì)其它天體的引力作用，于是兩主天體在相互引力的作用下作圓形運(yùn)動(dòng)，對(duì)于兩主天體引力下的第三體則作限制性運(yùn)動(dòng)，系統(tǒng)至多為3自由度，故不可積，但是仍然可以找到一個(gè)穩(wěn)定的積分.當(dāng)主天體作圓運(yùn)動(dòng)時(shí)，第三體的運(yùn)動(dòng)方程可記為：

(4)

(5)

(6)

該常數(shù)是旋轉(zhuǎn)坐標(biāo)系中的能量，將其等同于能量積分，應(yīng)用能量穩(wěn)定化方法，可以得到如下形式的改正因子：

(7)

3 結(jié)果分析

如前所述，由于模型不可積，只能用數(shù)值方法求解.本文以四階龍格庫(kù)塔方法(RK4)為基本的數(shù)值算法，根據(jù)程序設(shè)計(jì)原則，能量穩(wěn)定化程序應(yīng)加在基本算法之后.為考察穩(wěn)定化程序的作用，首先觀察雅克比積分隨時(shí)間的變化規(guī)律，如圖1所示.從理論上分析雅克比積分為常數(shù)，在積分運(yùn)動(dòng)方程時(shí)應(yīng)該保持不變，但是RK4所得到的誤差結(jié)果卻是直線(xiàn)增加，也就是說(shuō)誤差隨時(shí)間逐漸積累，即得不到有效控制，由此以來(lái)積分時(shí)間越長(zhǎng)，數(shù)值穩(wěn)定性越差.從計(jì)算精度上看，RK4方法只能獲得不超過(guò)10-6的結(jié)果，與理想值有一定差距.然而加入穩(wěn)定化程序之后，結(jié)論大不相同，無(wú)論是穩(wěn)定性還是計(jì)算精度都獲得了令人滿(mǎn)意的結(jié)果.容易得出，能量穩(wěn)定化方法在圓形限制性三體問(wèn)題中也有重要的應(yīng)用價(jià)值.

圖2即為龐加萊截面圖，此時(shí)將RK4方法和穩(wěn)定化程序作對(duì)比，可以很清晰的看出傳統(tǒng)低階數(shù)值方法獲得的相空間結(jié)構(gòu)并不真實(shí)，截面顯示該軌道是非周期的，也就是說(shuō)該條軌道可能是混沌的，該結(jié)果與實(shí)際情況相差甚大，加入穩(wěn)定化可以得到真實(shí)情況，軌道是有三個(gè)島鏈組成的擬周期軌道，沒(méi)有混沌行為，用RK4得到的混沌現(xiàn)象只是假象.

圖2 龐加萊截面圖

為驗(yàn)證上述結(jié)論是否正確，本文還用高階算法做出功率譜和偏差軌道的結(jié)構(gòu)演化圖，相鄰軌道的坐標(biāo)初始值改變量為0.001，如圖3和圖4所示.兩種動(dòng)力學(xué)分析方法均表明該軌道是周期的.

圖3 功率譜圖 圖4 偏差軌道演化圖

4 結(jié)論

基于平面圓形限制性三體問(wèn)題模型存在孤立的不變積分常數(shù)即雅克比積分，本文將n體問(wèn)題中廣泛應(yīng)用的能量穩(wěn)定化方法應(yīng)用到該模型，數(shù)值實(shí)驗(yàn)得出能量穩(wěn)定化方法對(duì)雅克比積分保持了較高的數(shù)值穩(wěn)定性和精度，而且該方法能有效的避免數(shù)值誤差產(chǎn)生的假混沌現(xiàn)象，與高階精度算法得出結(jié)論一致，從而提供了更為簡(jiǎn)單有效的討論高維動(dòng)力系統(tǒng)非線(xiàn)性行為的工具.

[1] Feng K,Sym.Diff.Geometry & Diff. Equation[M].Beijing:Science press,1985.

[2] Wu X, Huang T Y, Wan X S,et al.Comparison among correction methods of individual Kepler energies in n-body simulations[J].A J,2007,133(6):2643-2653.

[3] Ma D Z, Wu X, Zhu J F.Velocity scaling method to correct individual Kepler energies[J].New Astron,2008,13(5):216-223.

[4] Nacozy P E. The use of Integrals in Numerical Integrations of the n-body problem[J].Ap&SS,1971,14(11):40-51.

[5] Simo C,Stuchi T J.Central stable/unstable manifolds and the destruction of KAM tori in the planar hill problem[D].Physica D,2000,140(6):1-32.