Exercise_10 : Precession of the Perihelion of Mercury

Problem 4.10 & 4.11

冉峰 弘毅班 2014301020064

---

---

Abstract

• From this charpter, we learn the problem of Kepler's Laws. And with the help of the Kepler's Laws, first we can solve the motion of the simplest situation where a sun and a single planet get moved with interaction. We investigate a few of the properities of this model solar system. Additionally,the inverse-square law and the stability of plantery orbits will be discussed. Then we can learn the precession of the perihelion of Mercury and we can try to plot the figure of the precession of perihelion of Mercury to have a better learning. Also, just like the better we use the Euler-Cromer method to solve the trobule of the precession of perihelion of Mercury.

Background

Kepler's Laws

• According to the Newton's law of the gravitation the magnitude of this force is given by $$F_G=\frac{GM_SM_E}{r^2}$$
  From the Newton's secnd law, we have $$\begin{eqnarray*} \frac{d^2x}{dt^2}=\frac{F_{G,x}}{M_E}\quad\frac{d^2y}{dt^2}=\frac{F_{G,y}}{M_E} \end{eqnarray*}$$
  Thus $$F_{G,x}=-\frac{GM_SM_E}{r^2}cos\theta=-\frac{GM_SM_Ex}{r^2}$$
  We can follow our usual approach and write each of the second-order differential equations as two first-order differential equations.
  $$\begin{eqnarray*} \frac{dv_x}{dt}=-\frac{GM_S x}{r^3} \quad \frac{dx}{dt}=v_x \end{eqnarray*}$$
  $$\begin{eqnarray*} \frac{dv_y}{dt}=-\frac{GM_S y}{r^3} \quad \frac{dy}{dt}=v_y \end{eqnarray*}$$
  Worth mentioning, the selection of $SI$ is of great importance. We here chooe the distance $r$ between earth and sun to be $1AU$. And the time unit to be $1yr$. Thus we obtain:
  $$\begin{eqnarray*} GM_S=v^2r=4\pi^{2}(AU^3/yr^2) \end{eqnarray*}$$
  We can use Euler-Cromer method to solve the equations:
  $$r_i=(x_i^2+y_i^2)^{1/2}$$
  $$v_{x,i+1}=v_{x,i}-\frac{4\pi^{2}x_i}{r_i^{3}}\Delta t\qquad x_{i+1}=x_i+v_{x,i+1}\Delta t$$
  $$v_{y,i+1}=v_{y,i}-\frac{4\pi^{2}y_i}{r_i^{3}}\Delta t\qquad y_{i+1}=y_i+v_{y,i+1}\Delta t$$
  The orbital trajectory for a body of reduced mass $\mu$ is given in polar coordinates by $$\frac{d^2}{d\theta^2}(\frac{1}{r})+\frac{1}{r}=-\frac{\mu r^2}{L^2}F(r)$$
  Also the momentum L is a constant,
  Thus we have $$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})}$$
  

The inverse-square law and the stability of plantery orbits

• Also if a force: $$F_G=\frac{GM_SM_E}{r^{\beta}}$$
  Where $\beta$ is not always 2. So we can consider the problem with $\beta$ is different from 2.
  When $\beta$ isdifferent from 2, the elleptical orbits with this force law can be simulated with the planteary motion program given above by simply changing the exponent of $r$ in the equations for the velocity.
  

The precession of the perihelion of Mercury

• Perhaps you may be impressed with a simple superficial situation that all the precise movemet of planet can be calclated by the previous program precisely. Nevertheless, there are deviations from these laws.These deviations come from a number of sources, including the effects of the planets on each other. It turns out that this is a problem for which very few exact results are known.
  Of most importance is the precession of the perihelion of Mercury. The planets whose orbits deviate the most from circular are Mercury and Pluto.
  All we have to do is simulate the orbital motion using the force law predicted by the general relativity and measure the rate of precession of the orbit.
  $$\begin{eqnarray*} F_G\approx\frac{GM_S M_M}{r^2}(1+\frac{\alpha}{r^2}) \end{eqnarray*}$$
  And $$-\frac{GM_SM_M}{r_1}+\frac{1}{2}M_Mv_1^2=-\frac{GM_SM_M}{r_2}+\frac{1}{2}M_Mv_2^2$$
  Where $r_1v_1=bv_2$ and $b=a\sqrt{1-e^2}$
  

---

The Main Body

Code

