[Python] Why is the accuracy of scipy.integrate.solve_ivp (RK45) extremely poor compared to my...

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, scipy_ys, 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, my_ys, 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.

Continue reading...