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]
print(f(r)*np.gradient(results[1]/f(r))/dr +
np.gradient(f(r))/dr*results[1]/f(r) - g(r)*results[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])