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

郝海玲

(晉中職業(yè)技術(shù)學(xué)院　基礎(chǔ)部，山西省　晉中市　030600)

?

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

郝海玲

(晉中職業(yè)技術(shù)學(xué)院基礎(chǔ)部，山西省晉中市030600)

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

P穩(wěn)定；兩體問題；高階微商；截?cái)嗾`差

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

1　數(shù)值方法介紹

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

(x-ih))+…

(1)

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

O(h15)

(2)

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

O(h17)

(3)

(4)

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

2　Obrechkoff方法解兩體問題

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

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

(5)

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

(6)

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

(7)

其中，u和t的關(guān)系滿足

u-ecos(u)-t=0

(8)

式中，u表示做橢圓運(yùn)動(dòng)的轉(zhuǎn)動(dòng)的角度，t表示時(shí)間，e表示偏心率(0

本文嘗試用Obrechkoff方法的各種不同形式對(duì)兩體問題進(jìn)行了數(shù)值計(jì)算，希望能夠提高數(shù)值解的精度，下面介紹多次計(jì)算，結(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

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

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

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

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

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

圖4　不同偏心率的橢圓

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

3　結(jié)束語(yǔ)

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

[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