第十次作业

计算物理 郭忠智2014301020087

---

1 摘要

本次作业完成第4章的4.9题，观察行星运动轨道随着离心率e的变化而呈现的不同形状，研究在什么时间以后轨道发生明显的变化，文中展示了从运动开始到往后的一段时间内的运动轨迹的形状，并用Vpython画图展示了行星的位置随和时间的关系，得出运动轨道的离心率越大，轨道形状越偏离完美的圆形的结论。

2 背景介绍

In this section we saw that orbits are unstable for many value of $\beta$ that is not precisely 2 in (4.12).A related question,which we did not address is how unstable an orbit might be.That is,how long will it take for an unstable orbit to become obvious.Investigate this by studying orbits with the same value of $\beta$ and comparing the behavior with different values of ellipyicity of orbit.

3 正文

根据牛顿运动定律，太阳和地球之间的万有引力有如下表达式：

$$F_G=\cfrac{GM_SM_E}{r^2},$$
当我们忽略宇宙中其他物体的作用时，太阳和地球组成简单的两体系统，因为太阳的质量远比地球大，我们将地球绕着太阳的运动用以下微分方程来表达： $$\cfrac{d^2x}{dt^2}=\cfrac{F_{G,x}}{M_E}\\\cfrac{d^2y}{dt^2}=\cfrac{F_{G,y}}{M_E},$$
将$F_G$作矢量分解： $$F_{G,x}=-\cfrac{GM_SM_E}{r^2}\cos {\theta}=-\cfrac{GM_SM_Ex}{r^3}\\F_{G,y}=-\cfrac{GM_SM_E}{r^2}\sin {\theta}=-\cfrac{GM_SM_Ey}{r^3},$$
由此得到四个一阶微分方程： $$\cfrac{dv_x}{dt}=-\cfrac{GM_Sx}{r^3}\\\cfrac{dx}{dt}=v_x\\\cfrac{dv_y}{dt}=-\cfrac{GM_Sy}{r^3}\\\cfrac{dy}{dt}=v_y.$$
为了运算和绘图的方便，我们选取天文单位： $$1AU\approx1.5\times10^{11}m\\1yr\approx3.2\times10^7s.$$
由圆周运动的条件 $$\cfrac{M_Ev^2}{r}=F_G=\cfrac{GM_SM_E}{r^2},$$
得到
$$GM_S=v^2r=4\pi^2 AU^3/yr^2.$$
所以在天文单位制下$GM_S$是一个常数。
下面我们用数值算法来将运动图像展示出来，我们用Eular-Crommer方法改写以上得到四个一阶微分方程得到
$$v_{x,i+1}=v_{x,i}-\cfrac{4\pi^2x_i}{r_i^3} \Delta t\\x_{i+1}=x_i+v_{x,i+1}\Delta t\\v_{y,i+1}=v_{y,i}-\cfrac{4\pi^2y_i}{r_i^3} \Delta t\\y_{i+1}=y_i+v_{y,i+1}\Delta t,$$
这时，我们只需要给定初始条件就能得到需要的解。
两体运动可以转化为其中一体固定，另外一体以折合质量$\mu=\cfrac{m_1m_2}{m_1+m_2}$,半径$\overrightarrow{r}=\overrightarrow{r}_2-\overrightarrow{r}_1$绕着这个物体运动,如果 $$m_1\gt\gt m_2$$ 则有 $$\mu=m_2$$ ,且$\overrightarrow{r}$是行星相对于一个净值的太阳的位置。
轨道的极坐标方程为： $$\cfrac{d^2}{d\theta^2}(\cfrac 1r)+\cfrac 1r=-\cfrac{\mu r^2}{L^2}F(r),$$
其中$F(r)=\cfrac{GM_SM_P}{r^2}$,这个方程的解为： $$\cfrac 1r=(\cfrac{\mu GM_SM_P}{L^2})[1-e\cos(\theta+\theta_0)].$$

这是一个椭圆方程，其半长轴为$a=L^2/[\mu GM_SM_P(1-e^2)]$,其中L为角动量$L=\mu r^2\dot{\theta}.$
在椭圆的近地点有

$$r_{min}=a(1-e)\\v_{max}=\sqrt{GM_S}\sqrt{\cfrac{1+e}{a(1-e)}(1+\cfrac{M_P}{M_S})},$$
远地点有 $$r_{max}=a(1+e)\\v_{min}=\sqrt{GM_S}\sqrt{\cfrac{1-e}{a(1+e)}(1+\cfrac{M_P}{M_S})}.$$
当我们给定初始条件满足上述的远地点和近地点条件时，将会得到椭圆轨道。
以上的万有引力遵循平方反比律： $$F_G=\cfrac{GM_SM_P}{r^{\beta}},\beta=2.$$
程序代码如下：
problem4.9代码

4 结论

结果如下：
下面分别给出不符合平方反比律，并且离心率逐渐增大时的各运动情况：
e=0.1

e=0.3

e=0.5

e=0.7

e=0.9

地球的离心率为e=0.017,取$\beta=2$即符合平方反比律，代入得到如下结果，可以看到其轨道很好地近似为圆形：

通过以上画图可以看到随着离心率的变化，轨道偏离完美的圆形。
我们也可以观察当违反平方反比律时轨道的形状：

下面是$\beta=2.5$时的情况，此时的轨道已经变得杂乱无章：

4 致谢

代码并非自己独立完成，其中参考了13级学长臧之昊的代码，以及vpython中自带的orbits程序。