🌀
楕円軌道を計算する
数年に1度のペースで楕円軌道を計算したくなっては毎回コードをゼロから書いているので、いいかげん書き溜めておこうという記事。
準備
以前にJavascriptで書いたルンゲクッタのコードをPythonに移植
import math
def zip_lists(*args: list[list[float]]) -> list[list[float]]:
min_len = min(len(lst) for lst in args)
return [[lst[i] for lst in args] for i in range(min_len)]
# zip_lists([1,2],[10,20,30],[100,200]) -> [[1, 10, 100], [2, 20, 200]]
def adds(*args: list[list[float]]) -> list[float]:
def sum_list(lst: list[float]) -> float:
return sum(lst)
zipped = zip_lists(*args)
return [sum_list(group) for group in zipped]
# adds([1,2],[10,20,30],[100,200]) -> [111, 222]
def multiple(a: float, b: list[float]) -> list[float]:
return [a * v for v in b]
# multiple(5,[10,20,30]) -> [50, 100, 150]
def estimate_error(x: list[float], x_: list[float], delta: list[float], atol: list[float], rtol: list[float]) -> float:
err = 0.0
for i in range(len(delta)):
sc = atol[i] + max(abs(x[i]), abs(x_[i])) * rtol[i]
err += (delta[i] / sc) ** 2
err = math.sqrt(err / len(delta))
return err
def _ode(
f: callable,
x: list[float],
t: float,
h: float,
) -> tuple[list[float], float, float]:
x_, _ = f(x, t, h)
return x_, t + h, h
def _adaptive_ode(
f: callable,
x: list[float],
t: float,
h: float,
atol: list[float] | float,
rtol: list[float] | float,
maxiter: int,
facmax: float,
facmin: float,
fac: float,
) -> tuple[list[float], float, float]:
"""
f: 関数 f(x: list[float], t: float, h: float, rtol: list[float] | None = None, atol: list[float] | None = None) -> tuple[list[float], float | None]
x: 現在の状態ベクトル
t: 現在の時刻
h: 初期ステップサイズ
Reference
Ernst Hairer, Gerhard Wanner, Syvert P. Nørsett.
Solving Ordinary Differential Equations I: Nonstiff Problems.
Springer.
"""
if isinstance(atol, (int, float)):
atol = [atol] * len(x)
if isinstance(rtol, (int, float)):
rtol = [rtol] * len(x)
h_ = h
i = 0
while i < maxiter:
i += 1
x_, err = f(x, t, h_, atol, rtol)
if err is None:
break
t_ = t + h_
h_ *= min(facmax, max(facmin, fac * (1 / err) ** 0.2))
if err < 1:
return x_, t_, h_
raise RuntimeError("ODE solver failed to converge within max iterations.")
def rkf45(
f: callable,
x: list[float],
t: float,
h: float,
adaptive: bool = False,
atol = 1e-6,
rtol = 1e-3,
maxiter = 100,
facmax = 5.0,
facmin = 0.1,
fac = 0.9,
) -> tuple[list[float], float, float]:
"""
Runge-Kutta-Fehlberg 4(5)法によるODEステップ
f: 微分方程式の関数 f(x, t) -> dx/dt
x: 現在の状態ベクトル
t: 現在の時刻
h: ステップサイズ
"""
def calc(x: list[float], t: float, h: float, rtol: list[float] | None = None, atol: list[float] | None = None) -> tuple[list[float], float | None]:
k1 = f(x, t)
k2 = f(adds(x, multiple(h * 1/4, k1)), t + h * 1/4)
k3 = f(adds(x, multiple(h * 3/32, k1), multiple(h * 9/32, k2)), t + h * 3/8)
k4 = f(adds(x, multiple(h * 1932/2197, k1), multiple(h * -7200/2197, k2), multiple(h * 7296/2197, k3)), t + h * 12/13)
k5 = f(adds(x, multiple(h * 439/216, k1), multiple(h * -8, k2), multiple(h * 3680/513, k3), multiple(h * -845/4104, k4)), t + h)
k6 = f(adds(x, multiple(h * -8/27, k1), multiple(h * 2, k2), multiple(h * -3544/2565, k3), multiple(h * 1859/4104, k4), multiple(h * -11/40, k5)), t + h * 0.5)
x_ = adds(x, multiple(h * 25/216, k1), multiple(h * 1408/2565, k3), multiple(h * 2197/4104, k4), multiple(h * -1/5, k5))
if atol is None or rtol is None:
return x_, None
delta = adds(
multiple(h * 71/57600, k1),
multiple(h * -128/4275, k3),
multiple(h * -2197/75240, k4),
multiple(h * 1/50, k5),
multiple(h * 2/55, k6)
)
err = estimate_error(x, x_, delta, atol, rtol)
return x_, err
if not adaptive:
return _ode(calc, x, t, h)
return _adaptive_ode(calc, x, t, h, atol, rtol, maxiter, facmax, facmin, fac)
def dopri5(
f: callable,
x: list[float],
t: float,
h: float,
adaptive: bool = False,
atol = 1e-6,
rtol = 1e-3,
maxiter = 100,
facmax = 5.0,
facmin = 0.1,
fac = 0.9,
) -> tuple[list[float], float, float]:
"""
Dormand-Prince 5(4)法によるODEステップ
f: 微分方程式の関数 f(x, t) -> dx/dt
x: 現在の状態ベクトル
t: 現在の時刻
h: ステップサイズ
"""
def calc(x: list[float], t: float, h: float, rtol: list[float] | None = None, atol: list[float] | None = None) -> tuple[list[float], float | None]:
k1 = f(x, t)
k2 = f(adds(x, multiple(h * 1/5, k1)), t + h * 1/5)
k3 = f(adds(x, multiple(h * 3/40, k1), multiple(h * 9/40, k2)), t + h * 3/10)
k4 = f(adds(x, multiple(h * 44/45, k1), multiple(h * -56/15, k2), multiple(h * 32/9, k3)), t + h * 4/5)
k5 = f(adds(x, multiple(h * 19372/6561, k1), multiple(h * -25360/2187, k2), multiple(h * 64448/6561, k3), multiple(h * -212/729, k4)), t + h * 8/9)
k6 = f(adds(x, multiple(h * 9017/3168, k1), multiple(h * -355/33, k2), multiple(h * 46732/5247, k3), multiple(h * 49/176, k4), multiple(h * -5103/18656, k5)), t + h)
k7 = f(adds(x, multiple(h * 35/384, k1), multiple(h * 500/1113, k3), multiple(h * 125/192, k4), multiple(h * -2187/6784, k5), multiple(h * 11/84, k6)), t + h)
x_ = adds(x, multiple(h * 35/384, k1), multiple(h * 500/1113, k3), multiple(h * 125/192, k4), multiple(h * -2187/6784, k5), multiple(h * 11/84, k6))
if atol is None or rtol is None:
return x_, None
delta = adds(
multiple(h * 71/57600, k1),
multiple(h * -71/16695, k3),
multiple(h * 71/1920, k4),
multiple(h * -17253/339200, k5),
multiple(h * 22/525, k6),
multiple(h * -1/40, k7)
)
err = estimate_error(x, x_, delta, atol, rtol)
return x_, err
if not adaptive:
return _ode(calc, x, t, h)
return _adaptive_ode(calc, x, t, h, atol, rtol, maxiter, facmax, facmin, fac)
def dopri853(
f: callable,
x: list[float],
t: float,
h: float,
adaptive: bool = False,
atol = 1e-6,
rtol = 1e-3,
maxiter = 100,
facmax = 5.0,
facmin = 0.1,
fac = 0.9,
) -> tuple[list[float], float, float]:
"""
Dormand-Prince method order 8
Reference
dop853.f
http://www.unige.ch/~hairer/software.html
係数は参考文献の数字(30桁)をママ使っている
"""
def calc(x: list[float], t: float, h: float, rtol: list[float] | None = None, atol: list[float] | None = None) -> tuple[list[float], float | None]:
k1 = f(x, t)
k2 = f(adds(x, multiple(5.26001519587677318785587544488e-2 * h, k1)), t + 0.526001519587677318785587544488e-1 * h)
k3 = f(adds(x, multiple(1.97250569845378994544595329183e-2 * h, k1), multiple(5.91751709536136983633785987549e-2 * h, k2)), t + 0.789002279381515978178381316732e-1 * h)
k4 = f(adds(x, multiple(2.95875854768068491816892993775e-2 * h, k1), multiple(8.87627564304205475450678981324e-2 * h, k3)), t + 0.118350341907227396726757197510 * h)
k5 = f(adds(x, multiple(2.41365134159266685502369798665e-1 * h, k1), multiple(-8.84549479328286085344864962717e-1 * h, k3), multiple(9.24834003261792003115737966543e-1 * h, k4)), t + 0.281649658092772603273242802490 * h)
k6 = f(adds(x, multiple(3.7037037037037037037037037037e-2 * h, k1), multiple(1.70828608729473871279604482173e-1 * h, k4), multiple(1.25467687566822425016691814123e-1 * h, k5)), t + 0.333333333333333333333333333333 * h)
k7 = f(adds(x, multiple(3.7109375e-2 * h, k1), multiple(1.70252211019544039314978060272e-1 * h, k4), multiple(6.02165389804559606850219397283e-2 * h, k5), multiple(-1.7578125e-2 * h, k6)), t + 0.25 * h)
k8 = f(adds(x, multiple(3.70920001185047927108779319836e-2 * h, k1), multiple(1.70383925712239993810214054705e-1 * h, k4), multiple(1.07262030446373284651809199168e-1 * h, k5), multiple(-1.53194377486244017527936158236e-2 * h, k6), multiple(8.27378916381402288758473766002e-3 * h, k7)), t + 0.307692307692307692307692307692 * h)
k9 = f(adds(x, multiple(6.24110958716075717114429577812e-1 * h, k1), multiple(-3.36089262944694129406857109825 * h, k4), multiple(-8.68219346841726006818189891453e-1 * h, k5), multiple(2.75920996994467083049415600797e1 * h, k6), multiple(2.01540675504778934086186788979e1 * h, k7), multiple(-4.34898841810699588477366255144e1 * h, k8)), t + 0.651282051282051282051282051282 * h)
k10 = f(adds(x, multiple(4.77662536438264365890433908527e-1 * h, k1), multiple(-2.48811461997166764192642586468 * h, k4), multiple(-5.90290826836842996371446475743e-1 * h, k5), multiple(2.12300514481811942347288949897e1 * h, k6), multiple(1.52792336328824235832596922938e1 * h, k7), multiple(-3.32882109689848629194453265587e1 * h, k8), multiple(-2.03312017085086261358222928593e-2 * h, k9)), t + 0.6 * h)
k11 = f(adds(x, multiple(-9.3714243008598732571704021658e-1 * h, k1), multiple(5.18637242884406370830023853209 * h, k4), multiple(1.09143734899672957818500254654 * h, k5), multiple(-8.14978701074692612513997267357 * h, k6), multiple(-1.85200656599969598641566180701e1 * h, k7), multiple(2.27394870993505042818970056734e1 * h, k8), multiple(2.49360555267965238987089396762 * h, k9), multiple(-3.0467644718982195003823669022 * h, k10)), t + 0.857142857142857142857142857142 * h)
k12 = f(adds(x, multiple(2.27331014751653820792359768449 * h, k1), multiple(-1.05344954667372501984066689879e1 * h, k4), multiple(-2.00087205822486249909675718444 * h, k5), multiple(-1.79589318631187989172765950534e1 * h, k6), multiple(2.79488845294199600508499808837e1 * h, k7), multiple(-2.85899827713502369474065508674 * h, k8), multiple(-8.87285693353062954433549289258 * h, k9), multiple(1.23605671757943030647266201528e1 * h, k10), multiple(6.43392746015763530355970484046e-1 * h, k11)), t + h)
x_ = adds(x, multiple(5.42937341165687622380535766363e-2 * h, k1), multiple(4.45031289275240888144113950566 * h, k6), multiple(1.89151789931450038304281599044 * h, k7), multiple(-5.8012039600105847814672114227 * h, k8), multiple(3.1116436695781989440891606237e-1 * h, k9), multiple(-1.52160949662516078556178806805e-1 * h, k10), multiple(2.01365400804030348374776537501e-1 * h, k11), multiple(4.47106157277725905176885569043e-2 * h, k12))
if atol is None or rtol is None:
return x_, None
delta5 = adds(multiple(0.1312004499419488073250102996e-1, k1), multiple(-0.1225156446376204440720569753e+1, k6), multiple(-0.4957589496572501915214079952, k7), multiple(0.1664377182454986536961530415e+1, k8), multiple(-0.3503288487499736816886487290, k9), multiple(0.3341791187130174790297318841, k10), multiple(0.8192320648511571246570742613e-1, k11), multiple(-0.2235530786388629525884427845e-1, k12))
delta3 = adds(multiple(-0.189800754072407617468755659980, k1), multiple(4.45031289275240888144113950566, k6), multiple(1.89151789931450038304281599044, k7), multiple(-5.8012039600105847814672114227, k8), multiple(-0.422682321323791962932445679177, k9), multiple(-1.52160949662516078556178806805e-1, k10), multiple(2.01365400804030348374776537501e-1, k11), multiple(0.0226517921983608258118062039631, k12))
err5 = 0.0
err3 = 0.0
for i in range(len(x)):
sc = atol[i] + max(abs(x[i]), abs(x_[i])) * rtol[i]
err5 += (delta5[i] / sc) ** 2
err3 += (delta3[i] / sc) ** 2
denominator = err5 + 0.01 * err3 if err5 != 0 or err3 != 0 else 1
err = abs(h) * err5 * (1 / (len(x) * denominator)) ** 0.5
return x_, err
if not adaptive:
return _ode(calc, x, t, h)
return _adaptive_ode(calc, x, t, h, atol, rtol, maxiter, facmax, facmin, fac)
同じような精度の図が出せたので移植成功とする。
import matplotlib.pyplot as plt
gamma = 0.15
def f(x: list[float], t: float) -> list[float]:
return [
x[1],
-2 * gamma * x[1] - x[0]
]
def true_value(t: float) -> float:
return math.exp(-gamma * t) * math.cos(t * math.sqrt(1 - gamma ** 2))
sol_true = []
times = []
sol_rkf45 = []
times_rkf45 = []
sol_dopri5 = []
times_dopri5 = []
sol_dopri853 = []
times_dopri853 = []
x = [1, -0.15]
t = 0
h = 0.2
while t <= 20:
times.append(t)
sol_true.append(true_value(t))
t += h
x = [1, -0.15]
t = 0
h = 0.2
while t <= 20:
sol_rkf45.append(x)
times_rkf45.append(t)
calc = rkf45(f, x, t, h, True)
x = calc[0]
t = calc[1]
h = calc[2]
x = [1, -0.15]
t = 0
h = 0.2
while t <= 20:
sol_dopri5.append(x)
times_dopri5.append(t)
calc = dopri5(f, x, t, h, True)
x = calc[0]
t = calc[1]
h = calc[2]
x = [1, -0.15]
t = 0
h = 0.2
while t <= 20:
sol_dopri853.append(x)
times_dopri853.append(t)
calc = dopri853(f, x, t, h, True)
x = calc[0]
t = calc[1]
h = calc[2]
plt.figure(figsize=(10, 6))
plt.suptitle("Damped Harmonic Oscillator: Numerical vs Analytical Solutions")
plt.subplot(3,1,1)
plt.plot(times, sol_true, label="True Value")
plt.plot(times_rkf45, [x[0] for x in sol_rkf45], "o", linewidth=0, label="RKF45")
plt.grid(True)
plt.legend()
plt.xlim([0,20])
plt.ylabel("x(t)")
plt.subplot(3,1,2)
plt.plot(times, sol_true, label="True Value")
plt.plot(times_dopri5, [x[0] for x in sol_dopri5], "o", linewidth=0, label="DOPRI5")
plt.grid(True)
plt.legend()
plt.xlim([0,20])
plt.ylabel("x(t)")
plt.subplot(3,1,3)
plt.plot(times, sol_true, label="True Value")
plt.plot(times_dopri853, [x[0] for x in sol_dopri853], "o", linewidth=0, label="DOPRI853")
plt.grid(True)
plt.legend()
plt.xlim([0,20])
plt.ylabel("x(t)")
plt.xlabel("Time")
plt.show()
楕円軌道の真値
軌道の長半径、離心率からその時刻（またはその真近点離角での）座標を求める。
ケプラー方程式など軌道力学の式を解説した資料はごまんとあるが、微妙に係数が間違っていたり、記号が資料によって異なったりで、確認しながら正しくプログラムに落とし込むのに毎回大変苦労する。
今回もいくつかの資料[1][2][3]を見比べながら気をつけたつもり。
import math
mue = 398600.4418 # 重力定数 km³/s²
def true_anomaly(E, e):
# 真近点角 true anomaly
# 離心近点離角 eccentric anomaly
theta = 2 * math.atan(math.sqrt((1+e)/(1-e)) * math.tan(E/2)) # -pi/2 ~ pi/2
if theta < 0:
if math.sin(E) > 0:
theta += math.pi
else:
theta += 2*math.pi
elif theta > 0:
if math.sin(E) < 0:
theta += math.pi
return theta # 0 ~ 2pi
def r2d(rad):
return rad * 180 / math.pi
def d2r(deg):
return deg * math.pi / 180
def elliptical_orbit_state_from_time(t:float, a:float, e:float, mu=mue, eps=1e-6, maxIter=10000)-> tuple[float, float, float, float]:
"""
t: 近心点を通過した時刻からの経過時間[s]
a: 長半径[km]
e: 離心率 0<= e < 1
mu: 重力定数[km³/s²]
"""
h = math.sqrt(a * mu * (1 - e**2))
M = mu**2 / h**3 * math.sqrt((1 - e**2)**3) * t
cnt = 0
E = M
while(cnt<maxIter):
if cnt == maxIter:
raise RuntimeError("Kepler's equation failed to converge within max iterations.")
diff = (E - e * math.sin(E) - M) / (1 - e * math.cos(E))
if abs(diff) < eps:
break
E -= diff
cnt += 1
theta = true_anomaly(E, e)
return elliptical_orbit_state_from_theta(theta, a, e, mu)
def elliptical_orbit_state_from_theta(theta:float, a:float, e:float, mu=mue)-> tuple[float, float, float, float]:
"""
theta: 真近点離角 [rad]
a: 長半径[km]
e: 離心率 0<= e < 1
mu: 重力定数[km³/s²]
"""
h = math.sqrt(a * mu * (1 - e**2))
r = h**2 / (mu * (1 + e * math.cos(theta)))
return [r * math.cos(theta), r * math.sin(theta), -mu/h * math.sin(theta), mu/h * (e + math.cos(theta))]
確認のためにひまわり9号のTLEを持ってくる。
1 41836U 16064A 25236.84421777 .00000000 00000-0 00000-0 0 9997
2 41836 0.0245 140.8083 0001231 38.5062 238.7143 1.00267305 32232
ここから離心率と長半径をなんやかんや求めると、e=0.0001231, a=42165.99
e=0.0001231
a=42165.99
rs = []
ts = []
t = 0
while t <= 86400:
ts.append(t)
rs.append(elliptical_orbit_state_from_time(t, a, e))
t += 60*10
rs_ = []
theta = 0
dt = math.pi * 2 / 100
while theta <= math.pi * 2:
rs_.append(elliptical_orbit_state_from_theta(theta, a, e))
theta += dt
時刻を10分刻みで計算した図
真近点離角を2*PI/100radずつ計算した図
楕円軌道の数値解析
def f(x: list[float], t: float) -> list[float]:
r = math.sqrt(x[0]**2 + x[1]**2)
return [
x[2],
x[3],
-mu / r ** 3 * x[0],
-mu / r ** 3 * x[1]
]
def true_value(t: float) -> float:
return elliptical_orbit_state_from_time(t, a, e, mu=mue)
t = 0
h = 60 * 60
r0 = true_value(0)
times = [0]
sol_rkf45 = [r0]
sol_dopri853 = [r0]
sol_true = [r0]
while t <= 86400 * 100:
t += h
sol_rkf45.append(rkf45(f, sol_rkf45[-1], t, h, False)[0])
sol_dopri853.append(dopri853(f, sol_dopri853[-1], t, h, False)[0])
times.append(t)
sol_true.append(true_value(t))
plt.figure(figsize=(10, 6))
plt.plot([r[0] for r in sol_true], [r[1] for r in sol_true], "o", linewidth=0, color="green")
plt.plot([r[0] for r in sol_dopri853], [r[1] for r in sol_dopri853], ">", linewidth=0, color="red")
plt.plot([r[0] for r in sol_rkf45], [r[1] for r in sol_rkf45], "x", linewidth=0, color="blue")
plt.xlim([41250,42250])
plt.ylim([-7500,7500])
plt.grid(True)
plt.show()
軌道の一部を拡大した図だが、時間刻み1時間で100日分ぶんまわしても、真値（緑）から数値解析の値（赤）がほとんどズレていないのがわかる[4]。
離心率の大きな軌道、モルニヤ軌道で確認する。
1 23420U 94081A 25236.34247793 -.00000017 00000-0 00000-0 0 9999
2 23420 64.3239 292.7820 5883563 262.4580 31.1120 3.34495588152639
ここから離心率と長半径をなんやかんや求めると、e = 0.5883563, a = 18886.07
これは10分刻み、24時間分を真値（緑）、数値解析の値（赤）でプロット。
今後やること
fに推力をいられるように改良して強化学習で遊ぶ。
-
https://api.lib.kyushu-u.ac.jp/opac_download_md/4060170/eng2986.pdf ↩︎
-
青は4次のルンゲクッタだが、8次のドルマン＝プリンス法（赤）と比べると酷だが、それでもだいぶ頑張っているとは思う。 ↩︎
Discussion