郝海玲

(晉中職業(yè)技術學院　基礎部，山西省　晉中市　030600)

?

求解兩體問題數值方法的創(chuàng)新性研究

郝海玲

(晉中職業(yè)技術學院基礎部，山西省晉中市030600)

該文闡述了求解兩體問題的線性對稱多步數值方法——Obrechkoff法。N體問題是一個很難的問題，只有少數微分方程存在解析解，近似方法是解微分方程的主要手段，高精度的軌道問題需要長時間的數值積分。因此，選用線性對稱多步方法，在其主結構上增加高階微商，不僅有理想的精度和較好的穩(wěn)定性，而且可以大大減少截斷誤差，尤其在很大程度上減小了誤差系數。研究表明，該方法求解兩體問題的數值解具有高精度、高效率及穩(wěn)定性好的優(yōu)點。

P穩(wěn)定；兩體問題；高階微商；截斷誤差

N體問題是一個很難的問題，高精度的軌道問題需要長時間的數值積分。然而許多傳統(tǒng)的數值方法，如Runge-Kutta法、Adams法、Euler法和St?rmer-Cowell法等都不是很適用，因為這些方法都存在許多缺點。如它們的穩(wěn)定區(qū)域很小，長時間的積分會使誤差不斷增加，精度也很難通過增加節(jié)點來提高。因此需要一種高精度、穩(wěn)定性好的方法來解兩體問題。三角擬合的對稱多步法已經解決了很多問題[1-5]，所以采用Obrechkoff[6]多步方法來研究兩體問題。其利用高階微商來減少截斷誤差，很大程度上減小了誤差系數[7]。

1　數值方法介紹

Obrechkoff方法是利用增加高階微商來提高方法的精度。該方法在1942年被提出，但由于高階微商的使用計算非常復雜，并且當時的計算機有所限制，所以沒有得到廣泛的應用。近年來由于計算機的計算功能不斷提高，Obrechkoff方法已經被廣泛應用。如薛定鄂方程或Duffing方程。這兩種方法對線性多步方法增加高階微商，高階微商可以表示成一階微商和二階微商的線性組合。然而在計算過程中發(fā)現，一階微商的公式的精度也影響著該數值方法的精度[8]。Obrechkoff方法的計算方程為：

(x-ih))+…

(1)

為了保證方法的穩(wěn)定性，該方法的系數仍然通過三角擬合的方法。

O(h15)

(2)

可以看出四步方法的誤差系數非常小，并且具有較高的精度。

O(h17)

(3)

(4)

通過比較截斷誤差的誤差系數和h的階數，發(fā)現增加高階微商能夠明顯提高數值方法的精度，因此最后選擇用Obrechkoff數值方法來計算兩體問題。

2　Obrechkoff方法解兩體問題

本文用Obrechkoff方法求解著名的兩體問題[9]的數值解，為計算簡便，其中一個質點的質量遠遠大于另一個質點，所以小質點的質量可以忽略不計，大質點的質量等于1，引力常數為1，根據萬有引力定律，小質點的運動方程為：

-f(x(t),y(t))x(t)

(5)

-f(x(t),y(t))y(t)

(6)

顯然式(5)、式(6)是非線性常微分方程。兩體問題的解析解為：

(7)

其中，u和t的關系滿足

u-ecos(u)-t=0

(8)

