4-5次ルンゲ＝クッタ法とオイラー法の比較

ロボットアーム制御の勉強をしていたところ、ルンゲ＝クッタ法なる言葉がでてきました。
将来位置・速度の予想という事でオイラー法は何となく知ってましたが……ルンゲ＝クッタ法は言葉も思い出せず(実際やっていないかもしれない)。
というわけで、どうも一般的らしい4-5次ルンゲ＝クッタ法とオイラー法の比較したプロットを作ってみました。

実際のところ

オイラー法の方が速いらしいので、可視化するために時間も計測してみました。

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
import time

# 微分方程式を定義
def fun(t, y):
    return [y[1], -np.sin(y[0])]

# 初期値
y0 = [np.pi - 0.1, 0.0]

# 時間範囲を定義
t_span = [0, 10]

# solve_ivpを用いて問題を解く (4-5次ルンゲ＝クッタ法)
start = time.time()
sol = solve_ivp(fun, t_span, y0, method='RK45', dense_output=True)
end = time.time()
rk45_time = end - start

# オイラー法で同様の問題を解く
t = np.linspace(t_span[0], t_span[1], 1000)
y = np.empty((2, t.size))
y[:, 0] = y0

dt = t[1] - t[0]
start = time.time()
for i in range(1, t.size):
    y[:, i] = y[:, i-1] + dt * np.array(fun(t[i-1], y[:, i-1]))
end = time.time()
euler_time = end - start

# 結果をプロット
plt.figure(figsize=(12, 6))

plt.subplot(1, 2, 1)
plt.plot(sol.t, sol.y[0], label=f'RK45 (time: {rk45_time:.6f} sec)')
plt.plot(t, y[0], '--', label=f'Euler (time: {euler_time:.6f} sec)')
plt.title("y[0] comparison")
plt.legend()

plt.subplot(1, 2, 2)
plt.plot(sol.t, sol.y[1], label=f'RK45 (time: {rk45_time:.6f} sec)')
plt.plot(t, y[1], '--', label=f'Euler (time: {euler_time:.6f} sec)')
plt.title("y[1] comparison")
plt.legend()

plt.show()