关于水星的轨道方程以及运行周期等问题的探讨.doc

理学院07级计算2班*数学实验论文关于水星的轨道方程以及运行周期等问题的探讨摘要：运用万有引力定律和Newton第二定律，在复平面坐标系中建立了水星绕太阳运动轨迹的简单数学模型二阶微分方程组模型，然后通过解析法，求出了水星到太阳的最近距离为，再运用Kepler第二定律求出水星绕太阳运动的周期。又运用基于压缩映象的求根方法求出了以及从远日点开始的第50天（地球天）结束时水星的位置。接着又运用经典的四阶 Runge-Kutte法，通过MATLAB编程，计算出取不同的步长h时所对应的水星绕太阳运动的周期T以及从远日点开始的第50天（地球天）结束时水星的位置r，比较发现结果与前面方法所得的结果非常接近。最后利用MATLAB画出了水星绕太阳运动的轨迹曲线。关键字：水星 轨道 周期 基于压缩映象的求根方法 Runge-Kutte法一、模型建立：设太阳中心所在位置为复平面之原点O,在时刻t,行星位于 （4.1）所表示的点P.这里均为的函数,分别表示的模和辐角。 于是行星的速度为其加速度为 (4.2）而太阳对行星的引力依万有引力定律，大小为,方向由行星位置P指向太阳的中心O，故为,其中M太阳的质量，m为行星的质量，G为万有引力常数依Newton第二定律，我们得到 （4.3）将（4.2）带入（4.3），然后比较实部与虚部，就有 这是两个未知函数的二阶微分方程组。在确定某一行星轨道时，需要加上定解条件。我们设当t=0时，行星正处于远日点，而远日点位于正实轴上，距原点O为，行星的线速度为.那么就有如下的初值条件：问题（4.4）（4.9）就是行星绕太阳运行的轨迹的数学模型。将式（4.4）乘以r，即得从而 （4.10）其中为常数这样，有向线段在时间内扫过的面积等于 （4.11） 显然，这正是Kepler的第二定律：从太阳指向行星的线段在单位时间内扫过的面积相等。将式（4.10）改写后代入式（4.5）于是我们得到了行星运行的较为简单形式的数学模型：二、模型求解：1.行星的轨迹前面求得了行星的轨迹方程，要消去变量t，令，那么式（4.13）可以改写为 （4.17）从而将上式代入式（4.12），化简后为 （4.18）其中 式（4.18）是一个二阶常系数非齐次常微分方程，引进，立即可以求出这里A和是待定的常数。记上式可写为这就是行星的轨道方程，是一条平面二次曲线。由于行星绕太阳运行，故必有，这样我们得到了Kepler第一定律：行星的轨道是以太阳为一个焦点的椭圆上。由于r在时取得最大值（远日点），这意味着此时函数取得最大值1。于是就有从而行星轨迹的方程为 （4.19）对水星而言，我们已知由此可得显然，在式（4.19）中时，水星到达近日点，故其到太阳的距离为 带入上面的数据，从而计算得到水星到太阳的最近距离为2.运行周期设行星的周期为T，那么利用Kepler第二定律，即式（4.11），我们有 （4.20）上式左端为行星轨迹椭圆所围的面积，记为S我们已经求得椭圆的半长轴，半短轴，从而有将上式代入式(4.20),解得 (4.21)对于水星，将有关数据代入，易得3.行星的位置由于行星的运行满足Kepler第二定律，即式（4.11），而该式可改写为 （4.22）从而可得 （4.23）如果我们要求时相应的和，则意味着首先要解方程 （4.24）其中在求出了时的后，立即可以由式（4.19）得到相应的r.但是在实际问题中，即使可以积分出初等表达式，要由方程（4.24）来求出的解析解的可能性也是很小的。当然，我们可以利用计算编程利用迭代的方法，逐步增大的值来求得最符合条件的解，但这显然不是一个较为合理的方法。我们可以利用数值方法来对该问题进行求解，基于压缩映像的求根方法是其中较为有效的一种方法。其算法如下：我们已经求得椭圆的半长轴，半短轴，从而中心到焦点的距离为。因左焦点为原点，故椭圆中心位于，于是得到参数方程它们与的关系为故有当时解方程记，；那么上式即为如果令则任取，将有又由于所以当很小时，收敛的速度是很快的。当行星是水星时，我们已计算出，当时间从远日点开始的第50天结束时，意味着,从而不妨取，于是而由得 0.64891所以此时的距离r为4.下面我们再使用Runge-Kutte（四阶）来对刚才的几个问题进行求解令，那么我们可以得到一阶微分方程组：若记这个微分方程组中方程的右端依次为，则相应的四阶Runge-Kutte迭代格式法为其中 L,N的关系类似可得初值为 给定步长h的值，就可以逐步算得一系列的，又由于水星绕着太阳旋转，我们只要计算到而取得行星轨道上一系列点的近似坐标，继而可通过绘图来得到轨道的曲线。利用MATLAB编程计算可得：hTTrtheta0.058.00000092.59264.8893.9990.017.70000089.12044.8183.8850.0057.65000088.04174.7853.8260.0017.61000088.07874.7723.802其中前两列的单位为，第三列的单位为，第四列的单位为，最后一列的单位为。5.利用MATLAB画出水星绕太阳运动的轨迹曲线参考文献：1 乐经良 主编. 数学实验. 北京：高等教育出版社，20072 张学敏 主编. MATLAB基础及应用. 中国电力出版社, 20093 李庆扬，王能超，易大义编. 数值分析. 第五版. 北京：清华大学出版社，2009 附录：程序代码：MATLAB版本 R2008b（7.70） 公共函数function y1=Q(r)c1=2.7132e15;M=1.989e30;G=6.672e-11;y1=c12./(r.3)-M*G./(r.2);function y2=R(q)y2=q;function y3=S(r)c1=2.7132e15;y3=c1./(r.2);求周期Tclear;clc;format short h=0.05e7,0.01e7,0.005e7,0.001e7;%h为步长for i=1:4 q=0;r=0.6982e11;theta=0;T(i)=0; while theta=2*pi K1=Q(r); L1=R(q);N1=S(r); K2=Q(r+h(i)/2*L1);L2=R(q+h(i)/2*K1);N2=S(r+h(i)/2*L1); K3=Q(r+h(i)/2*L2);L3=R(q+h(i)/2*K2);N3=S(r+h(i)/2*L2); K4=Q(r+h(i)*L3);L4=R(q+h(i)*K3);N4=S(r+h(i)*L3); T(i)=T(i)+h(i); q=q+h(i)/6*(K1+2*K2+2*K3+K4); r=r+h(i)/6*(L1+2*L2+2*L3+L4); theta=theta+h(i)/6*(N1+2*N2+2*N3+N4); end;endday=T./(3600*24);disp(步长h对应取值为h=0.05e7,0.01e7,0.005e7,0.001e7时，其对应的周期T(单位s)为)disp(T)disp(步长h对应取值为h=0.05e7,0.01e7,0.005e7,0.001e7时，其对应的周期T(单位天)为)disp(day)求位置clear;clc;format short h=0.05e7,0.01e7,0.005e7,0.001e7;%h为步长for i=1:4 q=0;r(i)=0.6982e11;theta(i)=0;T=0; while T=50*24*3600 K1=Q(r(i); L1=R(q);N1=S(r(i); K2=Q(r(i)+h(i)/2*L1);L2=R(q+h(i)/2*K1);N2=S(r(i)+h(i)/2*L1); K3=Q(r(i)+h(i)/2*L2);L3=R(q+h(i)/2*K2);N3=S(r(i)+h(i)/2*L2); K4=Q(r(i)+h(i)*L3);L4=R(q+h(i)*K3);N4=S(r(i)+h(i)*L3); T=T+h(i); q=q+h(i)/6*(K1+2*K2+2*K3+K4); r(i)=r(i)+h(i)/6*(L1+2*L2+2*L3+L4); theta(i)=theta(i)+h(i)/6*(N1+2*N2+2*N3+N4); end;endformat longdisp(步长h对应取值为h=0.05e7,0.01e7,0.005e7,0.001e7时,)disp(从远日点开始的第50天结束时水星的位置（r,theta）为)disp( )disp(r,theta)画水星绕太阳运行的轨道曲线clear;clcformat short p=5.5472e10;e=0.2055h=0.01e7;theta=linspace(0,2*pi,100000);r=p./(1-e*cos(theta);x=r.*cos(theta);y=r.*sin(theta);plot(x,y,r)text(0,0,fontsize10太阳)text(0.6982e11,0,f