为什么 scipy.integrate.solve_ivp (RK45) 的准确性与我自制的 RK4 实现相比非常差？

Question

我想解决什么问题

我使用

scipy.integrate.solve_ivp

求解了常微分方程组，并将结果与我自制的 4 阶 Runge-Kutta (RK4) 实现进行了比较。

令人惊讶的是，

solve_ivp

（使用RK45）的准确性明显较差。

有人可以帮助我理解为什么会发生这种情况吗？

问题描述

我在 2D 平面中模拟了匀速圆周运动，初始速度为 (10, 0)，向心加速度为 10。 理论上，这个系统应该描述一个半径为 10 的圆。

我将

scipy.integrate.solve_ivp (method="RK45")

的结果绘制为蓝色，将我自制的 RK4 实现的结果绘制为红色。结果图如下所示：

由于RK4具有4阶精度，RK45具有5阶精度，因此我预计RK45会表现更好，并且更紧密地跟随圆圈。然而，RK4 的结果要优越得多。

相关源码

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

def get_a(t, r, v):
    # Accelerates perpendicular to the direction of motion with magnitude 10
    x, y = 0, 1
    
    direction_of_motion = math.atan2(v[y], v[x])
    direction_of_acceleration = direction_of_motion + math.pi / 2
    acceleration_magnitude = 10

    a_x = acceleration_magnitude * math.cos(direction_of_acceleration)
    a_y = acceleration_magnitude * math.sin(direction_of_acceleration)

    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])
    # returns    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)

    return np.concatenate([v, a])

# Simulation settings
initial_position = [0, 0]
initial_velocity = [10, 0]
end_time = 300
time_step = 0.01

# scipy.integrate.solve_ivp simulation (RK45)
initial_values = np.array(initial_position + initial_velocity)
time_points = np.arange(0, end_time + time_step, time_step)
result = solve_ivp(
    fun=get_v_and_a,
    t_span=[0, end_time],
    y0=initial_values,
    method="RK45",
    t_eval=time_points
)

scipy_ts = result.t
scipy_xs = result.y[0]
scipy_ys = result.y[1]

# Homemade RK4 simulation
my_ts, my_xs, my_ys = [], [], []
current_time = 0.0
step_count = 0
r = np.array(initial_position, dtype="float64")
v = np.array(initial_velocity, dtype="float64")
delta_t = time_step

while current_time <= end_time:
    my_ts.append(current_time)
    my_xs.append(r[0])
    my_ys.append(r[1])
    
    t = current_time + delta_t    
    r_and_v = np.concatenate([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)

    step_count += 1
    current_time = step_count * 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]

# Plot results
# 1. set up
size = max([abs(x) for x in scipy_xs] + [abs(y) for y in scipy_ys])*2.5
fig, ax = plt.subplots(figsize=(10, 10))
ax.set_xlim(-size/2, size/2)
ax.set_ylim(-size/2, size/2)
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_aspect("equal", adjustable="box")

# 2. plot
ax.plot(scipy_xs, scipy_ys, lw=0.1, color="blue",
    label="scipy.integrate.solve_ivp RK45")
ax.plot(my_xs, my_ys, lw=0.1, color="red", label="homemade code RK4")

# 3. time points
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. show
leg = ax.legend()
leg.get_lines()[0].set_linewidth(1)
leg.get_lines()[1].set_linewidth(1)
ax.grid(True)
plt.show()

关于自制代码

我尝试过的

我已阅读solve_ivp的

官方文档

，但无法弄清楚是什么导致了这种差异。

Answer 1

如果降低

solve_ivp()

的容差，我会得到更好的结果。

例如

result = solve_ivp(
    fun=get_v_and_a,
    t_span=[0, end_time],
    y0=initial_values,
    method="RK45",
    rtol=1e-5,
    atol=1e-6,
    t_eval=time_points
)

rtol默认值为1e-3，更改为1e-5可以使模拟更加准确。 rtol 值越小，solve_ivp 越准确，但代价是需要对函数进行更多评估。

这是 rtol=1e-5 的样子：

如果没有此更改，自制积分器会对您的函数进行 120,000 次调用，但 SciPy 仅进行 1,226 次调用。自制积分器以更小的时间步长模拟您的函数，因此最终会更加准确。将 rtol 降低到 1e-5 需要 3,512 次调用，因此您可以在更短的时间内获得与自制解决方案相当的精度。