计算物理第十次作业

计算物理

---

摘要

我们将研究行星的轨道运动，探究其影响因素，并应用于研究水星近日点进动，同时解答习题4.8和4.10

背景

　　
　　太阳系，作为我们探索外太空的第一步，一直以来备受我们关注，而对他的探索也从未停止。为何唯独这看似普通的第三颗行星上会出现智慧生命，还有没有未被探测到的新行星，太阳系内还有那些适宜生命生存的地方。而了解这一切的起源，都基于我们对太阳系最初的探索，即行星轨道的观测和计算。

1. 基本运动公式

　　根据我们之前学过的知识，即牛顿定律，我们知道地球和太阳间的引力遵循如下公式：
　　

$$F_G=\frac{GM_SM_E}{r^2}$$
　　其中$M_S$和$M_E$分别为太阳和地球的质量，$r$为日地距离，$G$为万有引力常数
　　因此我们有：
$$\begin{eqnarray*} \frac{d^2x}{dt^2}&=&\frac{F_{G,x}}{M_E}\\ \frac{d^2y}{dt^2}&=&\frac{F_{G,y}}{M_E}\\ \end{eqnarray*}$$
　　考虑到$F_{G,x}=-\frac{GM_SM_E}{r^2}cos\theta=-\frac{GM_SM_Ex}{r^3}$，则有：
$$\begin{eqnarray*} \frac{dv_x}{dt}&=&-\frac{GM_Sx}{r^3}\\ \frac{dx}{dt}&=&v_x\\ \frac{dv_y}{dt}&=&-\frac{GM_Sy}{r^3}\\ \frac{dy}{dt}&=&v_y \end{eqnarray*}$$
　　又因为：$GM_S=v^2r=4\pi AU^3/yr^2$，使用欧勒法则有：
$$\begin{eqnarray*} v_{x,i+1}&=&v_{x,i}-\frac{4\pi^{2}x_i}{r_i^{3}}\Delta t\\ x_{i+1}&=&x_i+v_{x,i+1}\Delta t\\ v_{y,i+1}&=&v_{y,i}-\frac{4\pi^{2}x_i}{r_i^{3}}\Delta t\\ y_{i+1}&=&y_i+v_{y,i+1}\Delta t \end{eqnarray*}$$
　　

2.开普勒三定律：

　　
　　1.所有行星绕太阳的轨道都是椭圆，太阳在椭圆的一个焦点上。
　　2.行星和太阳的连线在相等的时间间隔内扫过相等的面积。
　　3.所有行星绕太阳一周的恒星时间(T )的平方与它们轨道长半轴(a)的比值是一个常数，即$\frac{T^2}{a^3}=c$
　　

3. 平方反比定律和轨道稳定性　　

　　极坐标下的运动方程：

