# Quantify integration error of scipy ode / ODEPACK

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

• If you want a bound on the global error, I believe it can't be done with LSODA or most ODE solver codes. They only control local error. To bound global error you need a method that will go back and repeat the whole integration adaptively. – David Ketcheson Mar 6 '15 at 2:38
• Check out the Assimulo package. It has a nice interface and allows to use CVODE(S) as solver. For now only Python2 but on its way to Python3. – GertVdE Mar 9 '15 at 13:46
• @GertVdE thanks for the tip. I see that Assimulo also wraps the sensitivity options in CVODES. That will be useful for the third approach in Geoff's answer. – jensv Mar 10 '15 at 17:16