Created on Sat Nov 26 18:04:56 2016
@author: RanFeng
import math
import matplotlib.pyplot as plt
GM=4*(math.pi**2)
perihelion=0.39*(1-0.206)
class Precession:
    def __init__(self,e=0.206,time=1,dt=0.0001,vcoefficient=1.1,beta=2.1,alpha=0.01):
        self.alpha=alpha
        self.beta=beta
        self.vcoefficient=vcoefficient
        self.e=e
        self.a=perihelion/(1-e)
        self.x0=self.a*(1+e)
        #self.x0=1
        self.y0=0
        self.vx0=0
        self.vy0=self.vcoefficient*math.sqrt((GM*(1-self.e))/(self.a*(1+self.e)))
        #self.vy0=2*math.pi
        self.X=[self.x0]
        self.Y=[self.y0]
        self.Vx=[self.vx0]
        self.Vy=[self.vy0]
        self.T=[0]
        self.dt=dt
        self.time=time
        self.ThetaPrecession=[]
        self.TimePrecession=[]
    def run(self):
        while self.T[-1]<self.time:
            r=math.sqrt(self.X[-1]**2+self.Y[-1]**2)
            newVx=self.Vx[-1]-(GM*(1+self.alpha/r**2)*self.X[-1]/r**(1+self.beta))*self.dt
            newX=self.X[-1]+newVx*self.dt
            newVy=self.Vy[-1]-(GM*(1+self.alpha/r**2)*self.Y[-1]/r**(1+self.beta))*self.dt
            newY=self.Y[-1]+newVy*self.dt
            if abs(newX*newVx+newY*newVy)<0.0014 and r<self.a:
                theta=math.acos(self.X[-1]/r)*180/math.pi
                if (self.Y[-1]/r)<0:
                    theta=360-theta
                theta=abs(theta-180)
                self.ThetaPrecession.append(theta)
                self.TimePrecession.append(self.T[-1])
            self.Vx.append(newVx)
            self.Vy.append(newVy)
            self.X.append(newX)
            self.Y.append(newY)
            self.T.append(self.T[-1]+self.dt)
    def show_results(self):
        plt.plot(self.X,self.Y,'g')
        plt.legend(loc='upper right',frameon=False,fontsize='small')
        plt.grid(True)
        plt.title('Simulation of the precession of Mercury')
        plt.xlabel('x(AU)')
        plt.ylabel('y(AU)')
        ax = plt.gca()
        ax.set_aspect(1)
        #plt.xlim(-1,1)
        #plt.ylim(-0.6,0.6)
        #plt.scatter(0,0)
A = Precession()
A.run()
A.show_results()

Result

• First we can have a figure to see the normal result for Earth orbiting around the Sun.

• Then we can plot the figure that with elliptical oribit.
  When $\beta=2$
  

• Also we can chose different $\beta$, such as $\beta=2.01,2.1,2.5,3$.
  

  $$\beta=2.01$$

$$\beta=2.1$$

$$\beta=2.5$$

$$\beta=3$$

---

• Then we can investigate the precession of Mercury.
  First we can chose $\alpha=0.01$,$\beta=2$ and $dt=0.0001yr$
  
• Also we can chose different $\alpha$, as $\alpha=0.1,0.25,0.0008,0.05$
  $$\alpha=0.0008$$
  $$\alpha=0.05$$
  $$\alpha=0.25$$
  $$\alpha=0.1$$
• Just like former, we can chose the different $\beta$ when $\alpha=0.01$.
  $$\beta=2.01$$

$$\beta=2.1$$

$$\beta=2.5$$

$$\beta=3$$

• We can see that it's a little different from $\alpha=0$

---

• What's more, we can plot the orbit orientaion versus time
  when $\alpha=0.0008$, $\beta=2$
  
• Also with different $\alpha=0.01$, $\beta=2$
  

• With different $\alpha=0.01$, $\beta=2.01$
  

• It's easy to know the difference.

---

• Now we will try to investigate the different ellipyical oribit with different eccentricities. We shall see with bigger eccentricity, the trajecory becomes more elletical, as suggested by the definition of eccentricity.
  First we chose $\alpha=0.0008$, $\beta=2$ and $e=0.206$
  
• Also we chose $e=0.306,0.406,0.506$
  $$e=0.306$$
  $$e=0.406$$
  $$e=0.506$$

• we can plot them together.
  

• We can see with bigger eccentricity, the trajecory becomes more elletical, as suggested by the definition of eccentricity.

---

• We can chose $\alpha=0.01$, $\beta=2$ and $e=0.206,0.306,0.406,0.506$

---

• Finally,we can investigate the influence of different $v_0$.
  We chose $\alpha=0.0008$, $\beta=2$ and $e=0.206$
  and $v=1.2,1.4,1.5,1.7v_0$
  $$v=1.2v_0$$

$$v=1.4v_0$$

$$v=1.5v_0$$

$$v=1.7v_0$$