$$\begin{eqnarray*} \frac{d^2}{d\theta^2}(\frac{1}{r})+\frac{1}{r}&=&-\frac{\mu r^2}{L^2}F(r)\\ \end{eqnarray*}$$
　　考虑$L=\mu r^2\dot \theta$及$F(r)=-GM_SM_P/r^2$我们可以解出：
$$\begin{eqnarray*} \frac{1}{r}&=&（\frac{\mu GM_SM_P}{L^2})[1-ecos(\theta+\theta_0)]\\ r&=&(\frac{L^2}{\mu GM_SM_P})\frac{1}{1-ecos\theta} \end{eqnarray*}$$ 　　
　　其中$e$为偏心率，当$e=0$时为圆，$0<e<1$时为椭圆，$e=1$时为抛物线，$e>1$时为双曲线。
　　我们研究行星运动，故只考虑椭圆，则有$r_{min}=a(1-e)$,$r_{max}=a(1+e)$, 且$a=L^2/[\mu GM_SM_P(1-e^2)]$则有
　　 $$\begin{eqnarray*} 　　v_{max}&=&\sqrt{GM_S}\sqrt{\frac{(1+e)}{a(1-e)}(1+\frac{M_P}{M_S})}\\ 　　v_{min}&=&\sqrt{GM_S}\sqrt{\frac{(1-e)}{a(1+e}(1+\frac{M_P}{M_S})} 　　\end{eqnarray*}$$
　　
　　我们不禁会考虑，若引力不是平方反比会怎样，因此对于
　　 $$F_G=\frac{GM_SM_E}{r^{\beta}}$$
我们可以考虑改变$\beta$的值来观察轨道的变化

3. 水星近日点进动

　　根据相对论，我们修正万有引力

$$F_G\approx \frac{GM_SM_M}{r^2}(1+\frac{\alpha}{r^2})$$
　　
　　

正文

1. 主体程序

import pylab as pl
import math
class Solar_system :
    def __init__(self,i=0,initial_x=1,initial_y=0,initial_vx=0,initial_vy=2*math.pi,\
                 total_time=10, time_step=0.002, radius=1):
        self.x=[initial_x]
        self.y=[initial_y]
        self.vx=[initial_vx]
        self.vy=[initial_vy]
        self.r=[radius]
        self.time=total_time
        self.dt=time_step
        self.t=[0]
    def run(self):
        _time=0
        while(_time<self.time):
            self.vx.append(self.vx[-1]-4*pow(math.pi,2)*self.x[-1]/pow(self.r[-1],3)*self.dt)
            self.x.append(self.x[-1]+self.vx[-1]*self.dt)
            self.vy.append(self.vy[-1]-4*pow(math.pi,2)*self.y[-1]/pow(self.r[-1],3)*self.dt)
            self.y.append(self.y[-1]+self.vy[-1]*self.dt)
            self.t.append(_time)            
            _time += self.dt
    def show_results(self):
        pl.plot(self.x,self.y)
        pl.title('Solar system')
        pl.xlabel('x(AU)')
        pl.ylabel('y(AU)')
        pl.legend()
        pl.show()
a = Solar_system()
a.run()
a.show_results()  

2. 平方反比定律和轨道稳定性

　　我们为了研究非平方反比定律下行星轨道，因此这里需要修改万有引力公式，即

$$F_G=\frac{GM_SM_E}{r^{\beta}}$$
　　因此在程序中需要修改

    self.r.append(math.sqrt(pow(self.x[-1],2)+pow(self.y[-1],2)))
    self.vx.append(self.vx[-1]-4*pow(math.pi,2)*self.x[-1]/pow(self.r[-1],1+self.beta)*self.dt)
    self.x.append(self.x[-1]+self.vx[-1]*self.dt)
    self.vy.append(self.vy[-1]-4*pow(math.pi,2)*self.y[-1]/pow(self.r[-1],1+self.beta)*self.dt)
    self.y.append(self.y[-1]+self.vy[-1]*self.dt)

3.水星近日点进动

　　这里需要考虑相对论效应，修改万有引力公式：
　　

$$F_G \approx \frac{FM_SM_E}{r^2}(1+\frac{\alpha}{r^2})$$
　　故程序修改如下

    self.r.append(math.sqrt(pow(self.x[-1],2)+pow(self.y[-1],2)))
    self.vx.append(self.vx[-1]-4*pow(math.pi,2)*(1+self.alpha/pow(self.r[-1],2))*self.x[-1]/pow(self.r[-1],1+self.beta)*self.dt)
    self.x.append(self.x[-1]+self.vx[-1]*self.dt)
    self.vy.append(self.vy[-1]-4*pow(math.pi,2)*(1+self.alpha/pow(self.r[-1],2))*self.y[-1]/pow(self.r[-1],1+self.beta)*self.dt)
    self.y.append(self.y[-1]+self.vy[-1]*self.dt)

结论

1.开普勒定律

基本图像（圆轨道）

椭圆轨道

我们可以通过修改e的值来影响轨道形状：

即当$e=0$时为圆，$0<e<1$时为椭圆，$e=1$时为抛物线，$e>1$时为双曲线。

　　当然，我们也可以通过改变初始速度达到目的：
　　
　　这也是航天器变轨的基本规律
　　

2.平方反比定律

　　若我们将$r$的指数改变，会引起一系列变化
　　为了对比方便，我们将坐标轴范设为相同，且打点画图：

3.近日点进动

　　这里我们考虑实际的参量，即$a=0.39AU$,$e=0.206$,$v_1=8.2AU$,$r_1=0.47AU$，因此我们调整坐标系及参量绘图如下

　　分析$\theta$和$t$的关系：
　
　　我们可以改变$\alpha$的值，从而得到不同的结果：
　　
　　因此我们还可以做出$\alpha$和$d\theta/dt$的图像如下：
　　
　　
　　

4.Vpython

　　当然，我们还可以用Vpython对运动进行模拟，以更好的了解行星的运动

致谢

张梓桐同学帮忙指正了一些小bug