Solving ordinary differential equations with Python ~ Universal gravitation

Introduction

Use Python to solve the universal gravitational movement.

Equation of motion (3 bodies)

The equations of motion of three objects that move under universal gravitation are shown below. Objects are distinguished by the subscripts s (sun), e (earth), and l (luna).

mass

m_s \\
m_e \\
m_l \\

Position vector

\vec{r_s} = (r_{sx}, r_{sy}, r_{sz}) \\
\vec{r_e} = (r_{ex}, r_{ey}, r_{ez}) \\
\vec{r_l} = (r_{lx}, r_{ly}, r_{lz}) \\

Relative position vector

\vec{r_{se}} = \vec{r_e} - \vec{r_s} = - \vec{r_{es}} \\
\vec{r_{sl}} = \vec{r_l} - \vec{r_s} = - \vec{r_{ls}} \\
\vec{r_{el}} = \vec{r_l} - \vec{r_e} = - \vec{r_{le}} \\
r_{se} = | \vec{r_{se}} | = r_{es} \\
r_{sl} = | \vec{r_{sl}} | = r_{ls} \\
r_{el} = | \vec{r_{el}} | = r_{le} \\

Velocity vector

\vec{v_s} = (v_{sx}, v_{sy}, v_{sz}) = \dot{\vec{r_s}} \\
\vec{v_e} = (v_{ex}, v_{ey}, v_{ez}) = \dot{\vec{r_e}} \\
\vec{v_l} = (v_{lx}, v_{ly}, v_{lz}) = \dot{\vec{r_l}} \\

Equation of motion

\vec{F_s} = (F_{sx}, F_{sy}, F_{sz}) = m_s \ddot{\vec{r_s}} = m_s \dot{\vec{v_s}} \\
\vec{F_e} = (F_{ex}, F_{ey}, F_{ez}) = m_e \ddot{\vec{r_e}} = m_e \dot{\vec{v_e}} \\
\vec{F_l} = (F_{lx}, F_{ly}, F_{lz}) = m_l \ddot{\vec{r_l}} = m_l \dot{\vec{v_l}} \\

Universal gravitation

\vec{F_s} = \vec{F_{se}} + \vec{F_{sl}} \\
\vec{F_e} = \vec{F_{es}} + \vec{F_{el}} \\
\vec{F_l} = \vec{F_{ls}} + \vec{F_{le}} \\
\vec{F_{se}} = - \vec{F_{es}} = G * m_s * m_e \frac{\vec{r_{se}}}{r_{se}^3} \\
\vec{F_{sl}} = - \vec{F_{ls}} = G * m_s * m_l \frac{\vec{r_{sl}}}{r_{sl}^3} \\
\vec{F_{el}} = - \vec{F_{le}} = G * m_e * m_l \frac{\vec{r_{el}}}{r_{el}^3} \\

Equation of motion of three objects subject to universal gravitation

m_s * \ddot{\vec{r_s}} = G * m_s * m_e \frac{\vec{r_{se}}}{r_{se}^3} + G * m_s * m_l \frac{\vec{r_{sl}}}{r_{sl}^3} \\
\dot{\vec{v_s}} = \ddot{\vec{r_s}} = G * m_e \frac{\vec{r_{se}}}{r_{se}^3} + G * m_l \frac{\vec{r_{sl}}}{r_{sl}^3} \\

m_e * \ddot{\vec{r_e}} = G * m_e * m_s \frac{\vec{r_{es}}}{r_{es}^3} + G * m_e * m_l \frac{\vec{r_{el}}}{r_{el}^3} \\
\dot{\vec{v_e}} = \ddot{\vec{r_e}} = G * m_s \frac{\vec{r_{es}}}{r_{es}^3} + G * m_l \frac{\vec{r_{el}}}{r_{el}^3} \\
\dot{\vec{v_e}} = \ddot{\vec{r_e}} = - G * m_s \frac{\vec{r_{se}}}{r_{se}^3} + G * m_l \frac{\vec{r_{el}}}{r_{el}^3} \\

m_l * \ddot{\vec{r_l}} = G * m_l * m_s \frac{\vec{r_{ls}}}{r_{ls}^3} + G * m_l * m_e \frac{\vec{r_{le}}}{r_{le}^3} \\
\dot{\vec{v_l}} = \ddot{\vec{r_l}} = G * m_s \frac{\vec{r_{ls}}}{r_{ls}^3} + G * m_e \frac{\vec{r_{le}}}{r_{le}^3} \\
\dot{\vec{v_l}} = \ddot{\vec{r_l}} = - G * m_s \frac{\vec{r_{sl}}}{r_{sl}^3} - G * m_e \frac{\vec{r_{el}}}{r_{el}^3} \\

Format of scipy.integrate.ode

The equation of motion of universal gravitation is expressed in a format handled by scipy.integrate.ode.

\vec{y} = (\vec{r_{s}}, \vec{r_{e}}, \vec{r_{l}}, \vec{v_{s}}, \vec{v_{e}}, \vec{v_{l}}) \\
\dot{\vec{y}} = (\dot{\vec{r_{s}}}, \dot{\vec{r_{e}}}, \dot{\vec{r_{l}}}, \dot{\vec{v_{s}}}, \dot{\vec{v_{e}}}, \dot{\vec{v_{l}}}) \\
\dot{\vec{y}} = (\vec{v_{s}}, \vec{v_{e}}, \vec{v_{l}}, \dot{\vec{v_{s}}}, \dot{\vec{v_{e}}}, \dot{\vec{v_{l}}}) \\

