行星的轨道和位置

行星的轨道和位置高路（船舶海洋与建筑工程学院5120109107）一、 背景介绍16世纪以前，人们都认为行星绕太阳旋转的轨迹是圆。17世纪初，在丹麦天文学家T.Brache观察工作的基础上，Kepler提出了震惊当时科学界的行星运动三大定律：1.行星运行的轨道是以太阳为一个焦点的椭圆；2•从太阳指向某一行星的线段在单位时间内扫过的面积相等；3•行星运动周期的平方与其轨道椭圆长轴的立方之比值是不随行星而改变的常数。对这三条定律的分析和研究导致Newton发现了著名的万有引力定律，而同时，应用万有引力定律，Kepler的行星运动三大定律得到了理论上的推导。由于行星间引力的存在，基于万有引力定律的计算表明：行星的轨道应该是稍偏于以太阳为焦点的椭圆。计算结果与天文学家测得的实际结果在木星、土星等行星的轨道上相当吻合，然而在天王星的轨道上却存在着不容忽视的误差。当时人们只发现了太阳系的七大行星，天王星是其中最后发现的（1781年），于是科学家们猜想：还存在影响天王星运行轨道的其他行星。1864年，Adams（英）与Leverrier（法）分别推算出这颗可能存在的行星的位置，同年，天文学家就在他们推算的方位上找到了海王星。由于这颗行星的发现首先依赖于根据万有引力定律的计算，因此它被称为“铅笔尖上的行星”此后，仍是类似的猜想和推算导致了质量较小的冥王星被发现，这充分说明了Newton万有引力定律这样一个数学模型的正确性和重要性。二、 实际问题水星距太阳最远处（远日点）距离为6.982X1010m，此时地球绕太阳运动（公转）的速度为3.886X104m/s，试求：（1） 地球距太阳的最近距离；（2） 地球绕太阳运转的周期；（3） 在从远日点开始的第50天结束时，地球的位置与速度。三、 数学模型设太阳中心所在位置为复平面之原点0，在时刻t，行星位于

所表示的点p。这里丫=丫C）, ）均为才的函数，分别表示ZC）的模和辐角。于是行星的速度为dZdr. d。= ei。+zrei。dtdt dt=ei°dr.d。'=ei°一+ir dtdt丿其加速度为d2Z=d2Z=ei。dt2Idt丿+/r^!+2dLd^Vdt2 dtdt(2)mMG而太阳对行星的引力依万有引力定律,大小为 ，方向由行星位置p指•m2/kg2）为万有引力常数。向太阳的中心0，故为一mMGei。，其中M=1.989•m2/kg2）为万有引力常数。量，m是行星的质量，G二6.672x10-ii依Newton 第二定律，我们得到(3)mMG d2Z(3)-ei。=m—

r2 dt2将式（5.2）代入式将式（5.2）代入式（5.3），然后比较实部与虚部，就有d20 drd。八dt2 dtdt（d0、2d2r-rdt2MG(4,5)Vdt丿 r2这是两个未知函数的二阶微分方程组。在确定某一行星轨道时，需要加上定解条件。我们设当t=0时，行星正处于远日点，而远日点位于正实轴上，距原点0为［，行星的线速度为v0，那么就有初始条件:r=rt=0 0=0(6〜9)lt=(6〜9)v=_0.tv=_0.t=0 r0dt问题（4）〜（9）就是行星绕太阳运行的轨迹的数学模型。将式（4）乘以r,即得d( do)r2dtvdt丿从而dordor2_dt=C仙数)1(10)其中C=rv100ft+At1ft+At1r2d0Jt 2~di(11)”CAt

dt=1-(11)2显然，这正是Kepler的第二定律：从太阳指向行星的线段在单位时间内扫过的面积相等。将式（10）改写后代入式（5）得d2rC2MG —1=— dt2 r3 r2于是我们得到了行星运动的形式较为简单的数学模型：

