水星的運動規(guī)律摘要本文主要在已知水星的遠日點和繞日運行的線速度的條件下，通過建立微分方程模型，使用解析法和數(shù)值方法求解水星的軌道方程與位置。解析法的求解的過程中，結(jié)合了開普勒三大定律，準(zhǔn)確的給出了微分方程的精確解，求得水星到太陽的最近距離rx4.6016x1010(m)，水星繞太陽運行的周期約為88天。數(shù)值計算求解水星自m遠日點運行50天后的位置時，本文分別采用了Simpson求積法，基于壓縮映射的求根方法以及經(jīng)典的四階龍格一庫塔法，使用matlab數(shù)學(xué)軟件編程，得到了較為合理的行星運行模型的近似解，三種方法所得結(jié)果對應(yīng)分［=3.791，1x4.767x1010，92=3.791，丫?x4.767x1010及q=3.802，r3x4.779x101。。關(guān)鍵詞行星軌道微分方程Simpson法 四階龍格一庫塔法 matlab問題重述水星到太陽的最遠距離為0.6982x1011m，此時水星繞太陽運行的線速度為3.886x104m/$。試求問題一水星到太陽的最近距離問題二水星繞太陽運行的周期問題三從遠日點開始的第50天(地球天)結(jié)束時水星的位置并畫出軌道曲線問題分析求水星到太陽的最近距離以及水星繞太陽運行的周期等，需要先將水星軌道方程求出，因此可以根據(jù)Newton第二定律及萬有引力定律-mMGe9=m絲，建立微r2 dt2APAGEA1A分方程模型，將原問題轉(zhuǎn)化為求解帶有初值條件的微分方程問題，進而采用解析法或數(shù)值方法求解遠日點和周期。模型假設(shè).水星運行的軌道是以太陽為一個焦點的橢圓.從太陽指向水星的線段在單位時間內(nèi)掃過的面積相等.水星運行周期的平方與其運行軌道橢圓長軸的立方之比為常量符號系統(tǒng)v0 水星在遠日點的線速度M 太陽的質(zhì)量m 水星的質(zhì)量ro 水星在遠日點的距離T 周期五.建立模型與求解模型一水星的軌跡方程設(shè)太陽中心所在的位置為復(fù)平面的原點O,在時刻t,水星位于Z(t)=r&?所表示的點P。這里r=r(t),9=9(t)均為t的函數(shù)，分別表示Z(t)的模和輻角。于是TOC\o"1-5"\h\zdZdr d9 dr d9水星的速度為吆二空09+H空=e9(也+ir吧),加速度為dtdt dt dt dtd2Z 八 /d2r /d9、、,z d29 ddr 99、 /一、 口小,一日 上6+士口1=e9(--—r(-)2)+i(r——+2)(1.1),而太陽對仃星的引力依萬有引dt2 1dt2 dtdt2dtdt_力定律，大小為mMG，方向由行星位置P指向太陽的中心O,故為-mMGe9，其r2 r2中M=1.989x1030(kg)為太陽的質(zhì)量，m為水星的質(zhì)量，G=6.672x10-n(N?m2/kg2)為萬有引力常數(shù)。APAGEA1A

TOC\o"1-5"\h\z依Newton定律，我們得到-mMGe0=mdZ(1.2)，將(1.1)代入(1.2)，

r2 dt2然后比較實部與虛部，就有d20ddrd0八r——+2 =0dt2 dtdt1也-r(d0)2=-MG、dt2 dt r2需要加上定解距原點。為r°，這是兩個未知函數(shù)的二階微分方程組。在確定某一行星軌道時，條件。假設(shè)當(dāng)t=0時，行星正處于遠日點，而遠日點位于正實軸上，需要加上定解距原點。為r°，01dr<——dtt=dtt=0d0dtV

=-0-

t=0 r0d2rdt2rd2rdt2rIt=00L=0dr1dtt=0d0dtV

=—0-

t=0 r0因此問題轉(zhuǎn)化為求解帶初值問題的微分方程組TOC\o"1-5"\h\zd20 ddr00 八r——+2 =0dt2 dtdt/0、 MG一r(——)2= dtr2=r

0=0d0r2——=cdt1又將r"+2包膽=0兩邊同乘以r，即得—(r2d0r2——=cdt1其中c=rv，這樣有向線段加在時間At內(nèi)掃過的面積等于t1At1r2史dt=3，1 00 2dt2t這個正是Kepler的第二定律，從太陽指向水星的線段在單位時間內(nèi)掃過的面積相等。APAGEA1A

