Large residual when integrating 2nd order ode close to singularity with SciPy ode / ODEPACK

I am trying to integrate a 2nd order ODE with a singularity at close to the initial condition. Why do I get large residuals when I plug-in the result of my integration back into the ODE?

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. However, I even have trouble for simple cases, such as: $$f=r$$ $$g=-1 + r + \frac{1}{r}$$

The ODE is always singular at $r=0$ From Fobenius expansion the solution close to $r=0$ for this case should be $\xi \propto r^1 \approx 0$ and $\xi' \propto r^0 \approx 1$.

In the code below I set up the problem and integrate with scipy.integrate.ode using the lsoda solver. I plug the result back into the 2nd order ODE form: $$f \xi'' + f' \xi' -g \xi$$ I have to use numpy.gradient to calculate $\xi''$.

I get residuals on the edge order of order of $10^2$. Why are they so large? If I use constants for $f$ and $g$ I get residuals on the order of my step size.

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

# Print residual
dr = r[1] - r[0]


The output I get looks like the array below. There are large residuals at the edge.

array([  5.70419912e+02,  -1.59592700e+02,  -9.00242498e-01,
-1.50959245e-01,  -4.59057150e-02,  ...,
-1.13773852e-03,  -1.11566938e-03,  -1.07162899e-03,
-1.79645207e+00])

• You are computing the "residual" with a low order discretization and large steps. Your solution was computed with a high order discretization and small steps. – David Ketcheson Feb 5 '15 at 4:53
• You mean that lsoda is using an adaptive step size to achieve the default tolerance, probably much smaller than the 0.01 step size I am using to calculate the residual? How would you suggest checking the result? – jensv Feb 5 '15 at 5:29
• Why are you calculating the residual as $f (\xi'/f)'+f'\xi'-g\xi$? – Kirill Feb 5 '15 at 6:37
• results[1] is $f \xi'$. I am calculating $f (\frac{\xi' f}{f})'+f' \xi' - g \xi$ – jensv Feb 5 '15 at 7:06
• I think the edge residuals are explained by @DavidKetcheson, the first order derivatives used by np.gradient at the edge. As I decrease the step size the residual becomes smaller. I could change my question to large residual at edges and accept an answer about the step sizes. – jensv Feb 5 '15 at 15:25