水星的运动规律

PAGE1水星的运动规律摘要本文主要在已知水星的远日点和绕日运行的线速度的条件下，通过建立微分方程模型，使用解析法和数值方法求解水星的轨道方程与位置。解析法的求解的过程中，结合了开普勒三大定律，准确的给出了微分方程的精确解，求得水星到太阳的最近距离)(104.601610mrm，水星绕太阳运行的周期约为88天。数值计算求解水星自远日点运行50天后的位置时，本文分别采用了Simpson求积法，基于压缩映射的求根方法以及经典的四阶龙格—库塔法，使用matlab数学软件编程，得到了较为合理的行星运行模型的近似解，三种方法所得结果对应分13.791，1014.76710r，23.791，1024.76710r及33.802，1034.77910r。关键词行星轨道微分方程Simpson法四阶龙格—库塔法matlab一．问题重述水星到太阳的最远距离为110.698210m，此时水星绕太阳运行的线速度为43.88610m／s。试求问题一水星到太阳的最近距离问题二水星绕太阳运行的周期问题三从远日点开始的第50天（地球天）结束时水星的位置并画出轨道曲线二．问题分析求水星到太阳的最近距离以及水星绕太阳运行的周期等，需要先将水星轨道方程求出，因此可以根据Newton第二定律及万有引力定律222imMGdZemrdt，建立微PAGE1分方程模型，将原问题转化为求解带有初值条件的微分方程问题，进而采用解析法或数值方法求解远日点和周期。三．模型假设1．水星运行的轨道是以太阳为一个焦点的椭圆2．从太阳指向水星的线段在单位时间内扫过的面积相等3．水星运行周期的平方与其运行轨道椭圆长轴的立方之比为常量四．符号系统1．0v水星在远日点的线速度2.M太阳的质量3.m水星的质量4.or水星在远日点的距离5.T周期五．建立模型与求解模型一水星的轨迹方程设太阳中心所在的位置为复平面的原点O，在时刻t，水星位于()iZtre所表示的点P。这里(),()rrtt均为t的函数，分别表示()Zt的模和辐角。于是水星的速度为()iiidZdrddrdeireeirdtdtdtdtdt，加速度为2222222(())(2)idZdrdddrderirdtdtdtdtdtdt（1.1），而太阳对行星的引力依万有引力定律，大小为2mMGr，方向由行星位置P指向太阳的中心O,故为2imMGer，其中301.98910()Mkg为太阳的质量，m为水星的质量，11226.67210(/)GNmkg为万有引力常数。PAGE1依Newton定律，我们得到222imMGdZemrdt(1.2)，将（1.1）代入(1.2)，然后比较实部与虚部，就有22222220()ddrdrdtdtdtdrdMGrdtdtr这是两个未知函数的二阶微分方程组。在确定某一行星轨道时，需要加上定解条件。假设当t=0时，行星正处于远日点，而远日点位于正实轴上，距原点O为0r，行星的速度为0v。那么就有初值条件：000000000ttttrrdrdtvddtr因此问题转化为求解带初值问题的微分方程组222222000000020()00ttttddrdrdtdtdtdrdMGrdtdtrrrdrdtvddtr又将2220ddrdrdtdtdt两边同乘以r，即得2()0ddrdtdt，从而21drcdt（1.3）,其中100crv，这样有向线段OP在时间t内扫过的面积等于21122tttctdrdtdt，这个正是Kepler的第二定律，从太阳指向水星的线段在单位时间内扫过的面积相等。PAGE1将（1.3）代入2222()drdMGrdtdtr得221232cdrMGdtrr，于是我们可以得到水星运行的较为简单形式的数学模型：22123212000000tttcdrMGdtrrcddtrrrdrdt为了求得行星的轨迹方程，要消去变量t，令1ru，那么12cddtr可以改写为21dcudt从而222211122()drddududducccudtdtdddtd将上式代入221232cdrMGdtrr，化简后为221duudp(1.4)，其中21cpMG，引进1uup，立即可以求出01cos()uuAp，这里A和0是待定的常数。记eAp，上式可以写为01cos()pre这个就是水星的轨道方程，是一条平面二次曲线。由于水星绕太阳运行，故必有01e。由于r在t=0时取道最大值0r（远日点），这个就意味着此时函数0cos()取道最大值1.于是就有000,1per，从而轨迹方程为1cospre。对于水星而言，114000.69810(),3.88610(/)rmvms，又水星的近日点到太阳的距离1cos1mppree。依据已知数据，可知PAGE11521002.71310(/)crvms，21015.54710()cpmMG，010.2055per，从而计算水星到太阳的最近距离为)(104.601610mrm模型二水星的运行周期设水星的周期为T，那么利用Kepler第二定律，我们有2101122TdrdtCTdt（1.4）上式左端为水星轨迹椭圆所围的面积，记为S，由于椭圆的半长轴21epa，半短轴21epb，从而有2322)1(epabS将上式代入式（1.4）,解得23212)1(2eCpT（1.5）将有关数据代入，易得87.9919(d))(107.60256sT模型三水星的位置由于水星的运行满足Kepler第二定律，则该式可改写为tCdr12，从而可得0212)cos1(tdeCp如果我们要求1Tt时相应的和r，则意味着首先要解方程，211)(pTCF，其中deF02)cos1(1)(在求出了1Tt时的后，立即可以由1cospre得到相应的r。下面用数值方法求解水星的位置1.Simpson法PAGE1由被积函数21(1cos)e的恒正性可知()F单调，从而方程211)(pTCF的根必存在且唯一。取,(1,2,...)khkhk，记()kkFF。若1111122,nnCTCTFFpp，那么位于n与1n之间，在h适当小时，可取n。计算()F可采用不同的数值积分法，本文采用Simpson法，取步长h=0.001，具体求解过程见附录一，最后结果为3.791，104.76710r2.基于压缩映像的求根方法我们引入水星轨道椭圆的参数方程，由于椭圆的半长轴21epa，半短轴21epb，从而中心到焦点的距离为aeba22。因左焦点为原点，故椭圆中心位于（ae，0），于是得到参数方程(cos)sinxaeyb它们与,r的关系为tan,222xyryx此式可改写成1('')[sin()sin]Ctxyyxdabeeab当1Tt时解方程abTCe11sin记abTC11，sin)(eg，那么上式即()g，就是说要去求函数()g的不动点，求解方程不动点可以采用简单迭代法，对于水星，我们已计算出0.2055e，由于e很小，因此迭代收敛理论上可以很快，当时间从远日点开始的第50天结束时，意味着)(10432.071sT,从而3.5703)1(2322111epTCabTCPAGE1不妨取00，于是3.5703sin01e3.6557sin12e......54sin3.6747e65sin3.6747e故6747.3由式sin),cos(byeax，tan,222xyryx，可以计算出相应的，即由0.75849)cos(sintaneab得0.64891，而3.791此时的距离r为1022104.7668]sin[)]cos([bear（m）3.经典四阶Runge-Kutte法由我们将由最初的微分方程组求解水星的位置，方程组见下22123212000000tttcdrMGdtrrcddtrrrdrdt令dtdrq，那么我们可以得到一阶微分方程组：PAGE1000000212321tttdtdrrrrCdtdqdtdrrMGrCdtdq若记这个微分方程组中方程的右端依次为),,,(),,,(),,,,(rqtSrqtRrqtQ和，则相应的四阶Runge-Kutte迭代格式法为)22(643211KKKKhqqkk)22(643211LLLLhrrkk)22(643211NNNNhkk这里对于11234(22)6kkhqqKKKK，有),,,(1kkkkrqtQK)2,2,2,2(1112hNhLrhKqhtQKkkkk)2,2,2,2(2223hNhLrhKqhtQKkkkk),,,(3334hNhLrhKqhtQKkkkk初值为0,,00000rrq，则对于给定的步长值h，类似可以逐步计算一系列的,,kkkqr，由于行星绕着太阳运行，只需取12,2nn，而取得行星轨道上一系列点的近似坐标（,kkr），再通过极坐标与直角坐标的转换，继而可以绘出轨道曲线。通过matlab编程求解得3.802，104.77910r，轨道曲线如下PAGE1程序见附录二。六．模型推广本文建立的微分方程模型对于求解行星绕日运行轨道具有广泛的应用空间，只需给出行星的远日点和在远日点的运行线速度即可计算出轨道方程，用数学软件绘出近似的轨道曲线，对于研究天体运行有所帮助。此外，本文采用的求解微分方程的数值方法，具有较为快速且准确的收敛效果，可以用来求解其他类似的微分方程模型。七．参考文献【1】乐经良，数学实验，北京，高等教育出版社，1999年10月【2】周品，matlab数值分析，北京，机械工业出版社，2009年1月PAGE1八．附录附录一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;PA