TOC\o"1-5"\h\zd2r ,de、 MG d2r c2 MG將(1.3)代入———r(一)2=-一得一-+=-一，于是我們可以得到水dt2dtr2 dt2 r3r2星運行的較為簡單形式的數(shù)學(xué)模型：d2rc2 MG ——i—=- dt2r3r2d0 c =—1dtr2，二o

dr.為了求得行星的軌跡方程，要消去變量t為了求得行星的軌跡方程，要消去變量t，令r=1，那么也=:可以改寫為

u dtrd0 drd,du、 d0 drd,du、 d2ud0 d2u——=cu2從而一=-c—(―—)=-c- =-c2u2——dt1 dt1dtd0 1d02dt1d02將上式代入d2rc2MG ——i—=- dt2 r3 r2化簡后為黑+u=P其中p=M，引進u=u-1，立即可以求出u--即可以求出u--=u=Acos(0-0)，p 0這里A和00是待定的常數(shù)。記e=Ap，上式可以寫為r= 1-ecos(0-0)這個就是水星的軌道方程，是一條平面二次曲線。由于水星繞太陽運行，故必有0<e<1。由于r在t=0時取道最大值r)(遠日點)，這個就意味著此時函數(shù)TOC\o"1-5"\h\zcos(0-0)取道最大值1.于是就有0=0,e=1-P，從而軌跡方程為0 0 r0r=一P一。對于水星而言，r=0.698x10n(m),v=3.886x104(m/s)，又水星的1-eco0s 0 0近日點到太陽的距離r=P=P。依據(jù)已知數(shù)據(jù)，可知m1-ecos冗 1+eAPAGEA1A

c乙c=rve2.713X1015(m2/s),p=缶”5.547x10io(m),e=1_p氏0.2055，從而計i00 MG r0算水星到太陽的最近距離為re4.6016x1010(m)模型二水星的運行周期設(shè)水星的周期為T,那么利用Kepler第二定律，我們有』”r2d°dt=1CT(1.4)02dt21上式左端為水星軌跡橢圓所圍的面積，記為S，由于橢圓的半長軸a=Jp-，半短1_e2軸b軸b=P，從而有S=^ab=—— V1_e2 3(1_e2)2將上式代入式(1.4),解得T= %P23C(1_e2)21關(guān)數(shù)據(jù)代入，易得Te7.6025x106(s)e87.9919(d)模型三水星的位置由于水星的運行滿足Kepler第二定律，則該式可改寫為fe+Aer2d°=CAt，從°而可得j°―p2ade=t

0C(1_ecose)2CTFCTF(e)=—1P2如果我們要求t=T1時相應(yīng)的e和r，則意味著首先要解方程其中F(e)=』e7—1 de0(1_ecose)2在求出了t=T時的e=e后，立即可以由r=-P—得到相應(yīng)的r。

1 1-ecose下面用數(shù)值方法求解水星的位置.Simpson法APAGEA1A.、一 ，, 1 、，_ CT…一由被積函數(shù)——1——的恒正性可知F(6)單調(diào)，從而方程F(9)=J的根必

(1-ecos6)2 p2存在且唯一一。取存=h,6存在且唯一一。取存=h,6=kh(k=1,2,...),

