郝海玲

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

?

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

郝海玲

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

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

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

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

1　數(shù)值方法介紹

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

(x-ih))+…

(1)

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

O(h15)

(2)

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

O(h17)

(3)

(4)

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

2　Obrechkoff方法解兩體問題

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

-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表示做橢圓運動的轉(zhuǎn)動的角度，t表示時間，e表示偏心率(0

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

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

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線性化方法，雖然減少了誤差，提高了精度，但是計算工作量非常復雜，對現(xiàn)在的計算機內(nèi)存來說計算量過大，因此除了期待更加先進的計算機問世之外，我們還需要尋找更簡便且精度較高的數(shù)值方法。因此，本文還運用了兩步四階微商對稱方法進行了數(shù)值計算，得到的結(jié)果是最理想的，如圖2所示。

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

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

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

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

圖4　不同偏心率的橢圓

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

3　結(jié)束語

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

[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-)，女，碩士，講師，主要從事物理問題的數(shù)值解法方面的研究。

O411

Adoi:10.3969/j.issn.1672-4550.2016.04.014