水星的运动规律

1、水星的运动规律摘要本文主要在已知水星的远日点和绕日运行的线速度的条件下，通过建立微分方程模型，使用解析法和数值方法求解水星的轨道方程与位置。解析法的求解的过程中，结合了开普勒三大定律， 准确的给出了微分方程的精确解， 求得水星到太阳的最近距离 rm定 4.6016x101(m),水星绕太阳运行的周期约为 88 天。数值计算求解水星自远日点运行50 天后的位置时，本文分别采用了 Simpson 求积法，基于压缩映射的求根方法以及经典的四阶龙格一库塔法，使用 matlab 数学软件编程，得到了较为合理的行星运行模型的近似解，三种方法所得结果对应分用=3.791,1球 4.767 父 1010,也=

2、3.791,r2归 4.767 父 1010及司=3.802,r3划 4379M1010。关键词行星轨道微分方程 Simpson 法四阶龙格一库塔法 matlab一.问题重述水星到太阳的最远距离为0.6982父1011m 此时水星绕太阳运行的线速度为3.886父104m/s。试求问题一水星到太阳的最近距离问题二水星绕太阳运行的周期问题三从远日点开始的第 50天(地球大)结束时水星的位置并画出轨道曲线二.问题分析求水星到太阳的最近距离以及水星绕太阳运行的周期等，需要先将水星轨道方程2r求出，因此可以根据 Newton 第二定律及万有引力定律-嘤生七二m二|，建立微rdt分方程模型，将原问题转化为

3、求解带有初值条件的微分方程问题，进而采用解析法或数值方法求解远日点和周期。三.模型假设1.水星运行的轨道是以太阳为一个焦点的椭圆2 .从太阳指向水星的线段在单位时间内扫过的面积相等3 .水星运行周期的平方与其运行轨道椭圆长轴的立方之比为常量四.符号系统1.Vo水星在远日点的线速度2.M太阳的质量3.m 水星的质量4.ro水星在远日点的距离5.T周期五.建立模型与求解模型一水星的轨迹方程设太阳中心所在的位置为复平面的原点 O,在时刻 t,水星位于 Z(t)=rei0所表示的点 P。这里r=r(t),6=t)均为 t 的函数，分别表示Z(t)的模和辐角。于是水星的速度为乌=曳 89+舟9电一双曳+

4、ir吧),加速度为dtdt出dt出咚=9(曰-r(西)2)+i(r曰+2支妙J(1.1),而太阳对行星的引力依万有引dt:dtdtdtdtdt_力定律，大小为mMG,方向由行星位置 P 指向太阳的中心 O,故为-mMGe，其rr中 M=1.989M1030(kg)为太阳的质量， m 为水星的质量， G=6.672101(Nm2/kg2)为万有引力常数。行星的速度为 v0。那么就有初值条件:drdt由dt因此问题转化为求解带初值问题的微分方程组ro_d日drd8rr+2dt2dtdtir,er。8 匕=0_V0t-0一0又将rd|+2如此=0两边同乘以 r,即得色“2吧）=0,从而r23=c,d

5、t2dtdtdtdtdttt1cdic.t其中 C1=rv,这样有向线段 OP 在时间事内扫过的面积等于 1r2tdt=E,t2dt2依 Newton 定律，我们得到一吗。rd27=md4(1.2),将 1.1)代入(1.2),dt然后比较实部与虚部，就有dSr-dtd2r羽cdrd?,2 二 0dtdt,du、2MG-r()=dtr2这是两个未知函数的二阶微分方程组。条件。假设当 t=0 时，行星正处于远日点,在确定某一行星轨道时，需要加上定解而远日点位于正实轴上，距原点。为 r0,t=0tz0二r0二0tZ0=0_V0t=0=一d2rdt2-注)dtdrdtt=0=0史dt这个正是 Kep

6、ler 的第二定律，从太阳指向水星的线段在单位时间内扫过的面积相等。,222将（1.3）代入d4r（变）2=-学得町-c3=-嘤，于是我们可以得到水dt2出r2dt2r3r2星运行的较为简单形式的数学模型：22drG2 一dtrddtdrd,du、d2ud122d2u二C|()=-C1TV=-C1u77dtdtd?ddtd12/2,与，化简后为 d4+u=-(1.4),其中 p=u,引进 u=u-1，立r2d12pMGp即可以求出u-，=u=Acos（8-日0）,这里 A 和日。是待定的常数P以写为p1-ecos(1-2)这个就是水星的轨道方程，是一条平面二次曲线。由于水星绕太阳运行，故必有0

7、e1。由于 r 在 t=0 时取道最大值 r0（远日点），这个就意味着此时函数cos。-备）取道最大值 1.于是就有为=0,e=1-2，从而轨迹方程为r0r=p0对于水星而言，r=0.698M1011（m）,v。=3.886M104（m/s）,又水星的1-ecos近日点到太阳的距离rm=p=p。依据已知数据，可知1-ecos二1eMGdrdtt=o=0为了求得行星的轨迹方程，要消去变量t,令 r=1,那么变=与可以改写为udtrd2一二c1udt,22drG2dtr将上式代入记e=Ap,上式可152G210G=r0v0ft2.713x10(m2/s),p5.547=10(m),MG模型二水星的

8、运行周期(1.5)将有关数据代入，易得 T：7.6025106(s):87.9919(d)模型三水星的位置由于水星的运行满足 Kepler 第二定律，则该式可改写为 Er2d8=CQt,从而可得p2dl-t)C1(1-ecos?)如果我们要求 t=TI时相应的日和 r,则意味着首先要解方程，FR)=CT1,p其中在求出了 t=T1 时的日后，立即可以由 r=-J 得到相应的 r。1 -ecos-卜面用数值方法求解水星的位置1 .Simpson 法e=1-fc0.2055,从而计算水星到太阳的最近距离为rm4.60161010(m)设水星的周期为Kepler第二定律，我们有T1r2d%tJCIT

