scipy の solve_ivp の精度が不自然に悪い (rk45)
解決したいこと
常微分方程式を用意し、 scipy の solve_ivp と 自作の 4次ルンゲクッタコードを比較したのですが、 solve_ivp の精度が極端に悪いです。この原因について分かる方いたら教えていただけますでしょうか。
発生している問題
2次元空間で初速度 (10,0), 向心加速度の大きさ10 で等速円運動のシミュレーションを行いました。
理論的には、これは半径10の円を描くはずです。
solve_ivp (method="RK45") による結果を青、
自作コード (RK4) による結果を赤でプロットしたところ、次の画像のようになりました。
rk4は4次精度、rk45は5次精度なので、rk45のほうがより半径10の円に近づくはずかと思うのですが、実際には rk4 のほうがはるかに良い結果を得ております。
該当するソースコード
from scipy.integrate import solve_ivp
import math
import matplotlib.pyplot as plt
import numpy as np
def get_a(t, r, v):
# 進行方向に垂直な方向に10だけ加速し続ける
x,y = 0,1
進行方向 = math.atan2(v[y], v[x])
加速度の方向 = 進行方向 + math.pi / 2
加速度の大きさ = 10
a_x = 加速度の大きさ * math.cos(加速度の方向)
a_y = 加速度の大きさ * math.sin(加速度の方向)
return np.array([a_x, a_y])
def get_v_and_a(t, r_and_v):
# r_and_v == np.array([r_x, r_y, v_x, v_y])
# v_and_a == np.array([v_x, v_y, a_x, a_y])
r = r_and_v[0:2]
v = r_and_v[2:4]
a = get_a(t, r, v)
v_and_a = np.concat([v,a])
return v_and_a
# シミュレーションの設定
初期位置 = [0, 0]
初速度 = [10, 0]
終了時刻 = 300
時刻幅 = 0.01
# シミュレーションの実行 (scipy.integrate.solve_ivp rk45)
初期値 = np.array(初期位置 + 初速度)
時点たち = np.arange(0, 終了時刻+時刻幅, 時刻幅)
結果 = solve_ivp(
fun = get_v_and_a,
t_span = [0, 終了時刻],
y0 = 初期値,
method = "RK45",
t_eval = 時点たち
)
scipy_ts = 結果.t
scipy_xs = 結果.y[0]
scipy_ys = 結果.y[1]
# シミュレーションの実行 (自作コード rk4)
my_ts, my_xs, my_ys = [], [], []
時刻 = 0.
ステップ = 0
r = np.array(初期位置, dtype="float64")
v = np.array(初速度, dtype="float64")
delta_t = 時刻幅
while 時刻 <= 終了時刻:
my_ts += [時刻]
my_xs += [r[0]]
my_ys += [r[1]]
t = 時刻 + delta_t
r_and_v = np.concat([r,v])
k0 = get_v_and_a(t - delta_t, r_and_v )
k1 = get_v_and_a(t - 1/2 * delta_t, r_and_v + 1/2 * k0 * delta_t)
k2 = get_v_and_a(t - 1/2 * delta_t, r_and_v + 1/2 * k1 * delta_t)
k3 = get_v_and_a(t , r_and_v + k2 * delta_t)
ステップ += 1
時刻 = ステップ * delta_t
r_and_v += (1/6 * k0 + 1/3 * k1 + 1/3 * k2 + 1/6 * k3)*delta_t
r = r_and_v[0:2]
v = r_and_v[2:4]
# 結果のグラフ
# 1. 用意
幅 = max([abs(x) for x in scipy_xs] + [abs(y) for y in scipy_ys])*2.5
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111)
plt.xlim(-幅/2, 幅/2)
plt.ylim(-幅/2, 幅/2)
plt.xlabel("X")
plt.ylabel("Y")
ax.set_aspect("equal", adjustable="box")
# 2. 軌跡
ax.plot(scipy_xs, scipy_ys, lw = 0.1, color = "blue")
ax.plot(my_xs, my_ys, lw = 0.1, color = "red")
# 3. 時刻
for i, t in enumerate(scipy_ts):
if t % 5 == 0:
ax.text(
scipy_xs[i], scipy_ys[i], str(round(t)), fontsize=8,
ha = "center", va = "center", color = "blue"
)
if t % 1 == 0 and (t <= 5 or 295 <= t):
ax.text(
my_xs[i], my_ys[i], str(round(t)), fontsize=8,
ha = "center", va = "center", color = "red"
)
# 4. 表示
leg = ax.legend(["scipy.integrate.solve_ivp rk45", "homemade code rk4"])
leg.get_lines()[0].set_linewidth(1)
leg.get_lines()[1].set_linewidth(1)
ax.grid(True)
plt.show()
自作コードについて
$$\begin{array}{l}
\overrightarrow{{\rm r\_and\_v}(t)}=\overrightarrow{{\rm r\_and\_v}(t-\Delta t)}
+ \left(
\frac{1}{6} \overrightarrow{k_0} +
\frac{1}{3} \overrightarrow{k_1} +
\frac{1}{3} \overrightarrow{k_2} +
\frac{1}{6} \overrightarrow{k_3}
\right)\Delta t\\
\overrightarrow{k_0} = \overrightarrow{{\rm get\_v\_and\_a}\left(
t-\Delta t,
\overrightarrow{{\rm r\_and\_v}(t-\Delta t)}
\right)}\\
\overrightarrow{k_1} = \overrightarrow{{\rm get\_v\_and\_a}\left(
t-\frac{1}{2}\Delta t,
\overrightarrow{{\rm r\_and\_v}(t-\Delta t) + \frac{1}{2}\overrightarrow{k_0}\Delta t}
\right)}\\
\overrightarrow{k_2} = \overrightarrow{{\rm get\_v\_and\_a}\left(
t-\frac{1}{2}\Delta t,
\overrightarrow{{\rm r\_and\_v}(t-\Delta t) + \frac{1}{2}\overrightarrow{k_1}\Delta t}
\right)}\\
\overrightarrow{k_3} = \overrightarrow{{\rm get\_v\_and\_a}\left(
t,
\overrightarrow{{\rm r\_and\_v}(t-\Delta t) + \overrightarrow{k_2}\Delta t}
\right)}\\
\end{array}$$
のように更新式を定義しています。ただし$\overrightarrow{{\rm r\_and\_v}(t)}$ は、時刻 $t$ における位置と速度をまとめた4次元ベクトルであり、 $\overrightarrow{{\rm get\_v\_and\_a}\left(
t,
\overrightarrow{{\rm r\_and\_v}(t)}
\right)}$ は、$\frac{\rm d}{{\rm d}t}\overrightarrow{{\rm r\_and\_v}(t)}$ であり、位置と速度をまとめた4次元ベクトルを、速度と加速度をまとめた4次元ベクトルに変換する関数です。
(自作コードについて 終わり)
自分で試したこと
公式ドキュメント を読んでみました。