はじめに
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()
実行結果
