水星的运动规律

水星的运动规律摘要本文主要在已知水星的远日点和绕日运行的线速度的条件下，通过建立微分方程 模型，使用解析法和数值方法求解水星的轨道方程与位置。解析法的求解的过程中，结合了开普勒三大定律，准确的给出了微分方程的精确解，求得水星到太阳的最近 距离rm 4.6016 1010(m)，水星绕太阳运行的周期约为88天。数值计算求解水星自 远日点运行50天后的位置时，本文分别采用了 Simpson求积法，基于压缩映射的 求根方法以及经典的四阶龙格一库塔法，使用matlab数学软件编程，得到了较为合理的行星运行模型的近似解，三种方法所得结果对应分-3.791，A 4.767 1010，10 1023.791，D 4.767 10 及 33.802，B 4.779 10。关键词行星轨道微分方程Simpson法四阶龙格一库塔法matlab问题重述水星到太阳的最远距离为0.6982 1011m此时水星绕太阳运行的线速度为3.886 104 m/ s。试求问题一水星到太阳的最近距离问题二水星绕太阳运行的周期问题三 从远日点开始的第50天(地球天)结束时水星的位置并画出轨道曲线问题分析求水星到太阳的最近距离以及水星绕太阳运行的周期等，需要先将水星轨道方程求出，因此可以根据Newton第二定律及万有引力定律mMG i2 erd2Zm dt2，建立微分方程模型，将原问题转化为求解带有初值条件的微分方程问题，进而采用解析法 或数值方法求解远日点和周期。模型假设1 水星运行的轨道是以太阳为一个焦点的椭圆2从太阳指向水星的线段在单位时间内扫过的面积相等3水星运行周期的平方与其运行轨道椭圆长轴的立方之比为常量四. 符号系统1 . Vo水星在远日点的线速度2. M太阳的质量3. m水星的质量4. r。水星在远日点的距离5. T周期五. 建立模型与求解模型一水星的轨迹方程设太阳中心所在的位置为复平面的原点0,在时刻t,水星位于Z(t) rei所表示的点P。这里r r(t),(t)均为t的函数，分别表示Z(t)的模和辐角。于是水星的速度为手評心鑫ei吟吟)，加速度为d2Zdt2ei (4dt2r(j) i(rd2dtdt22dr Jdt dt(1.1)，而太阳对行星的引力依万有引力定律，大小为mMG2 ， r方向由行星位置P指向太阳的中心O,故为mMG i 厂e r，其中M 1.989 1030(kg)为太阳的质量，m为水星的质量，G 6.672 10 11(N m2/kg2)为万有引力常数。依Newton定律，我们得到2mMG id Z2 e m 2(1.2)，将(1.1)代入(1.2)，rdt然后比较实部与虚部，就有d2,rd小r 220dtdt dtd2rzd、2MGdt2 r(dt)r2这是两个未知函数的二阶微分方程组。在确定某一行星轨道时，需要加上定解 条件。假设当t=0时，行星正处于远日点，而远日点位于正实轴上， 距原点0为ro， 行星的速度为vo。那么就有初值条件：r t 0ro0Vqr°drdtd dt因此问题转化为求解带初值问题的微分方程组dt22 dr d_dt dtd2rdt2/d、20MGrr t 0rodrdtddtVqd2又将 rr，即得其中& r°Vo，这样有向线段t内扫过的面积等于t t丄rJt t 2 dt个正是Kepler的第二定律，从太阳指向水星的线段在单位时间内扫过的面积相等。将（3）代入亨吟2啤得岁纟r dt rMGT ，r于是我们可以得到水星运行的较为简单形式的数学模型:.2 2d rg2-3dtrdqdtr2MGr2r t 0r0drdt为了求得行星的轨迹方程，要消去变量，那么ddt2可以改写为rd dtc1u2从而生沁（空dt dt dd2uC1 d 2 dt2 2 d u c1 ud2rdt22C13r2MG-，化简后为理rd1-(1.4)，其中 p P2C1MG，引进U1，立P即可以求出1-u Acos( 0)P，这里A和°是待定的常数。记e Ap，上式可以写为r1 ecos( 0)这个就是水星的轨道方程，是一条平面二次曲线。由于水星绕太阳运行，故必 由于r在t=0时取道最大值ro （远日点），这个就意味着此时函数cos( 0)取道最大值1.于是就有00,e1卫，从而轨迹方程为 r°pr1 ecos。对于水星而言，r0 0.6981110 (m), V043.886 10 (m/s)，又水星的近日点到太阳的距离rm1 ecos 1 e依据已知数据，可知Ci r°Vo 2.713 1015(m2/s) , p 5.547 1010(m) , e 1 - 0.2055，从而计 MGr0算水星到太阳的最近距离为rm 4.6016 1010(m)模型二水星的运行周期设水星的周期为T，那么利用 Kepler第二定律，我们有T 1 2 d1-r2dt GT( 1.4)0 2dt2上式左端为水星轨迹椭圆所围的面积，记为S,由于椭圆的半长轴a轴bp一-，从而有S胡e2ab 一(12L飞将上式代入式(1.4),解得T2 p_3C1 (1(1.5 )将有关数据代入,易得 T7.6025106(s)87.9919(d)模型三水星的位置由于水星的运行满足Kepler第二定律,则该式可改写为r2dC1t，从而可得0 C1(1 ecos )2如果我们要求t T1时相应的和r，则意味着首先要解方程，F(C1T1)p2其中1F() 0(ncor?d在求出了 t 时的后，立即可以由rp一得到相应的r1 ecos下面用数值方法求解水星的位置1. Simps on 法由被积函数(12的恒正性可知F()单调，从而方程ecos )F( )2PC1T1的根必存在且唯一。取h, k kh(k 1,2,.)， FkF( k)。若 FnC1T1f2 , Fn 1pC1T1T，P那么位于n与n 1之间，在h适当小时，可取计算F()可采用不同的数值积分法，本文采用Simps on法，取步长h=0.001.具体求解过程见附录一，最后结果为3.791，r4.767 10102.基于压缩映像的求根方法我们引入水星轨道椭圆的参数方程,由于椭圆的半长轴，半短轴b 一P-，从而中心到焦点的距离为1 e2ae。因左焦点为原点，故椭圆中心位于(ae，0)，于是得到参数方程a(e cos bsi n它们与r,的关系为r2tan x此式可改写成G(xy'yx')d abes in()esin ab当t T1时解方程esinC1T1abCi”，g(esin，那么上式即 g()，就是说要去求函数g()的不动点，求解方程不动点可以采用简单迭代法，对于水星，我们已计算出e 0.2055，由于e很小，因此迭代收敛理论上可以很快，当时间从远日点开始的第50天结束时，意味着£0.432 107(s),从而也孚(1 e2)弓 3.5703ab p即由不妨取0由式X0，于是1esi n03.57032esi n13.65575esin 43.67476esin 53.6747故_3.6747cos ),y bsi n2 2，x ya(er2,- tan ,x可以计算出相应的tanbsin0.75849a(e cos )0.64891，而3.791此时的距离r为r、a(e cos )2 bsin 24.76681010( m3.经典四阶Runge-Kutte法由我们将由最初的微分方程组求解水星的位置，方程组见下MGdc,dtr2r t o rodrdt令q，那么我们可以得到一阶微分方程组:dq dt dr dt drdrdtCi2r3Cl2rr°若记这个微分方程组中方程的右端依次为MG2rQ(t,q, r,),R(t,q,r,)和S(t,q,r,),则相应的四阶Runge-Kutte迭代格式法为qk ihqk(K162K2 2K3 K4)rk 1rk 7 (L162L2 2L3 L4)k 1hk 6(N12N2 2N3 N4)这里对于h 2qk iqk(Ki62K2 2K3K4），有K1Q(tk , qk ,rk, k)K2Q(tk 列hK1hL1k，rk，22hN.k 2 )K3hhK2hL2Q(tk 2，qk2，rk 2，hN2k 2 )K4Q(tk h,qkhK3,rk hL3, khZ)初值为 q°0, r°r°, °0，则对于给定的步长值h，类似可以逐步计算一系列的qk , rk , k ，由于行星绕着太阳运行，只需取n 2 , n 12 ，而取得行星轨道上一系列点的近似坐标（rk, k），再通过极坐标与直角坐标的转换，继而可以绘出轨道曲线。通过matlab编程求解得 3.802， r 4.779 1010，轨道曲线如下程序见附录二。六. 模型推广本文建立的微分方程模型对于求解行星绕日运行轨道具有广泛的应用空间，只 需给出行星的远日点和在远日点的运行线速度即可计算出轨道方程，用数学软件绘 出近似的轨道曲线，对于研究天体运行有所帮助。此外，本文采用的求解微分方程 的数值方法，具有较为快速且准确的收敛效果，可以用来求解其他类似的微分方程 模型。七. 参考文献【1】 乐经良，数学实验，北京，高等教育出版社，1999年10月【2】 周品，matlab数值分析，北京，机械工业出版社，2009年1月八.附录附录一function q1=y2(x)q1= (1-0.2055*cos(x)4-2;h=0.001;k=1;x=h*k;f=quad('y2',0,x)while f(k)<3.8091k=k+1;x=k*h;f(k)=quad('y2',0,x);endx附录二format longc1=2.7132e15;M=1.989e30;G=6.672e-11;Q=i nlin e('2.7132e15A2/(rA3)-1.989e30*6.672e-11/(rA2)');R=i nlin e('q');S=i nli ne('2.7132e15/(rA2)');q=0;r=0.6982e11;theta=0;t=0;k=1;h=0.001e7;while theta<=2*piK仁 Q(r);L1=R(q);N 仁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*(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)*si n(ee(k);%水星任意位置的纵坐标k=k+1;en d;plot(xx,yy)%画出水星的轨道曲线text(O,O,' 太阳')text(0.6982e11,0,' 远日点')text(-4.6078e+010,0,' 近日点')hold on ;plot(0,0,'r.','MarkerSize',20);hold offhold on ;plot(0.6982e11,0,'r.','MarkerSize',20);hold offhold on ;plot(-4.6078e+010,0,'r.','MarkerSize',20);hold offtitle('水星绕太阳运行的轨道曲线')clcq=0;r=0.6982e11;theta=0;t=0;k=1;h=0.001e7;while t<=50*24*3600% 求水星自远日点开始第50天的位置K仁 Q(r);L1=R(q);N 仁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*(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)*si n( ee(k);k=k+1;endr=rr(k)theta=theta% 水星自远日点开始第50天的位置(半径及角度)text(xx(k),yy(k),' 自远日点运行 50天后的位置')hold onplot(xx(k),yy(k),'g.','MarkerSize',20);hold off