Quantify integration error of scipy ode / ODEPACK

Question

I am trying to integrate a 2nd order ODE with potential several singularities using the lsoda solver wrapped in scipy.integrate.ode(). I would like to put an error bar on the solution or at least estimate of upper bound to the error.

For understand how to calculate truncation errors for simple finite difference methods. But LSODA seems to be more complicated from the scipy documentation, the solver uses a adaptive stepsize and automatically switches between a stiff and non-stiff method.

Options to set the absolute tolerance, relative tolerance and maximum step size exist but I don't understand exactly how they would relate to the error in my solution from this document.

There are a few cases where analytic solutions are possible. e.g. a case that I debugged in a previous question.

The equation is Newcomb's Euler-Lagrange equation from the field of plasma physics.

$$\frac{d}{dr}(f \frac{d\xi}{dr}) - g \xi = 0$$ or as a set of first order ODE's: $$y_0 = \xi $$ $$y_1 = f \xi' $$ $$y_0' = \frac{y_1}{f} $$ $$y_1' = y_0 g $$

$f$ and $g$ are complicated expressions of magnetic fields and pressure gradients. For the more complicated cases, I use splines generated with scipy.interpolate.InterpolateUnivariateSpline() to describe the magnetic fields and pressure gradients. $f$ and $g$ can vary rapidly and $f$ can be zero at several locations resulting in singularities. I use Frobenius expansions to find the solutions close to singularties.

Here is a short example code for a simple case of $f$ and $g$, where there is only a singularity at $r=0$.

import numpy as np
from scipy.integrate import ode

# Setup ODE system
def f(r):
   return r

def g(r):
    return -1 + r + 1./r

def der(r, y):
    y_der = np.zeros(2)
    y_der[0] = y[1]/f(r)
    y_der[1] = g(r)*y[0]
    return y_der

#Integrate
integrator = ode(der)
integrator.set_integrator('lsoda') 
r_init = 1E-3
init = [r_init, 1.]
integrator.set_initial_value(init, t=r_init)
r = np.linspace(r_init, 1., 100)
results = np.zeros((2, r.size))
results[:,0] = init
for i, position in enumerate(r[1:]):
    integrator.integrate(position)
    results[:, i+1] = integrator.y

Answer 1

If I remember correctly, you have some options:

• You could use an integrator that calculates rigorous bounds via interval arithmetic. For instance, Taylor models are such an approach, and you could look at work by Berz and Makino. You'll need to use a different integrator, and probably change your code substantially to enable use of interval arithmetic data types (and operations).

• You could estimate the error using a posteriori approaches, of which there are several. This preprint by Emil Constantinescu has a comprehensive literature review. One of these approaches (a reference by Robert Skeel in Emil's paper reviews thirteen methods, so there are others) is to derive a defect equation that approximates the error, and integrate that in time, so you have to solve your ODE at least twice.

• Based on global error estimates, you could actually try to control the global error in the integration. Cao and Petzold describe such an approach using DASPK. In order to use this method, you need to use a solver capable of adjoint sensitivity analysis. Cao and Petzold use DASPK, but you could also use CVODES. This is one instance of the approach DavidKetcheson remarked on.