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