9、02dt2(1.4)上式左端为水星轨迹椭圆所围的面积，记为 S,由于椭圆的半长轴a，半短1一e2轴 b=,p,从而有 S=nab=3231e(1-e2)2将上式代入式(1.4),解得 T=27Tp23C1(1-e2)2C1T1CIT,2abp3(1-e2)23.5703由被积函数12 的恒正性可知F(e)单调， 从而方程 F(e)=C3 的根必(1-ecosi)2p存在且唯一。取 Ae=h,ek=kh(k=1,2,.),记 Fk=F(%)。若Fn3 理，七书之理，PP那么 3 位于斗与斗卅之间，在 h 适当小时，可取 6&en0nn-in计算F(6)可采用不同的数值积分法，本文采用 Simps

10、on 法，取步长 h=0.001,具体求解过程见附录一，最后结果为三=3.791,r，4.76710102 .基于压缩映像的求根方法我们引入水星轨道椭圆的参数方程，由于椭圆的半长轴a=P2,半短轴1-eb=,从而中心到焦点的距离为 la2-b2=ae。因左焦点为原点，故椭圆中1 -e2心位于(ae,0),于是得到参数方程,Lx;a(ecos)y-bsin它们与r,e的关系为x2y2=r2,-=tanx一.、_W44 一“4此式可改与成 C/t=(xy-yx)d=abesin(：小)-esinabG当 t=TI时解方程esin-CarCT记九=C-,g*=esin、那么上式即邛=gF),就是说要

11、去求函 ab数gT)的不动点，求解方程不动点可以采用简单迭代法，对于水星，我们已计算出e*0.2055,由于 e 很小，因此迭代收敛理论上可以很快，当时间从远日点开始的第 50 天结束时，意味着 T1=0.432 父 107(s),从而不妨取中0=0,于是1=-esin0=3.57032-esin1=3.65575=-esin4=3.67476=-esin5=3.6747故一二 3.6747由式x=a(e+cosy=bsin中，x2+y2=r2,y=tan6,可以计算出相应的 0,x即由一 bsin:tarn=0.75849a(ecos)得 8=0.64891,而 6=8+n=3.791此时的

12、距离 r 为r=Ja(e+cos)2十bsin/2=4.7668父1010(m3.经典四阶 Runge-Kutte 法由我们将由最初的微分方程组求解水星的位置，方程组见下亚上 MGdt2-r3一一 r2dc,dt-r2drt田=0drI 八出t=0=0令qd,那么我们可以得到一阶微分方程组:dt_2_dqCiMGdrqdtd-Cidtr2t_0drdtt=0日,t卫二0若记这个微分方程组中方程的右端依次为 Q(t,q,r,0),R(t,q,r,8)和 S(t,q,r,6),则相应的四阶 Runge-Kutte 迭代格式法为h,qk1=qk(KI2K22K3K4)6和久2h八1一人一(N12叫2

13、N3N4)6h,这里对于qk=qk+(K1+2K2+2K3+K4),有6初值为 q。=0,r。=0,=0,则对于给定的步长值 h,类似可以逐步计算一系列的 qkJk,由于行星绕着太阳运行，只需取W2 兀,5 书之 2 兀，而取得行星轨道系列点的近似坐标(rk,ek),再通过极坐标与直角坐标的转换，继而可以绘出轨道曲线。通过 matlab 编程求解得日=3.802,r比4.779父101,轨道曲线如下KI=Q(tk,qk,r-k)K2hhKi=Q(tk-,qk-/k22.hLi-T,A.吗2K4h二Q(tk-,qkhK22，rk2)-Q(tkh,qkhK3,rk也,，h)程序见附录二。六.模型推

14、广本文建立的微分方程模型对于求解行星绕日运行轨道具有广泛的应用空间，只需给出行星的远日点和在远日点的运行线速度即可计算出轨道方程，用数学软件绘出近似的轨道曲线，对于研究天体运行有所帮助。止匕外，本文采用的求解微分方程的数值方法，具有较为快速且准确的收敛效果，可以用来求解其他类似的微分方程模型。七.参考文献11 乐经良，数学实验，北京，高等教育出版社，1999 年 10 月2周品，matlab 数值分析，北京，机械工业出版社，2009 年 1 月附录附录一functionq1=y2(x)q1=(1-0.2055*cos(x)A-2;h=0.001;k=1;x=h*k;f=quad(y2,0,x)

15、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.7132e15A2/(rA3)-1.989e30*6.672e-11/(rA2);R=inline(q);S=inline(2.7132e15/(rA2);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

16、=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;q=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);%水星任意位置的横坐标yy(k)=rr(k)*sin(ee(k);%水星任意位置的纵坐标k=k+1;end;plot(xx,yy)%画

17、出水星的轨道曲线text(0,0,太阳)text(0.6982e11,0,远日点)text(-4.6078e+010,0,近日点)holdon;plot(0,0,r.,MarkerSize,20);holdoffholdon;plot(0.6982e11,0,r，MarkerSize,20);holdoffholdon;plot(-4.6078e+010,0,r.,MarkerSize,20);holdofftitle(水星绕太阳运行的轨道曲线)clcq=0;r=0.6982e11;theta=0;t=0;k=1;h=0.001e7;whilet=50*24*3600%求水星自远日点开始第 50 天的位置K1=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;q=q+h/6