4
$\begingroup$

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])
| cite | improve this question | |
$\endgroup$
  • $\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$ – David Ketcheson Feb 5 '15 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 Feb 5 '15 at 5:29
  • $\begingroup$ Why are you calculating the residual as $f (\xi'/f)'+f'\xi'-g\xi$? $\endgroup$ – Kirill Feb 5 '15 at 6:37
  • $\begingroup$ results[1] is $f \xi'$. I am calculating $f (\frac{\xi' f}{f})'+f' \xi' - g \xi$ $\endgroup$ – jensv Feb 5 '15 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 Feb 5 '15 at 15:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.