Cumulative interpolation error of ODE solvers vs. regular error

interpolation_error.py

from scipy.integrate import solve_ivp
from scipy.integrate._ivp.ivp import METHODS

import numpy as np

method = "Radau"
times = np.linspace(0,500,1000)
atol = 1e-10
rtol = 1e-5

analytical_solution = np.vstack([
		np.sin(times),
		np.cos(times),
	])

ODE = lambda t,y: np.array([y[1],-y[0]])
initial = [0,1]

normal_solution = solve_ivp(
		fun = ODE,
		t_span = [times[0],times[-1]],
		y0 = initial,
		method = method,
		dense_output = True,
		atol=atol, rtol=rtol,
	).sol(times)

stepwise_solution = [initial]
for i in range(1,len(times)):
	integrator = METHODS[method](
			fun = ODE,
			t0 = times[i-1],
			y0 = stepwise_solution[i-1],
			t_bound = np.inf,
			atol=atol, rtol=rtol,
		)
	while integrator.t < times[i]:
		integrator.step()
	stepwise_solution.append(integrator.dense_output()(times[i]))

stepwise_solution = np.vstack(stepwise_solution).T

import matplotlib.pyplot as plt

plt.plot(
		times,
		normal_solution[0]-analytical_solution[0],
		label='regular error'
	)
plt.plot(
		times,
		stepwise_solution[0]-analytical_solution[0],
		label='regular error and cumulative interpolation error'
	)
plt.xlabel('t')
plt.ylabel('error')
plt.legend()
plt.show(block=False)