Ex_10：Problem4.9 orbits are unstable for any value of β（β=2）

Author： 刘永杰 材料 2014301020094

背景资料

行星运动是天文学的重要研究方向。对于太阳系中行星的运动，我们既有解析解，又有数值解，是发展数值计算方法的良好实验室。本文使用Euler-Cromer方法，对开普勒定律以及行星运动进行了研究。同时回答了4.7题。首先考察一个行星绕太阳的运动的动力学方程。由开普勒定律可知，行星会做圆周运动，或更一般的椭圆运动。由万有引力定律和牛顿第二定律可以写出如下的动力学方程：

在这个方程中，我们选择了太阳为静止状态。这是由于行星的质量与太阳相比太小了。在后文中我们将会看到两者质量可比时的运动状态。通过查询资料可知太阳系八大行星的性质如下：

太阳系的部分3D模拟

在下面我使用Vpython模拟了太阳系中水星、金星、地球、火星的运动,,所用的数据均来自上表。之所以没有加入其它行星，是由于从上表可以看出，其它行星距离太阳的距离太远，如果加入在图中就看不到前几颗行星了。

平方反比定律失效的情况（4.9）

当平方反比定律失效时，行星的轨道将会发生进动。这里我们考察当beta等于2.1时，对于同一个beta，当离心率不同时，轨道的去向的变化过程。

代码

import matplotlib.pyplot as plt
import math
# this function gives trac
def trac(vyy):
    x=[1]
    y=[0]   
    vx=[0]
    vy=[vyy]
    t=[0]
    dt=0.002
    while t[-1]<=20:
        r2=(x[-1])**2+(y[-1])**2
        r=math.sqrt(r2)
        ox=-4*math.pi*math.pi*x[-1]/(r**3.1)
        oy=-4*math.pi*math.pi*y[-1]/(r**3.1)
        vx.append(vx[-1]+ox*dt)
        vy.append(vy[-1]+oy*dt)
        x.append(x[-1]+vx[-1]*dt)
        y.append(y[-1]+vy[-1]*dt)
        t.append(t[-1]+dt)
    return x,y
fig=plt.figure(figsize=[8,8])
x1=trac(2.5*math.pi)[0]
y1=trac(2.5*math.pi)[1]
plt.plot(x1,y1)
plt.scatter([0],[0],s=1000,color='yellow')
plt.xlabel('x(AU)')
plt.ylabel('y(AU)')
plt.title('beta=2.1,vy0=2.5pi')
plt.show()

演示效果

当初始速度为2pi时，轨道为圆形，轨道方向不转动。初始速度为2.2pi时，轨道为椭圆，轨道方向开始转动。当初始为2.5pi时，轨道椭圆的程度更大，轨道方向的转动速度增加。

总结：

当我们将经典力学运用到宏观物体特别是简单天体物理时，都会得到很好的验证。

参考文献

计算物理 Nicholas J.Giordano, Hisao Nakanishi 清华大学出版社