d2rC2MGTOC\o"1-5"\h\z —1=— dt2 r3 r2d0 C=idtr2vr| =r00(12〜(12〜16)dt四、解析方法行星的轨迹1为求得行星的轨迹方程，要消去变量t，令r=—,那么式（13）可以写为udodtdodt=CU2 1(17)从而dr= 从而dr= 1du= 1dud0dtu2dt u2dtdu—C1dodt2=—C1ddt2=—C1d(du、

dtvd0丿=cd2ud01d02dt=—C2u21d2u将上式代入式（12）,简化后为d2u 1(18)+u=一(18)do2 pC2

其中p=Mg•-1式（18）是一个二阶常系数非齐次微分方程，引进u=u—，立即可以求P

u—一=u=acoSe—0)

p 0这里A和00是待定的常数。记e=AP，上式可写为(19)r= P(19)1+ecoS0—0）0这就是行星的轨道方程，是一条平面二次曲线。由于行星绕太阳运行，故必有0<e<1，这样我们得到了Keple第一定律：行星的轨道是以太阳为一个焦点的椭圆。由于r在t=0时取到最大值rQ（远日点），而0=0，这意味着此时函数COS（-00）取到最小值-1，于是就有e二兀,0从而行星轨迹的方程为r现在回到开始的问题，对水星而言，r=6.982x现在回到开始的问题，对水星而言，r=6.982x1010（m）,0v=3.886x104

0显然，水星的近日点到太阳的距离r由在式（19）中取0=兀得到，即mr=丄m1+e依已知数据，可知沁2.713沁2.713x1015p=MG〜5.54771010s)沁0.205498259）9040r0从而计算得地球到太阳的最近距离为ru4.602x1010（m）m（在计算过程中多保留一位有效数字）。II、行星的周期

