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

integrator = ode(der)
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:]):
    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,
  • $\begingroup$ 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. $\endgroup$ Commented Feb 5, 2015 at 4:53
  • $\begingroup$ 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? $\endgroup$
    – jensv
    Commented Feb 5, 2015 at 5:29
  • $\begingroup$ Why are you calculating the residual as $f (\xi'/f)'+f'\xi'-g\xi$? $\endgroup$
    – Kirill
    Commented Feb 5, 2015 at 6:37
  • $\begingroup$ results[1] is $f \xi'$. I am calculating $f (\frac{\xi' f}{f})'+f' \xi' - g \xi$ $\endgroup$
    – jensv
    Commented Feb 5, 2015 at 7:06
  • 1
    $\begingroup$ 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. $\endgroup$
    – jensv
    Commented Feb 5, 2015 at 15:25


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.