\dot{\vec{v_{s}}} = G * m_e \frac{\vec{r_{se}}}{r_{se}^3} + G * m_l \frac{\vec{r_{sl}}}{r_{sl}^3} \\
\dot{\vec{v_{e}}} = - G * m_s \frac{\vec{r_{se}}}{r_{se}^3} + G * m_l \frac{\vec{r_{el}}}{r_{el}^3} \\
\dot{\vec{v_{l}}} = - G * m_s \frac{\vec{r_{sl}}}{r_{sl}^3} - G * m_e \frac{\vec{r_{el}}}{r_{el}^3} \\

Specific value

The masses, distances, velocities, and physical constants of the Sun, Earth, and Moon are:

Python Program

test.py

import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from scipy.integrate import ode

G = 6.67E-11
ms = 1.99E+30
me = 5.97E+24
ml = 7.35E+22
Gms = G*ms
Gme = G*me
Gml = G*ml
re0 = 1.50E+11
rl0 = 3.84E+08
ve0 = 2.98E+04
vl0 = 1.02E+03

# initial condition
# [0]:rsx,  [1]:rsy,  [2]:rsz,  [3]:rex,  [4]:rey,  [5]:rez, [6]:rlx,  [7]:rly,  [8]:rlz,  
# [9]:vsx, [10]:vsy, [11]:vsz, [12]:vex, [13]:vey, [14]:vez,[15]:vlx, [16]:vly, [17]:vlz,  
y0 = np.array([
  0,    # rsx
  0,    # rsy
  0,    # rsz
  re0,  # rex
  0,    # rey
  0,    # rez
  re0,  # rlx
  0,    # rly
  rl0,  # rlz
  0,    # vsx
  0,    # vsy
  0,    # vsz
  0,    # vex
  ve0,  # vey
  0,    # vez
  vl0,  # vlx
  ve0,  # vly
  0,    # vlz
])

# y
# [0]:rsx,  [1]:rsy,  [2]:rsz,  [3]:rex,  [4]:rey,  [5]:rez, [6]:rlx,  [7]:rly,  [8]:rlz,  
# [9]:vsx, [10]:vsy, [11]:vsz, [12]:vex, [13]:vey, [14]:vez,[15]:vlx, [16]:vly, [17]:vlz,  
# return
# [0]:rsx',  [1]:rsy',  [2]:rsz',  [3]:rex',  [4]:rey',  [5]:rez', [6]:rlx',  [7]:rly',  [8]:rlz',  
# [9]:vsx', [10]:vsy', [11]:vsz', [12]:vex', [13]:vey', [14]:vez',[15]:vlx', [16]:vly', [17]:vlz',  
def f(t, y):
    rse = np.array([y[3]-y[0], y[4]-y[1], y[5]-y[2]]) # sun to earth
    nse = np.linalg.norm(rse)
    rsl = np.array([y[6]-y[0], y[7]-y[1], y[8]-y[2]]) # sun to luna
    nsl = np.linalg.norm(rsl)
    rel = np.array([y[6]-y[3], y[7]-y[4], y[8]-y[5]]) # earth to luna
    nel = np.linalg.norm(rel)
    vsdot =   Gme * rse / nse**3 + Gml * rsl / nsl**3
    vedot = - Gms * rse / nse**3 + Gml * rel / nel**3
    vldot = - Gms * rsl / nsl**3 - Gme * rel / nel**3
    return [y[9], y[10], y[11], y[12], y[13], y[14], y[15], y[16], y[17],
            vsdot[0], vsdot[1], vsdot[2],
            vedot[0], vedot[1], vedot[2],
            vldot[0], vldot[1], vldot[2]]

solver = ode(f)
solver.set_integrator(name="dop853")
solver.set_initial_value(y0)

dt = 60 * 60 * 24
tw = dt * 365 * 2
ts = []
rsx = []
rsy = []
rsz = []
rex = []
rey = []
rez = []
rlx = []
rly = []
rlz = []
while solver.t < tw:
    solver.integrate(solver.t+dt)
    m = ms + me + ml;
    # center of mass
    com = [1/m * (ms * solver.y[0] + me * solver.y[3] + ml * solver.y[6]),
           1/m * (ms * solver.y[1] + me * solver.y[4] + ml * solver.y[7]),
           1/m * (ms * solver.y[2] + me * solver.y[5] + ml * solver.y[8])]
    ts += [solver.t]
    rsx += [solver.y[0] - com[0]]
    rsy += [solver.y[1] - com[1]]
    rsz += [solver.y[2] - com[2]]
    rex += [solver.y[3] - com[0]]
    rey += [solver.y[4] - com[1]]
    rez += [solver.y[5] - com[2]]
    rlx += [solver.y[6] - com[0]]
    rly += [solver.y[7] - com[1]]
    rlz += [solver.y[8] - com[2]]

plt.figure(0)
plt.axes(projection='3d')
plt.plot(rsx, rsy, rsz, "r-+")
plt.plot(rex, rey, rez, "b-+")
plt.plot(rlx, rly, rlz, "g-+")
plt.show()

Execution result