设行星的周期为t，那么利用Kepler第二定律即式（11）,我们有JT1r2dLdt=1CT02dt21上式左端为行星轨迹椭圆所围的面积，记为S,由于椭圆的半长轴a=幺，半短轴b任，从而有S—兀ab—将上式代入式（20）,解得T=_C1a=幺，半短轴b任，从而有S—兀ab—将上式代入式（20）,解得T=_C1对于水星，将有关的数据代入，易有T沁7.602x106d87.99）(21)III、行星的位置由于行星的运行满足KePler第二定律即式（5.11）,而该式可改写为f0+A0r2do二CAt9 i(22)从而可得J0z Pd0=t 0C（1-eco0）21如果我们要求t=t时相应的0和r，则意味着首先要解方程(23)其中F(0)=CZP2(24)F(0)=f0z 1 七d00\1-eco0)2在求出了t=T时的0=0后，立即可以由式（19）得到相应的丫，再用式（13）求出此时行星的角速度继而得到线速度。五、数值方法1.一般的方法由于被积函数的恒正性知F（0）单调，从而方程（24）的根仔必唯一存在.取Ao二h,0二kh(k二Ao二h,0二kh(k二1,2,…),记F<k亠L,F1〉一」,那么0位于0np2 n+1 p2 n与0L之间，在h适当小时，可取0-0.下面我们分别用矩形法、梯形法、n+1 n辛普森法计算F（0）.⑴矩形法左矩形法相应的MATLAB 程序为：functionm1_11(h)ep=0.20550;C1=2.713e15;p=5.547e10;Tl=50*24*3600;f=C1*Tl/p入2;theta(1)=0;F(1)=0;fori=2:1e6F(i)=F(i-1)+h*(1-ep*cos(theta(i-1)))入-2;theta(i)=theta(i-1)+h;ifF(i)>fbreak;endendn=i-2t=n*hr=p/(1-ep*cos(t))dtheta=C1/r入2v=r*dtheta分别取不同步长得到的结果:hn0"r0'v0.05753.75004.74661.20425.71570.013783.78004.76131.19685.69810.0057573.78504.76381.19555.69510.00137903.79004.76631.19425.6920和104m/s.•右矩形法相应的MATLAB程序为：functionm1_12(h)ep=0.20550;C1=2.713e15;p=5.547e10;Tl=50*24*3600;f=C1*Tl/p入2;theta(1)=0;F(1)=0;fori=2:1e6theta(i)=theta(i-1)+h;F(i)=F(i-1)+h*(1-ep*cos(theta(i)))入-2;ifF(i)>fbreak;endendn=i-2t=n*hr=p/(1-ep*cos(t))dtheta=C1/r入2v=r*dtheta分别取不同步长得到的结果:hn0"ro'v0.05763.80004.77141.19175.68590.013793.79004.76631.19425.69200.0057583.79004.76631.19425.69200.00137913.79104.76681.19405.6914其中，h，n,0，r，0和v的单位分别是107s，次，rad(弧度)，lOiom,10-6rad/s和104m/s.•中矩形法相应的MATLAB程序为：functionm1_13(h)ep=0.20550;C1=2.713e15;p=5.547e10;Tl=50*24*3600;f=C1*Tl/p入2;theta(1)=0;F(1)=0;fori=2:1e6theta(i)=theta(i-1)+h;F(i)=F(i-1)+h*(1-ep*cos((theta(i)+theta(i-1))/2))入-2;ifF(i)>fbreak;endendn=i-2t=n*hr=p/(1-ep*cos(t))dtheta=C1/r入2v=r*dtheta分别取不同步长得到的结果:hn0"r0'v0.05753.75004.74661.20425.71570.013793.79004.76631.19425.69200.0057583.79004.76631.19425.69200.00137903.79004.76631.19425.6920其中，h,n，厂，厂,0和v的单位分别是107S,次，rad(弧度),10iom,10-6rad/s和104m/s.⑵梯形法相应的MATLAB 程序为：functionm1_2(h)ep=0.20550;C1=2.713e15;p=5.547e10;Tl=50*24*3600;f=C1*Tl/p入2;theta(1)=0;F(1)=0;fori=2:1e6theta(i)=theta(i-1)+h;F(i)=F(i-1)+h*((1-ep*cos(theta(i-1)))入-2+(1-ep*cos(theta(i)))入-2)/2;ifF(i)>fbreak;endendn=i-2t=n*hr=p/(1-ep*cos(t))dtheta=C1/r入2v=r*dtheta分别取不同步长得到的结果：hn0r0v0.05753.75004.74661.20425.71570.013793.79004.76631.19425.69200.0057583.79004.76631.19425.69200.00137903.79004.76631.19425.6920其中，h,n，厂,r,0和v的单位分别是107s,次，rad(弧度),10iom,10-6rad/s和104m/s.⑶辛普森法相应的MATLAB 程序为：functionm1_3(h)ep=0.20550;C1=2.713e15;p=5.547e10;Tl=50*24*3600;f=C1*Tl/p入2;theta(1)=0;F(1)=0;fori=1:1e6theta(2*i)=theta(2*i-1)+h;theta(2*i+1)=theta(2*i)+h;F(2*i+1)=F(2*i—1)+h*((1—ep*cos(theta(2*i—1)))入—2+4*(1—ep*cos(theta(2*i)))入-2+(1-ep*cos(theta(2*i+1)))入-2)/3;ifF(i)>fbreak;endendn=i-1t=n*hr=p/(1-ep*cos(t))dtheta=C1/r入2v=r*dtheta分别取不同步长得到的结果:hn0r0v0.05763.80004.77141.19175.68590.013803.80004.77141.19175.68590.0057603.80004.77141.19175.68590.00137923.79204.76731.19375.69080.0001379083.79084.76671.19405.6915

其中，h,n,0,r,0和v的单位分别是107s,次，rad(弧度),10iom,10-6rad/s和104m/s.II.基于压缩映射的求根方法首先，回到行星的轨道曲线，我们将引进轨道椭圆的参数方程。由于椭圆的半p 一］p长轴a= ，半短轴b= ，从而中心到焦点的距离为1-e2 1-e2W2-b2=ae。因左焦点为原点，故椭圆中心位于(ae,0),于是得到参数方程(25)它们与r，0的关系为x=a(e+cos屮)，y(25)它们与r，0的关系为y(26)x2+y2=r2, =tan0(26)空:竺d空:竺d^X2+y2 ，从而式(22)又可以改写成d0这样由于CAt=卜+旳(xy'-yx')d^申=J®+A申la(e+cos^)bcosQ-bsinQ(-asinQ)l/Q=fQ+AQab(ecosQ+1)dQ=ablesinG+Aq)-esinq]+abAQ于是，我们要求t=t时的行星位置就意味着要解方程esin叶一CTG)二九一esin申，那么上式即(27)这就是说要求函数gC)的不动点，由于|g'(p)=-ecos<e<1

