楕円軌道上の時刻ごとの位置を Python で数値的に求める

「楕円軌道上の時刻ごとの位置を求めるための下ごしらえ」で下ごしらえした式を Python で数値的に解く。

ライブラリの import

In [1]:
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}

In [2]:
# 微分方程式の右辺 独立変数は 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() で常微分方程式を数値的に解く

In [3]:
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$$

In [4]:
# [-1] は最後の項

phii[-1] - 2*np.pi
Out[4]:
-2.5348612098241574e-12

一定時間間隔ごとの位置をグラフに

一定時間間隔 $\displaystyle \Delta T = \frac{1}{N_{\rm div}}$ ごとの位置：

In [5]:
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() 関数で数値的に解きます。

In [6]:
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)

一定時間間隔ごとの位置

In [7]:
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組ほど。

In [8]:
# 2つの数値解 X(phii), Xu(ui) の差
print((X(phii) - Xu(ui))[10])
print((X(phii) - Xu(ui))[20])
print((X(phii) - Xu(ui))[30])

3.0744295997919835e-12
-4.904965322793942e-13
-1.5972778655282127e-12

念のために2つの数値解をグラフにすると，当然ながら重なっていることがわかる。

In [9]:
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();