A recent question in Astronomy SE Numerical Programming using odeint takes more than 17 minutes got me interested in looking closer at SciPy's odeint.
The problem is a modified orbital mechanical problem in the solar system. I'd used odeint with great success in orbital mechanics without really understanding the details of how it works under the hood. I chose a reasonably small rtol, enough evaluation times for the output to make nice plots and if necessary rescaled from billions of kilometers and seconds to distances and times of order unity in AU and years.
I had always assumed that the integrator chose its own step sizes internally, and when it was all finished re-evaluated the results at the requested times by interpolation.
odeint kept the integration results and interpolation coefficients hidden from the user, and the newer solve_ivp exposes these if you set
dense_output=True. This will allow one to choose a specific region of interest and do a fine-grid interpolation at a later time.
It also allows for a bigger choice of integration routines.
In the script below, I get the same basic answer (1E-08 variation in final positions and velocities) whether I request the output evaluate at 10 time steps or a million.
In all cases there's about 500 function evaluations, which makes sense considering the final results are the same.
- But the number of nonzero values in
info['hu']varies from 9 to 274, and I haven't a clue what that means.
- Total process time increases from 5 to 750 milliseconds. Is that due to the time spent interpolating being 100x longer than the integration for 1 million output values?
I wrote a short script (included below) for a 1D strongly anharmonic oscillator, here are the results:
import numpy as np import matplotlib.pyplot as plt from scipy.integrate import solve_ivp, odeint import time def deriv(y, t, k1, k2): x, v = y a = -k1 * x - k2 * x**3 return np.hstack((v, a)) k1, k2 = 0.5, 2 y0 = np.array([1, 0], dtype=float) args = (k1, k2) n_steps = [10**i + 1 for i in range(1, 8)] rtol = 1E-10 atol = None t_evals, results, infos, process_times = , , ,  for n_step in n_steps: t_eval = np.linspace(0, 10, n_step) time_start = time.process_time() result, info = odeint(deriv, y0, t_eval, args=args, full_output=True, rtol=rtol, atol=atol) results.append(result) infos.append(info) process_time = time.process_time() - time_start process_times.append(process_time) t_evals.append(t_eval) print(n_step, process_time * 1E+06) if True: fig, (row1, row2) = plt.subplots(2, 2) for t_eval, result in zip(t_evals, results): x, v = result.T row1.plot(t_eval, x) row2.plot(t_eval, v) row1.set_ylabel('position') row2.set_ylabel('velocity') x_finals, v_finals = list(zip(*[result[-1] for result in results])) row1.plot(x_finals - x_finals[-1]) row2.plot(v_finals - v_finals[-1]) plt.show() tot_func_evals = [info['nfe'][-1] for info in infos] eval_times = [np.array(info['tcur']) for info in infos] d_eval_times = [et[1:] - et[:-1] for et in eval_times] milli_seconds = [1000 * t for t in process_times] non_zero_d_eval_times = [[d for d in d_eval_time if d != 0] for d_eval_time in d_eval_times] if True: fig, ax = plt.subplots(1, 1) things = non_zero_d_eval_times, n_steps, tot_func_evals, milli_seconds for thing, n_step, n_func_eval, ms in zip(*things): non_zeros = len(thing) ax.plot(thing, label=str((n_step, non_zeros, n_func_eval, round(ms, 1)))) ax.set_title('non-zero changes in evaluation times') ax.legend() plt.show()