式中，u表示做橢圓運動的轉動的角度，t表示時間，e表示偏心率(0

本文嘗試用Obrechkoff方法的各種不同形式對兩體問題進行了數值計算，希望能夠提高數值解的精度，下面介紹多次計算，結果都比較理想。

P穩(wěn)定四步方法的主要結構如下：

q0(h)=y(t+2h)+y(t-2h)+c0y(t)+

c1(y(t+h)+y(t-h))+h2[c2(y″(t+2h)+

y″(t-2h))+c3(y″(t+h)+y″(t-h))+

c4y″(t)]+h4[c5(y(4)(t+h)+y(4)(t-h))+

c6y(4)(t)]

(10)

圖1　ω=1,e=0.1,h=0.1,t=100 000

本文還對四步高階微商方法加上最外點，需要利用Newton線性化方法，雖然減少了誤差，提高了精度，但是計算工作量非常復雜，對現在的計算機內存來說計算量過大，因此除了期待更加先進的計算機問世之外，我們還需要尋找更簡便且精度較高的數值方法。因此，本文還運用了兩步四階微商對稱方法進行了數值計算，得到的結果是最理想的，如圖2所示。

圖2　ω=1,e=0.1,h=0.1,t=150 000

而用八步線性對稱方法得到了兩體問題的數值解[3]，如圖3所示。

圖3　ω=1,e=0.1,h=0.1,t=100 000

通過對圖1～圖3進行比較，可以看出新的方法不僅精度更高且穩(wěn)定性較好。在計算過程中發(fā)現，對于兩體問題，當e<0.5時，該方法比較理想。當1>e>0.5時，該方法的誤差就會隨著時間的增加而增大。但是這種現象不是數值方法引起的，可以由它的物理意義，即開普勒第二定律來解釋：行星和太陽之間的連線在相等的時間內所掃過的面積相等。質點的運動方程表明坐標中心取在橢圓焦點上，如圖4所示。

圖4　不同偏心率的橢圓

當經過相等的時間，e值大的橢圓較扁，并且右焦點離質點運動軌道很近，而e值小的橢圓，質點的軌道離焦點較遠。因此，在相等的時間內，各個橢圓右焦點附近的質點的速度有很大的不同，運動質點在e值大的橢圓上速度較快，在e值小的橢圓上速度較慢，這樣就表現為在e值大的軌道上運動的誤差隨時間的加長而有可能增加得更快。

3　結束語

這種新的P穩(wěn)定Obrechkoff結構求解兩體問題，目前達到了理想的精度和較好的穩(wěn)定性，但還需進一步研究。若在穩(wěn)定性良好且縮小步長的情況下能提高精度，那么該方法將可以運用到三體問題，從而為研究N體問題提供精度較高的數值方法，能夠比較精確地了解天體中的某個星球經過長時間的運動后的運動情況，預知該星球是否會影響人們賴以生存的地球。

[1]WANGZhongcheng.Anewtrigonometrically-fittingtechniquetoconstructasymmetriclinearmulti-stepmethodforthenumericalsolutionofanorbitalproblem[J].NewAstronomy，2005, 11(2)：90-102.

[2]SIMOSTE.Exponentially-fittedandtrigonometrically-fittedmethodsforlong-termintegrationoforbitalproblems[J].NewAstronomy，2003, 8(5)：391-400.

[3]SIMOSTE.Exponentially-fittedandtrigonometrically-fittedmethodsforthenumericalsolutionoforbitalproblems[J].NewAstronomy，2003(8);391-400.

[4]SIMOSTE.Trigonometrically-fittedpartitionedmultistepmethodsfortheintegrationoforbitalproblems[J].NewAstronomy，2004, 9(6)：409-415.

[5]SIMOSTE.Dissipativetrigonometricallyfittedmethodsforthenumericalsolutionoforbitalproblems[J].NewAstronomy，2004, 9(1)：59-68.

[6]NAUKA.Obrechkoffn,onmechanicalquadrature(bulgarianfrenchsummary)[J].SpisanieBulgar，1942(65)：191-289.

[7]ZHAODeyin，WANGZhongcheng，DAIYongming.Importanceofthefirst-orderderivativeformulaintheObrechkoffmethod[J].ComputerPhysicsCommunications，2005, 167(2)：65-75.

[8]SIMOSTE.Ap-stablecompleteinphaseObrechkofftrigonometricfittedmethodforperiodicinitial-valueproblem[J].ProceedingsoftheRoyalSocietyA，1993, 441(1912)：283-289.

[9]PSIHOYIOSG，SIMOSTE.Trigonometrically-fittedsymmetricmultistepmethodsfortheapproximatesolutionoforbitalproblems[J] .NewAstronomy，1939, 14(2)：175-184.

[10]WANGZhongcheng.Trigonometrically-fittedmethodwithfourierfrequencyspectrumforundampedduffingequation[J].ComputerPhysicsCommunications，2006，175(4)：241-249.

Innovative Research of Obrechkoff Method for the Numerical Solution of Orbital Problems

HAO Hailing

(DepartmentofFoundation,JinzhongVocationalTechnicalCollege,Jinzhong030600,China)

WefocusonthenewkindofP-stableObrechkoffmethodforthenumericalsolutionoforbitalproblems.However,onlyafewofthesedifferentialequationscanbesolvedexactly.Approximatemethodsarethemainmeansforsolving,analyzingandunderstandingphysicsproblems.ThroughimprovingtheWang’smethod,wedevelopanewkindofP-stablefour-stepObrechkoffmethodbyaddingthehigher-orderderivatives.Thisproposedmethodisveryeffectivebuthasveryhighlocaltruncationerror.ThenumericalexperimentsforthenumericalsolutionoforbitalproblemshastheadvantageovertheWang’smethodinaccuracyandefficiency.

P-stable；orbitalproblem；higher-orderderivatives；localtruncationerror

2014-10-11；修改日期: 2016-06-22

郝海玲(1979-)，女，碩士，講師，主要從事物理問題的數值解法方面的研究。

O411

Adoi:10.3969/j.issn.1672-4550.2016.04.014