Need help to fully understand SciPy's odeint's reported step sizes, eval times, # of funct calls & total proc. time (re. question in Astronomy SE)

Question

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.

1. But the number of nonzero values in info['hu'] varies from 9 to 274, and I haven't a clue what that means.
2. 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[0].plot(t_eval, x)
        row2[0].plot(t_eval, v)
    row1[0].set_ylabel('position')
    row2[0].set_ylabel('velocity')
    x_finals, v_finals = list(zip(*[result[-1] for result in results]))
    row1[1].plot(x_finals - x_finals[-1])
    row2[1].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()

Answer 1

1. The reported data are all only from the last internal step. There is no internal step performed for the initial point at t_eval[0]. So the output time t_eval[i+1] is between tcur[i]-hu[i] and tcur[i]. There is no data for the other internal steps that do not have an output point inside their interval. So if multiple outputs fit into the same internal step interval, then they all get the same hu and tcur value. You can also check the nst array, the number of internal steps performed so far to get the current output value.

2. Yes. Some basis time for the internal steps plus a time proportional to the number of output values. There might be some non-computing overhead like cache misses that would account for the actual times not exactly conforming to this model.

I do not know why the computation of the internal steps and thus the number of function evaluations is different between the runs. There could be some difference in the translation of the given parameters to the internal parameters of the stepper (DLSODA+RWORK+IWORK), especially from the first output step size to the initial step size of the stepper, or some data is lost and re-initialized between consecutive stepper calls, but logically there should be no difference.