• Plot them together
  

• We chose $v=1.1v_0,\beta=2.1,\alpha=0.0008$
  

• We chose $v=1.1v_0,\beta=2.1,\alpha=0.01$

---

VPython code

from visual import *
#from Numeric import *
print '''
Simple solar system model showing the orbits of the planets
Note: In this version the orbits are not tilted so
      the longitude of the ascending node is not correct.
created by Lensyl Urbano, Jan. 2006.
'''
class orbit:
    def __init__(self,e=0.0,rad=1.0, mu=10.0, planet_size=1.0, period=365,
                 inclination=0.0, lap=0.0, frate=10000000):
        self.orbit_path = curve(color=color.green,radius=0.075)
        self.planet = sphere(color=color.yellow,radius=planet_size)
        if e >= 1.0:
            e = 0.99
        if e < 0.0:
            e = 0.0
        fx = 2.0 - 2.0*min(abs((1-e*e)/(1-e)), abs(-(1-e*e)/(1+e)))
        self.fx = fx
        self.e = e
        self.rad = rad
        self.dangle = 360./period
        self.inclination = inclination
        self.lap = lap
        for i in range(int(period)+1):
            rate(frate)
            theta = radians(float(360.*i/period))
            #h_sq = a*(1-e*e)*mu
            r = (1-e*e)/(1-e*cos(theta))
            nx = rad*(r*cos(theta) - fx)
            ny = rad*r*sin(theta) 
            #print i, f, r, cos(theta), sin(theta), nx, ny
            #nz = rotate(vector(nx,ny,0.0), angle=inclination, 
            posi = rotate(vector(nx,ny,0.0), angle=radians(inclination), axis=(0,1,0))
            posi = rotate(posi, angle=radians(lap), axis=(0,0,1))
            self.orbit_path.append(pos=posi)
            self.planet.pos=posi
        self.angle = 0.0
    def re_draw(self,e=0.0,rad=1.0,mu=10.0, frate=10000000):
        #redraw last path
        if e >= 1.0:
            e = 0.99
        if e < 0.0:
            e = 0.0
        fx = 2.0 - 2.0*min(abs((1-e*e)/(1-e)), abs(-(1-e*e)/(1+e)))
        self.fx = fx
        self.e = e
        self.rad = rad
        for i in range(361):
            rate(frate)
            theta = radians(float(i))
            #h_sq = a*(1-e*e)*mu
            r = (1-e*e)/(1-e*cos(theta))
            nx = rad*(r*cos(theta) - fx)
            ny = rad*r*sin(theta) 
            self.orbit_path.pos[i]=vector(nx, ny, 0.0)
    def move_planet(self):
        self.angle = self.angle + self.dangle
        #to show the planet orbiting the Sun
        theta = radians(float(self.angle))
        r = (1-self.e*self.e)/(1-self.e*cos(theta))
        nx = self.rad*(r*cos(theta) - self.fx)
        ny = self.rad*r*sin(theta) 
        posi=rotate(vector(nx,ny,0.0), angle=radians(self.inclination), axis=(0,1,0))
        posi = rotate(posi, angle=radians(self.lap), axis=(0,0,1))
        self.planet.pos = posi
#class proton:
origx = 0
origy = 0
w = 704 #+4+4
h = 576 #+24+4
scene.width=w
scene.height=h
scene.x = origx
scene.y = origy
mu = 10.0 #gravitational acceleration
#e = 0.6 #eccentricity
#a = 10.0 #distance of sun from center of elipse
# to draw an elipse
nucleus = frame()
p1 = sphere(color=color.white, frame=nucleus)
#orbitor = sphere()
e1 = orbit(e=0.50, rad=15.0, mu=mu, 
              planet_size=0.5, period=365,
              lap=90)
##e2 = orbit(e=0.0, rad=15.0, mu=mu, 
##              planet_size=1.0, period=365,
##              lap=0, inclination=90)
runtime = 0.0
framerate = 60
pick = None
while 1:
    rate(framerate)
    runtime += 1.0/framerate
    e1.move_planet()
##    e2.move_planet()

Result

---

Conclusion

• With the help of Kepler's Laws and the Euler-Cromer method, we can solve the motion of the simplest situation where a sun and a single planet get moved with interaction. We investigate a few of the properities of this model solar system. Then we can add the inverse-square law and the stability of plantery orbits will be discussed. Thus e can learn the precession of the perihelion of Mercury. And we can try to plot the figure of the precession of perihelion of Mercury. Thus we can know the precession of perihelion of Mercury more clear.

---

Thanks For

• The book of Computational Physics, chapter 4
• Wang Shixing
• Wikipedia
• Vpython from ZZT