「楕円軌道上の時刻ごとの位置を求めるための下ごしらえ」で下ごしらえした式を Python で数値的に解く。
ライブラリの import
import numpy as np
import matplotlib.pyplot as plt
# 以下はグラフを SVG で Notebook にインライン表示させる設定
%config InlineBackend.figure_formats = ['svg']
無次元化された式
「楕円軌道上の時刻ごとの位置を求めるための下ごしらえ」で下ごしらえした式（無次元化済み）。
scipy.integrate.solve_ivp() は $\displaystyle \frac{dy}{dt} = f(t, y)$ の形の式を解くので，変数を揃える。
$ \phi \rightarrow y, \ T \rightarrow t$
\begin{eqnarray}
\frac{dy}{dt} &=& f(t, y) = 2 \pi \frac{(1+e \cos y)^2}{(1-e^2)^{3/2}} \tag{1}\\
R(y) &\equiv& \frac{r}{a} = \frac{1-e^2}{1+e\cos y} \tag{2}\\
X(y) &\equiv& \frac{x}{a} = R(y) \cos y \tag{3}\\
Y(y) &\equiv& \frac{y}{a} = R(y) \sin y \tag{4}
\end{eqnarray}
# 微分方程式の右辺 独立変数は t, y 決め打ち
def f(t, y):
global e
return 2*np.pi*(1+e*np.cos(y))**2/((1-e**2)**(3/2))
def R(y):
global e
return (1-e**2)/(1+e*np.cos(y))
def X(y):
global e
return R(y)*np.cos(y)
def Y(y):
global e
return R(y)*np.sin(y)
# 離心率
e = 0.6
solve_ivp() で常微分方程式を数値的に解く
from scipy.integrate import solve_ivp
# T の初期値
t0 = 0
# T の終了値
t1 = 1
t_span = [t0, t1]
# phi の初期値
y_ini = [0]
# 分割数
Ndiv = 36
t_list = np.linspace(t0, t1, Ndiv+1)
answer = solve_ivp(f, t_span, y_ini,
t_eval = t_list,
# 計算精度の設定
rtol = 1.e-13, atol = 1.e-15)
# 計算結果は .y として取り出します。
phii = answer.y[0]
ちょうど1周期 $T = T_{\rm end}$ たてば $\phi = 2 \pi$ になっているはずだから，数値誤差を確認する。
$$\phi(T_{\rm end}) – 2 \pi$$
# [-1] は最後の項
phii[-1] - 2*np.pi
一定時間間隔ごとの位置をグラフに
一定時間間隔 $\displaystyle \Delta T = \frac{1}{N_{\rm div}}$ ごとの位置：
fig = plt.figure(figsize=[6,6])
# 縦横のアスペクト比を等しく
plt.axes().set_aspect('equal')
# 軸の設定
plt.axis([-2, 1, -1.5, 1.5])
plt.scatter(X(phii), Y(phii), c="blue", s=20);
# 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法則を視覚的に確認する」にまとめているように，楕円軌道の $x, y$ を軌道長半径 $a$ で規格化した $X_u, Y_u$ を，離心近点離角 $u$ を使って以下のようにあらわすことができる。
\begin{eqnarray}
X_u &\equiv& \frac{x}{a} = \cos u – e\\
Y_u &\equiv& \frac{y}{a} = \sqrt{1-e^2} \sin u
\end{eqnarray}
離心近点離角 $u$ と周期 $P$ で規格化された時間 $T$ との関係は，以下のケプラー方程式を数値的に解いて得られる。
\begin{eqnarray}
2\pi T &=& u – e \sin u, \quad T \equiv \frac{t}{P}
\end{eqnarray}
ケプラー方程式の数値解
$$g(u_i) \equiv u_i – e \sin u_i – 2\pi T_i =0, \quad T_i = \frac{i}{N_{\rm div}}$$
を $u_i$ について root_scalar() 関数で数値的に解きます。
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つの解の比較
$\phi$ に関する1階常微分方程式を数値的に解いた解である X(phii), Y(phii) と，ケプラー方程式を数値的に解いた解である Xu(ui), Yu(ui) は，当然ながら数値誤差の範囲で一致してるはずである。
2つの配列の引き算をして差を調べる。37組全て表示すると冗長なので，代表して3組ほど。
# 2つの数値解 X(phii), Xu(ui) の差
print((X(phii) - Xu(ui))[10])
print((X(phii) - Xu(ui))[20])
print((X(phii) - 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(phii), Y(phii), c="blue", s=20, 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();
