行星的轨道和位置的数学解法

行星的轨道和位置的数学解法

作者：石磊a ，林川b 指导教师：乐经良C 教授

a : 上海交通大学电子信息与电气工程学院F0303032班(5030309885) , 电话：54740807

b : 上海交通大学电子信息与电气工程学院F0303032班(5030309880) , 电话：54741769

c : 上海交通大学理学院数学系

摘要：本文主要涉及常微分方程及对微分方程的建模与求解，数值积分的计算；利用多种微分方程的数值方法求解得到行星运行的参数和位置。研究基于压缩映象的求根方法和微分方程的Runge-Kutte 法。特别对Runge-Kutte 法进行较深入的讨论。并通过数值方法解微分方程得到的行星位置演示水星和冥王星的运行轨道，编制软件。

关键词：微分方程 数值方法 Runge – Kutte 法

问题的重述

17世纪初，在丹麦天文学家T.Brache 观察工作的基础上，Kepler 提出了震惊当时科学界的行星运行三大定律：

1. 行星运行的轨道是以太阳为一个焦点的椭圆；

2. 从太阳指向某一行星的线段在单位时间内扫过的面积相等；

3. 行星运行周期的平方与其运行轨道椭圆长轴的立方之比值是不随行星而改变的常数。

对这三条定律的分析和研究导致Newton 发现了著名的万有引力定律，而同时，应用万有引力定律，Kepler 的行星运行三大定律得到了理论上的推导。

数学模型

设太阳中心所在位置为复平面之原点O ，在时刻t ，行星位于

θi re t Z =)( （4.1）

所表示的点P 。这里)(),(t t r r θθ==均是t 的函数，分别表示)(t Z 的模和辐角。 于是行星的速度为

⎪⎭

⎫ ⎝⎛+=+=dt d ir dt dr

e dt d ire e dt dr dt dZ i i i θθθθθ 其加速度为

⎥⎥⎦⎤⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛++⎪⎪⎭⎫ ⎝⎛⎪⎭⎫ ⎝⎛-=dt d dt dr dt r d r i dt d r dt r d e dt Z d i θθθ22222222 (4.2） 而太阳对行星的引力依万有引力定律，大小为2

r

mMG

，方向由行星位置P 指向太阳的中心O ，故为θi e r

mMG 2

-，其中)(10989.130

kg M ⨯=为太阳的质量,m 是行星的质量，)/(10672.62211kg m N G ⋅⨯=-为万有引力常数。

依Newton 第二定律，得到

222dt

Z

d m

e r mMG i =-θ （4.3）

将（4.2）代入（4.3），然后比较实部和虚部，有

⎪⎪⎩⎪

⎪⎨⎧-=⎪⎭⎫ ⎝⎛-=+22

222

202r MG dt d r dt

r d dt d dt dr dt d r θθθ (4.4)(4.5) 如设0=t 时，行星正处于远日点，远日点位于正实轴上，距原点O 为0r ,行星的线速度为0v ,那么就有初值条件：

⎪⎪⎪⎪⎩⎪⎪⎪⎪

⎨⎧========000000

000r v dt

d dt dr r r t t t t θθ (4.6）(4.7）(4.8）(4.9）

方程(4.4) ~ (4.9)就是行星绕太阳运行的轨迹的数学模型。

将式（4.4）乘以r ，即得

02=⎪⎭

⎫ ⎝⎛dt d r dt d θ 常数）( 12

C dt

d r =⇒θ

（4.10） 其中 001v r C =

这样，有向线段在时间t ∆内扫过的面积等于

2

2112t C dt dt d r t

t t

∆=⎰

∆+θ

(4.11)

而这正是Kepler 第二定律。 将式(4.10)改写后代入式(4.5)

2

32

122r

MG

r dt r d C -=- 由此得到行星运动的较为简单形式的数学模型：

⎪⎪⎪⎪⎪⎪⎩

⎪

⎪⎪

⎪⎪⎪⎨⎧====-=-===00000021

2

321

22t t t dt dr

r r r C dt d r MG r dt

r d C θθ

实验任务

1． 在求解方程(4.24)时,试用矩形法、梯形法和Simpson 法来计算数值积分，

并对所得的结果加以比较？

解答：由于行星的运动满足Kepler 第二定律式(4.11),改写该式为

t C d r ∆=⎰∆+12θ

θθ

θ

()

t d e C p =-⇒⎰

θ

θθ0

2

12

cos 1 如果要求1T t =时相应的θ和r ，则意味着首先要解方程

()21

1p

T C F =

θ 其中

()()

⎰

-=θ

θθθ0

2

cos 11

d e F (4.24)

利用矩形法计算次积分，并带入水星数据，得

Clear r0,v0,h,c1,p,e,F,T1,k,1,2,r,,v; h0.0001;

r00.69821011;

v0 3.886104;

M 1.9891030;

G 6.6721011;

c1r0v0;

p c12M G;

e1p r0;

T150243600;

k0;

For F0,F c1T1p2,k,1k h;

F F h1e Cos12;

r p1e Cos1;

c1r2;

v r;

Print h,"",k1,"",1,"",r,"", ,"",v,"";

计算得到数据