Por que a precisão do scipy.integrate.solve_ivp (RK45) é extremamente baixa em comparação com minha implementação caseira do RK4?

O que eu quero resolver

Resolvi um sistema de equações diferenciais ordinárias usando scipy.integrate.solve_ivpe comparei os resultados com minha implementação caseira de Runge-Kutta de quarta ordem (RK4).

Surpreendentemente, a precisão solve_ivp(usando RK45) é significativamente pior.

Alguém poderia me ajudar a entender por que isso pode estar acontecendo?

Descrição do problema

Simulei um movimento circular uniforme em um plano 2D com uma velocidade inicial de (10, 0) e uma aceleração centrípeta de magnitude 10. Teoricamente, esse sistema deveria descrever um círculo com raio de 10.

Plotei os resultados scipy.integrate.solve_ivp (method="RK45")em azul e aqueles da minha implementação caseira do RK4 em vermelho. O gráfico resultante é mostrado abaixo:

Como o RK4 tem precisão de 4ª ordem e o RK45 tem precisão de 5ª ordem, eu esperava que o RK45 tivesse um desempenho melhor e seguisse o círculo mais de perto. No entanto, os resultados do RK4 são muito superiores.

Código fonte relevante

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()

Sobre o Código Caseiro

O que eu tentei

Li a documentação oficial , solve_ivpmas não consegui descobrir o que pode estar causando essa discrepância.