第十次作业

行星运动

---

摘要

在本次作业中，我们把目光投向太阳系，考虑在太阳系中行星运动的相关问题。与课本所述相互印证，探究对于变形的万有引力定律和水星进动的相关规律。

背景

考虑牛顿的万有引力定律

$$F_G=\frac{GM_SM_E}{r^2}$$
从牛顿第二定律有 $$\frac{d^2x}{d^2t}=\frac{F_{G，x}}{M_E}$$ $$\frac{d^2y}{d^2t}=\frac{F_{G，y}}{M_E}$$
利用AU（天文单位，天文学中测量距离，特别是测量太阳系内天体之间的距离的基本单位，地球到太阳的平均距离为一个天文单位。一天文单位约等于1.496亿千米），方程可以写为： $$GM_S=v^2r=4π^2\frac{AU^3}{yr^2}$$
故可得到如下微分方程组
$$\begin{equation} \left\{ \begin{aligned} \overset{.}v_{x，i+1}=v_{x，i}-\frac{4π^2x_i}{r_i^3}Δt \\ x_i+1=x_i+v_{x，i+1}Δt \\v_{y，i+1}=v_{y，i}-\frac{4π^2y_i}{r_i^3}Δt\\y_i+1=y_i+v_{y，i+1}Δt \end{aligned} \right. \end{equation}$$
可以假设万有引力的形式如下 $$F_G=\frac{GM_SM_E}{r^β}$$
对轨道偏离最大的水星，其引力可近似为 $$F_G≈\frac{GM_SM_E}{r^2}（1+\frac{α}{r^2}）$$

正文

一.在不同β下的椭圆轨道模拟

import matplotlib.pyplot as pl
class PlanetOrbit:
    def __init__(self,time_step=0.001,init_x=1,init_y=0,init_vx=0,init_vy=1.3*math.pi,
        beta=2.10,T_a=0.998,eccentricity=0.017,total_time=2):
        self.Kep = 4*pow(math.pi,2)/T_a
        self.e = eccentricity
        self.dt = time_step
        self.r = []
        self.x = [init_x]
        self.y = [init_y]
        self.vx = [init_vx]
        self.vy = [init_vy]
        self.t = total_time
        self.power = beta + 1
    def Run(self):
        loop = True
        i = 0
        while(loop):
            self.r.append(pow(((self.x[i])**2+(self.y[i])**2),0.5))
            self.vx.append(self.vx[i]-(self.Kep*self.x[i]*self.dt)/(self.r[i]**self.power))
            self.x.append(self.x[i]+self.vx[i+1]*self.dt)
            self.vy.append(self.vy[i]-(self.Kep*self.y[i]*self.dt)/(self.r[i]**self.power))
            self.y.append(self.y[i]+self.vy[i+1]*self.dt)
            i += 1
            if (i>=(self.t/self.dt-1)):
                loop = False
    def DrawOrbit(self):
        pl.plot(self.x,self.y,label="beta = $self.beta$")
        pl.title('planet orbiting the Sun')
        pl.xlabel('x(AU)')
        pl.ylabel('y(AU)')
        pl.axis('equal' )
        pl.legend()
        pl.show()
a=PlanetOrbit()
a.Run()
a.DrawOrbit()

结果如下




β=3.00时，如下

二.在不同α下的水星轨道模拟（在进动条件下远日点的变化）

import matplotlib.pyplot as plt
import numpy as np
import math
alfa=0.01
x0=0
y0=0
xp=[]
yp=[]
tp=[]
theta=[]
class elliptical:
    def __init__(self,GMs=4*math.pi**2,dt=0.0001,time=2):
        self.GMs=GMs
        self.x=[0.47]
        self.y=[0]
        self.vx=[0]
        self.vy=[8.2]
        self.dt=dt
        self.time=time
        self.r=[math.sqrt(self.x[0]**2+self.y[0]**2)]
        self.t=[0]
    def calculate(self):
        for i in range(int(self.time//self.dt)):
            self.vx.append(self.vx[i]-self.GMs*self.x[i]*self.dt/self.r[i]**3-alfa*self.GMs*self.x[i]*self.dt/self.r[i]**5)
            self.vy.append(self.vy[i]-self.GMs*self.y[i]*self.dt/self.r[i]**3-alfa*self.GMs*self.y[i]*self.dt/self.r[i]**5)
            self.x.append(self.x[i]+self.vx[i+1]*self.dt)
            self.y.append(self.y[i]+self.vy[i+1]*self.dt)
            self.t.append(self.t[i]+self.dt)
            self.r.append(math.sqrt(self.x[i+1]**2+self.y[i+1]**2))
    def precession(self):
        for j in range(len(self.r)-2):
            if (self.r[j+1]**2-self.r[j]**2>0 and self.r[j+1]**2-self.r[j+2]**2>0):
                xp.append(self.x[j+1])
                yp.append(self.y[j+1])
                tp.append(self.t[j+1])
                if self.x[j+1]>=0 and self.y[j+1]/self.x[j+1]>=0:
                    theta.append(math.atan(self.y[j+1]/self.x[j+1])*180/math.pi)
                if self.x[j+1]>=0 and self.y[j+1]/self.x[j+1]<0:
                    theta.append(math.atan(self.y[j+1]/self.x[j+1])*180/math.pi+360)
                if self.x[j+1]<0 and self.y[j+1]/self.x[j+1]>=0:
                    theta.append(math.atan(self.y[j+1]/self.x[j+1])*180/math.pi+180)
                if self.x[j+1]<0 and self.y[j+1]/self.x[j+1]<0:
                    theta.append(math.atan(self.y[j+1]/self.x[j+1])*180/math.pi+180)
    def show_results(self):
        for j in range(len(xp)):
            plt.plot([0,xp[j]],[0,yp[j]])
        plt.plot(self.x,self.y,'g',label=r'$\alpha$=0.0008'+' trajectory')
        plt.scatter(x0,y0)
        plt.title(r'$\alpha$=0.0008'+'  Simulation of the precession of Mercury ',fontsize=14)
        plt.xlabel(u'x(AU)',fontsize=14)
        plt.ylabel(u'y(AU)',fontsize=14)
        plt.legend(fontsize=14,loc='best')
a=elliptical()
a.calculate()
a.precession()
a.show_results()
plt.show()

结果如下
远日点


3.Vpython

结论

1.β越大，轨道的进动越明显
2.α越大，远日点的变化越明显

致谢

Professor Cai
The book of Computational Physics, chapter 4
Wikipedia