kt己F=F(6)。若F<kkCT尸CT—ii,F >—ii,P2 n+1 p2那么9位于6與6+1之間，在h適當(dāng)小時，可取6。6n。計算F(6)可采用不同的數(shù)值積分法，本文采用Simpson法，取步長h=0.001,具體求解過程見附錄一，最后結(jié)果為6=3.791,r氏4.767x1010.基于壓縮映像的求根方法我們引入水星軌道橢圓的參數(shù)方程,由于橢圓的半長軸1J,2，半短軸b=-JL=，從而中心到焦點的距離為db?=ae。因左焦點為原點，故橢圓中v1—e2心位于(ae,0),于是得到參數(shù)方程(x=a(e+cos①)[y=bsin中它們與r,6的關(guān)系為x2+y2=r2,—=tan6x此式可改寫成 CAt=』中+△"xy'-yx')d3=ab[esin(p+Ap)-esin①]+abApp當(dāng)t=T1時解方程.一一CTesinp+p=—ab記九=￡i,g(p)=入-esinp，那么上式即p=g(p)，就是說要去求函ab數(shù)g(p)的不動點，求解方程不動點可以采用簡單迭代法，對于水星，我們已計算出e“0.2055，由于e很小，因此迭代收斂理論上可以很快，當(dāng)時間從遠日點開始的第50天結(jié)束時，意味著T=0.432x107(s)，從而CTCT3X= =—M1-e2)2x3.5703abp2APAGEA1A

不妨取①=0，于是q二九-esin30二3.57033二九一esinb=3.65573二九-esin3=3.67473二九-esin3=3.6747故6=3.6747由式x=a(e+cos6),y=bsin6,x2+y2=r2,'=tan0，可以計算出相應(yīng)的0,x即由即由tan0= bsin__=tan0= bsin__=0.758490=0.64891,而0=0+兀=3.791此時的距離r為r=[[a(e+cos[)]2+[bsin共]2=4.7668x1010(m)3.經(jīng)典四階Runge-Kutte法由我們將由最初的微分方程組求解水星的位置，方程組見下d2rdt2d0MGdt r2r=rt=0 0drdtt=00ldtt=00lt=0令q=上，那么我們可以得到一階微分方程組:dtAPAGEA1Adq C2MG =-3—— dt r3r2dr~r=qdtTOC\o"1-5"\h\zd6 C<一二-dt r2rL=r0包=0dtt=061to = 0若記這個微分方程組中方程的右端依次為Q(t,q,r,6),R(t,q,r,6)和S(t,q,r,6),則相應(yīng)的四階Runge-Kutte迭代格式法為h. 、q=q+-(K+2K+2K+K)k+1 k6 1 2 3 4hr=r+-(L+2L+2L+L)k+1 k61 2 3 46k6k+i=6+-(N+2N+2N+N)k61 2 3 4h這里對于q =q+-(K+2K+2K+K)，有k+1 k6 1 2 3 4K1=Q(tk，qk，rk，Sk)hL+hL+ 1,62K=Q(t+,q+ 1,rk2k2k=Q=Q(tk+—,q+ 3-,r+2,6+ 3-)2k2k2k2K=Q(t+h,q+hK,r+hL,6+hN)初值為q0=0,r0=r0,6°=0，則對于給定的步長值h，類似可以逐步計算一系列的qk,r,之，由于行星繞著太陽運行，只需取6,42兀,6n+1>2兀，而取得行星軌道上一系列點的近似坐標(biāo)(r,6k)，再通過極坐標(biāo)與直角坐標(biāo)的轉(zhuǎn)換，繼而可以繪出軌道曲線。通過matlab編程求解得6=3.802，r氏4.779x1010,軌道曲線如下APAGEA1A

水星繞太陽運行的軌道曲線珀遠日點運行印天后的位置程序見附錄二。水星繞太陽運行的軌道曲線珀遠日點運行印天后的位置程序見附錄二。六. 模型推廣本文建立的微分方程模型對于求解行星繞日運行軌道具有廣泛的應(yīng)用空間，只需給出行星的遠日點和在遠日點的運行線速度即可計算出軌道方程，用數(shù)學(xué)軟件繪出近似的軌道曲線，對于研究天體運行有所幫助。此外，本文采用的求解微分方程的數(shù)值方法，具有較為快速且準(zhǔn)確的收斂效果，可以用來求解其他類似的微分方程模型。七. 參考文獻【1】樂經(jīng)良，數(shù)學(xué)實驗，北京，高等教育出版社，1999年10月【2】周品，matlab數(shù)值分析，北京，機械工業(yè)出版社，2009年1月APAGEA1A八.附錄附錄一functionq1=y2(x)q1=(1-0.2055*cos(x))「-2;h=0.001;k=1;x=h*k;f=quad('y2',0,x)whilef(k)<3.8091k=k+1;x=k*h;f(k)=quad('y2',0,x);endx附錄二formatlongc1=2.7132e15;M=1.989e30;G=6.672e-11;Q=inline(，2.7132e15"2/(r"3)-1.989e30*6.672e-11/(r"2)');R=inline(，q，);S=inline('2.7132e15/(r"2)');q=0;r=0.6982e11;theta=0;t=0;k=1;h=0.001e7;whiletheta<=2*piK1=Q(r);L1=R(q);N1=S(r);K2=Q(r+h/2*L1);L2=R(q+h/2*K1);N2=S(r+h/2*L1);K3=Q(r+h/2*L2);L3=R(q+h/2*K2);N3=S(r+h/2*L2);K4=Q(r+h*L3);L4=R(q+h*K3);N4=S(r+h*L3);t=t+h;APAGEA1Aq=q+h/6*(K1+2*K2+2*K3+K4);r=r+h/6*(L1+2*L2+2*L3+L4);theta=theta+h/6*(N1+2*N2+2*N3+N4);rr(k)=r;ee(k)=theta;xx(k)=rr(k)*cos(ee(k));%水星任意位置的橫坐標(biāo)yy(k)=rr(k)*sin(ee(k));%水星任意位置的縱坐標(biāo)k=k+1;end;plot(xx,yy)%畫出水星的軌道曲線text(0,0,‘太陽')text(0.6982e11,0,'遠日點')text(-4.6078e+010,0,'近日點')holdon;plot(0,0,'r.','MarkerSi