Pythonで常微分方程式を解く~万有引力

# はじめに
Pythonを使って万有引力の運動を解きます。

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

質量

```math
m_s \\
m_e \\
m_l \\
```

位置ベクトル

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

相対位置ベクトル

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

速度ベクトル

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

運動方程式

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

万有引力

```math
\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物体の運動方程式

```math
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で扱う形式で表現します。

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

```math
\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

```py3: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()
```

実行結果

![figure_0.png](https://qiita-image-store.s3.amazonaws.com/0/84205/364fbd33-560e-9655-7cb8-0a12412f79db.png)