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.