从而gG）有性质:(28)gG）-gG）(28)121由于e<1,满足式（28）的函数（或映射）g被称为压缩的。不难得到gG）有唯一的不动点孑，即式（27）有唯一的根犷。而且如果令(29)¥=gG）k=1,2,(29)k k-i则任取申0，将有由于则任取申0，将有由于¥T®(kk二g©)-gGJk k-11因此，当e很小时，¥7收敛的速度是很快的。k当行星是水星时，我们已经算出e沁0.20549869250904，当时间为从远日点开始的第50天结束时，意味着T=0.432X107（s），从而CTCT（ \九=11= 111—e2召abp2沁3.57034不妨取¥0=0，于是依照迭代格式（29），就有Q=九一esinQ=3.57034TOC\o"1-5"\h\z1 0Q=九―esinQ=3.57034-0.20550sin3.57034=3.655771Q=九―esinQ=3.57034—0.20550sin3.65577=3.671412Q=九—esinQ=3.57034—0.20550sin3.67141=3.674193Q=九—esinQ=3.57034—0.20550sin3.67419=3.674694Q=九—esinQ=3.57034—0.20550sin3.67469=3.674775Q=九—esinQ=3.57034—0.20550sin3.67477=3.674796Q=九—esinQ=3.57034—0.20550sin3.67479=3.674797这表明gC）的不动点即式（27）唯一的根为/-3.67479。注意这迭代格式用了八次就收敛到不动点，因此这是一种快速的收敛格式。据式（25），（26）可以将犷转化到相应的厂，即由tan伎=?sinQ_）=0.75866a\+cosQ丿计算得 arctan（0.75866）匕0.64902而实际上e＞兀所以，得「=3.79061此时的距离r和线速度v分别为r=、Ja（+cosqX+bsinQ】=4.766806342759505x1010mCv=r0=—L=5.6919x104m/sr与前面的数据加以比较，容易看出，这里基于不动点的快速收敛迭代格式具有不容置疑的优势，它不但快速，而且精确。Runge—Kutta法Runge一尺加加法中用得最多的格式是经典的（四阶）Runge一Kutta法。其迭代计算格式为

