はじめに

Pythonを使って万有引力の運動を解きます。

運動方程式(3体)

万有引力を受けて運動する3物体の運動方程式を以下に示します。
物体をs(sun;太陽), e(earth;地球), l(luna;月)の添え字で区別します。

質量

m_s \\
m_e \\
m_l \\

位置ベクトル

\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}) \\

相対位置ベクトル

\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} \\

速度ベクトル

\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}} \\

運動方程式

\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}} \\

万有引力

\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} \\

万有引力を受けた3物体の運動方程式

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} \\

scipy.integrate.odeの形式

万有引力の運動方程式を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} \\

具体値

太陽、地球、月の質量、距離、速度、および、物理定数は次の通りです。

項目	記号	値	単位
万有引力定数	G	6.67E-11	m^3 kg^-1 s^-2
太陽質量	$m_s$	1.99E+30	kg
地球質量	$m_e$	5.97E+24	kg
月質量	$m_l$	7.35E+22	kg
地球公転半径	$r_{se0}$	1.50E+11	m
月公転半径	$r_{el0}$	3.84E+08	m
地球公転速度	$v_{e0}$	2.98E+04	m/s
月公転速度	$v_{l0}$	1.02E+03	m/s

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()

実行結果