「万有引力の2体問題の運動方程式を数値的に解く前の下ごしらえ ver. 2」で下ごしらえした式を Python で数値的に解く。
無次元化された運動方程式
「万有引力の2体問題の運動方程式を数値的に解く前の下ごしらえ ver. 2」にまとめたように，軌道長半径 $a$ と周期 $P$ で無次元化した変数に対する運動方程式は以下のようになるのであった。
\begin{eqnarray}
\frac{d^2 X}{dT^2} = – 4\pi^2\frac{X}{\left(X^2 + Y^2\right)^{\frac{3}{2}}} \\
\frac{d^2 Y}{dT^2} = – 4\pi^2\frac{Y}{\left(X^2 + Y^2\right)^{\frac{3}{2}}}
\end{eqnarray}
無次元化された初期条件
無次元化された変数に対する初期条件は $T = 0$ で
\begin{eqnarray}
X(0) &=& 1-e \\
Y(0) &=& 0\\
V_x(0) &=& 0 \\
V_y(0) &=& 2\pi\sqrt{\frac{1+e}{1-e}}
\end{eqnarray}
この初期条件で数値計算すると，（規格化された）長半径 $1$，離心率 $e$ の楕円の軌道が得られるハズ。
じゃあ，そこまでわかっているのなら，なぜわざわざ数値計算するのかというと，それは楕円軌道の解が時間 $t$ の陽関数として与えられているわけではないから。時刻 $t$ のとき，どこにいるかが解析的にわかっていないので，それを知りたいから数値計算することになる。
ライブラリの import
import numpy as np
import matplotlib.pyplot as plt
# 以下はグラフを SVG で Notebook にインライン表示させる設定
%config InlineBackend.figure_formats = ['svg']
Runge-Kutta 用に連立1階微分方程式系に
\begin{eqnarray}
\frac{dX}{dT} &=& F_1(V_x) = V_x \\
\frac{dY}{dT} &=& F_2(V_y) = V_y \\
\frac{dV_x}{dT} &=& F_3(X, Y)
= -4\pi^2 \frac{X}{\left(X^2 + Y^2\right)^{\frac{3}{2}}} \\
\frac{dV_y}{dT} &=& F_4(X, Y)
= -4\pi^2 \frac{Y}{\left(X^2 + Y^2\right)^{\frac{3}{2}}}
\end{eqnarray}
from scipy.integrate import solve_ivp
# scipy.integrate.solve_ivp() は
# dy/dt = Fun(t, y) の形を解く。変数名は決め打ち。
def Fun(t, y):
# solve.ivp() に渡す関数の引数の順番に注意。t が先
X = y[0]
Y = y[1]
Vx = y[2]
Vy = y[3]
F1 = Vx
F2 = Vy
F3 = -4*np.pi**2 * X/(X**2 + Y**2)**(3/2)
F4 = -4*np.pi**2 * Y/(X**2 + Y**2)**(3/2)
return [F1, F2, F3, F4]
# 離心率 e = 0.6
e = 0.6
t0 = 0
tend = 1
t_span = [t0, tend]
# 初期条件 X0, Y0, Vx0, Vy0
X0 = 1-e
Y0 = 0
Vx0 = 0
Vy0 = 2*np.pi*np.sqrt((1+e)/(1-e))
y_ini = [X0, Y0, Vx0, Vy0]
Ndiv = 36
t_list = np.linspace(t0, tend, Ndiv+1)
sol = solve_ivp(Fun, t_span, y_ini,
t_eval = t_list,
rtol = 1.e-12,
atol = 1.e-14)
ちょうど1周期 $T=t/P = 1$ までいくと，$T=0$ の時の初期値になっているはずだから，数値誤差を確認する。
数値解 sol から X すなわち y[0] を取り出すには，sol.y[0]。最後の項は sol.y[0][-1]。最初の項（初期条件）sol.y[0][0] との差を確認する。
X = sol.y[0]
Y = sol.y[1]
# X の誤差
print(X[-1]-X[0])
# Y の誤差
print(Y[-1]-Y[0])
一定時間間隔ごとの位置をグラフに
# グラフのサイズの設定
fig = plt.figure(figsize=[6,6])
# 縦横のアスペクト比を等しくして
plt.axes().set_aspect('equal')
# 軸の設定
plt.axis([-2, 1, -1.5, 1.5])
plt.scatter(sol.y[0], sol.y[1], s=20, c='blue')
# x軸 y軸
plt.axhline(0, color='black', dashes=(3, 3), linewidth=0.6)
plt.axvline(0, color='black', dashes=(3, 3), linewidth=0.6);
ケプラー方程式を数値的に解いた場合の解との比較
「ケプラー方程式を数値的に解いてケプラーの第2法則を視覚的に確認する」にまとめているように，
\begin{eqnarray}
x &=& r \cos\varphi = a(\cos u – e)\\
y &=& r \sin\varphi = a \sqrt{1-e^2} \sin u\\
\frac{2\pi}{P} t &=& u – e \sin u
\end{eqnarray}
だから，無次元化した変数では
\begin{eqnarray}
X_u(u) &=& \frac{x}{a} = \cos u – e\\
Y_u(u) &=& \frac{y}{a} = \sqrt{1-e^2} \sin u\\
2\pi T &=& u – e \sin u, \quad T = \frac{t}{P}
\end{eqnarray}
ケプラー方程式の数値解
$$g(u) \equiv u – e \sin u$$
として $T_i = \frac{i}{N_{\rm div}}$ のときにケプラー方程式 $g(u_i) = 2\pi T_i $ を数値的に解いて $u_i$ を求める。
from scipy.optimize import root_scalar
# root_scalar を使うときは，変数名は x 決め打ち
def g(x):
global e, T
return x - e*np.sin(x) - 2*np.pi*T
Tlist = np.linspace(0, 1, Ndiv+1)
ui = []
for T in Tlist:
ans = root_scalar(g, bracket = [0, 2*np.pi], xtol = 1.e-14)
ui.append(ans.root)
一定時間間隔ごとの位置
def Xu(u):
global e
return np.cos(u) - e
def Yu(u):
global e
return np.sqrt(1-e**2) * np.sin(u)
2つの解の比較
運動方程式を数値的に解いた解である X, Y とケプラー方程式を数値的に解いた解である Xu(ui), Yu(ui) は，当然ながら数値誤差の範囲内で一致しているはずである。
2つのリストの引き算をして差を調べる。37組全て表示すると冗長なので，代表して3組ほど。
# 2つの数値解 X, Xu(ui) の差
print((X - Xu(ui))[10])
print((X - Xu(ui))[20])
print((X - Xu(ui))[30])
念のために2つの数値解をグラフにすると，当然ながら重なっていることがわかる。
# グラフのサイズの設定
fig = plt.figure(figsize=[6,6])
# 縦横のアスペクト比を等しくして
plt.axes().set_aspect('equal')
# 軸の設定
plt.axis([-2, 1, -1.5, 1.5])
plt.scatter(X, Y, s=20, c='blue',label = "運動方程式の数値解")
plt.scatter(Xu(ui), Yu(ui), c="red", s=2, label = "ケプラー方程式の数値解")
# x軸 y軸
plt.axhline(0, color='black', dashes=(3, 3), linewidth=0.6)
plt.axvline(0, color='black', dashes=(3, 3), linewidth=0.6)
plt.legend();