(K+2K+2K+K)TOC\o"1-5"\h\z6 1 2 3 4这里=fG,y)kk\o"CurrentDocument"( h 1 \K=fT+-,y+—hK( k 2 k 2 2丿\o"CurrentDocument"( h 1 \K=fT+-,y+—hK(31)( k 2 k 2 2(31)K=fG +h,y+hK)k k 3Runge-Kutta法是-种单步法，所谓单步法是指在迭代计算过程中第k+1步的值仅仅依赖第k步的值，因此单步法在使用时较为方便、简单。而且可以在计算过程中改变步长。在软件MATLAB 中有现成的用Runge-Kutta方法求解常微分方程（或方程组）的专用程序，只要调用这个程序，使用是十分方便的。现在我们将Runge―知加法应用到水星绕太阳运行的数学模型（12）〜（16）,dr令dr令q=冠，那么我们得到一阶微分方程组:TOC\o"1-5"\h\zdq C2 MG__=_1_— \o"CurrentDocument"dt r3 r2dr丿=qdtd0 C =_Ldtr2 (32)r=rt=0 0q=0t=0e=0Jt=0若记这个微分方程组中方程的右端项依次为QC，q，r，0）,R(t,q,r,0)和S(t,q,r,0),则相应的四阶Runge-R(t,q,代格式为q=q+h(K+2K+2K+K)TOC\o"1-5"\h\zk+i k6 1 2 3 4r=r+h(L+2L+2L+L)k+1 k61 2 3 40 =0+h(N+2N+2N+N)k+1 k6 1 2 3 4其中K，L和nC=^2,3,4)类似于单个方程时的形式(31)。由初始值i i iq=0,r=r,0=00 000给定步长h的值，就可以逐步算得一系列的qr，0，由于行星绕着太阳运kkk行，我们只要计算到0<2兀,0 >2兀n+1而取得行星轨道上一系列点的近似坐标c,%），继而可通过绘图来得到轨道kk的曲线。同时可以取相应于0的时间作为周期，即取n行星到太阳的最近距离r则可取行星到太阳的最近距离r则可取t=mn为2的整数部分）相应的近似坐标中的r。m使用MATLAB ，先建立常微分方程组函数m2_1_fun.m:functiondy=m2_1_fun(t,y)C1=2.713e15;MG=1.989e30*6.672e-11;dy=zeros(3,1);dy(1)=C"2/y(2)入3_MG/y(2)入2;dy(2)=y(1);dy(3)=C1/y(2)入2;然后，通过如下方式调用Runge―Kutta法的专有程序，建立m文件m2_1.m:functionT=m2_1(h)[t,y]=ode45(@m2_1_fun,[0:h:100*24*3600],[0,6.982e10,0]);n=max(find(y(:,3)<2*pi));T=t(n)rm=y(round(n/2),2)polar(y(:,3),y(:,2))theta=y(round(50*24*3600/h),3)r=y(round(50*24*3600/h),2)取不同的步长就可以求得水星轨道运行的周期T、到太阳的最近距离r和从远m日点开始第50天结束时水星的位置，同时得到水星的公转轨道。hTTrmr00.01750000086.814.60594.69693.64600.005755000087.384.60104.72243.70710.0005759500087.914.59924.75923.78560.0001759700087.934.59924.76163.7904其中，h，T，丫，r和0的单位分别是s,d，1010m,1010m和rad（弧度）.m下面是步长h=0.0001X107s时，由MATLAB 所作的在极坐标下的水星公转轨迹.9080000000000270六、实验总结这一次我选做的数学实验是“行星的轨道与位置”这个实验是以Kepler三大定律和Newton第二定律为理论依托，实验任务主要是求解行星绕太阳的公转轨迹、近日点、远日点和在给定时刻的位置并完成由KePler三定律到万有引力定律的推导。我主要完成了教材所给的任务2、3,其中任务3以附录的形式给出，而任务1由于与任务2中的完全类似，故将其略去不做。经过上一次的实验及后续几节课老师的讲解，我对MATLAB有了更深认识，运用起来也相对更熟练一些。这次我自己编写的几个MATLAB程序都得到了令我满意的结果，较上次有了很大的进步。总之，我希望在接下来的一个月的时间里，能够继续认真学习MATLAB的相关应用，争取更好地完成最后一个数学实验。七、参考文献⑴乐经良，向万隆，李世栋，等.数学实验.第二版.北京:高等教育出版社，2011.⑵张志涌，杨祖樱，等.MATLAB教程R2010a.北京:北京航空航天大学出版社,2010.附：从开普勒定律到牛顿万有引力定律第一部分准备―、极坐标中的椭圆方程椭圆定义为到定点的距离与到定直线的距离之比为常数e的点的集合。如图1所示，在极坐标中，Ox为极轴l是垂直于极轴的定直线，它与O点r=ep+rcos0整理可得:(1)pe(1)r= 1-ecos0二、极坐标中的位置矢量极坐标中，r表示从原点到曲线上一点的距离，如果我们以原点0为参考，则r实际上只表示出了位置矢量的大小。为了明确其方向，我们沿着r所在的直(2)线做出单位矢量i作为径向单位向量。另外，将i旋转|得到j作为横向单位向量。显然物体的位置矢量可表示为：(2)r=ri上式中等号右边的r表示的是位矢的大小，i表示的位矢的方向。但是应当注意的是，不管是r还是i，都不一定是常量。这和直角坐标系中的单位向量是常量是有区别的。

图2图2另外，r和i都是e的函数，在运动学中e又是时间t的函数。所以，r和i都是时间t的函数，所以我们也可以说位置矢量r是时间的函数。在这里，我们必须清楚的是，极坐标中的矢量表示和用极坐标表示函数关系并不完全是一回事。若用极坐标表示数量关系，我们只需要用标量式r=rG）即可，在表示矢量时，我们不得不在这个基础上加上了单位向量i。三、极坐标中的速度和加速度下面我们先求单位向量对时间的导数。图3图3在图3中，以Ox方向为x轴，O为原点，垂直Ox向上为y轴建立直角坐标系，用g、耳表示沿x轴、y轴的单位向量，则i、j可分别表示为：i=cos随+sin的j=cos-sin +j=cos因此di=d(cos+sin旳)=d(cos+sin的)dB=($込gg+cos的)dt dt dB dt dt对比j的表达式有,d=gjdt其中o'表示e对时间的导数：同理可知：dtF面我们对位矢函数r=ri求导，这样可以得到在极坐标系中的速度公式。(5)dr di -(5)v二二ri+r=ri+rOj

dt dtdv dv dridrOja= = +—dt dtdt.di'ri+r—+ri+rOj+rOj+rOdjdt丿jOj+rOj—rO2i丿yO+2rO)jri+其中a=r—rO2 (6)ra二rO+2rO (7)O分别表示径向加速度和横向加速度。第二部分推导开普勒定律的内容是：开普勒第一定律，也称椭圆定律:每一个行星都沿各自的椭圆轨道环绕太阳,而太阳则处在椭圆的一个焦点上。开普勒第二定律，也称面积定律：在相等时间内，太阳和运动着的行星的连线所扫过的面积都是相等的。开普勒第三定律，也称调和定律：各个行星绕太阳公转周期的平方和它们的椭圆轨道的半长轴的立方成正比。由（7）式可知：a=rO'+2rO=1CrrO+r2O' —""刃r rdt由开普勒第二律可知:r29'=常数dr29丿小故上式中—— =0dt这就是说,r9r9+2r9=0 由椭圆方程r=pen可得:1-e由椭圆方程r=pen可得:1