# 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. 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. 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. Mar 10 '15 at 17:16

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.

• Thank you. These links are good starting points. I think I'll accept your answer and revisit this question when I have completed my error estimates. Mar 10 '15 at 17:14