我正在尝试将二阶 ODE 与接近初始条件的奇点进行积分。当我将积分结果插入 ODE 时，为什么会得到很大的残差？

该方程是来自等离子体物理领域的纽康的欧拉-拉格朗日方程。

$$\frac{d}{dr}(f \frac{d\xi}{dr}) - g \xi = 0$$ 或作为一组一阶 ODE： $$y_0 = \xi$$ $$y_1 = f \xi'$$ $$y_0' = \frac{y_1}{f}$$ $$y_1' = y_0 g$$

$f$和$g$是磁场和压力梯度的复杂表达式。但是，对于简单的情况，我什至遇到麻烦，例如：

$$f=r$$ $$g=-1 + r + \frac{1}{r}$$

ODE 总是单数的$r=0$ 从 Fobenius 扩展解决方案接近$r=0$对于这种情况应该是$\xi \propto r^1 \approx 0$和$\xi' \propto r^0 \approx 1$.

在下面的代码中，我设置了问题并使用 lsoda 求解器与 scipy.integrate.ode 集成。我将结果插入二阶 ODE 形式：

$$f \xi'' + f' \xi' -g \xi$$ 我必须使用 numpy.gradient 来计算$\xi''$.

我得到边缘顺序的残差$10^2$. 为什么它们这么大？如果我使用常量$f$和$g$我按照步长的顺序得到残差。

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])

我得到的输出看起来像下面的数组。边缘有很大的残差。

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])