Why is the accuracy of scipy.integrate.solve_ivp (RK45) extremely poor compared to my homemade RK4 implementation?

What I want to solve

I solved a system of ordinary differential equations using scipy.integrate.solve_ivp and compared the results with my homemade 4th-order Runge-Kutta (RK4) implementation.

Surprisingly, the accuracy of solve_ivp (using RK45) is significantly worse.

Could someone help me understand why this might be happening?

Problem Description

I simulated uniform circular motion in a 2D plane with an initial velocity of (10, 0) and a centripetal acceleration of magnitude 10. Theoretically, this system should describe a circle with a radius of 10.

I plotted the results from scipy.integrate.solve_ivp (method="RK45") in blue and those from my homemade RK4 implementation in red. The resulting plot is shown below:

Since RK4 has 4th-order accuracy and RK45 has 5th-order accuracy, I expected RK45 to perform better and follow the circle more closely. However, the results from RK4 are far superior.

Relevant Source Code

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

About the Homemade Code

What I Tried

I have read the official documentation for solve_ivp but couldn't figure out what might be causing this discrepancy.

I get substantially better results if I lower the tolerances on solve_ivp().

e.g.

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
)

The default value of rtol is 1e-3, and changing the value to 1e-5 makes the simulation more accurate. Smaller values of rtol make solve_ivp more accurate, at the cost that it takes more evaluations of your function.

Here's what rtol=1e-5 looks like:

Without this change, the homemade integrator makes 120,000 calls to your function, but SciPy makes only 1,226 calls. The homemade integrator is simulating your function at a much smaller timestep, so it ends up being more accurate. Solving this problem with an rtol of 1e-5 requires 3,512 calls, so you can get accuracy comparable to the homemade solution in less time.