计算科学 - 量化 scipy ode / ODEPACK 的积分误差 - 吾爱随笔录

我正在尝试使用包装在 scipy.integrate.ode() 中的 lsoda 求解器将二阶 ODE 与潜在的几个奇点集成。我想在解决方案上放置一个错误栏，或者至少估计错误的上限。

了解如何计算简单有限差分方法的截断误差。但是从scipy 文档来看， LSODA 似乎更复杂，求解器使用自适应步长并自动在刚性和非刚性方法之间切换。

存在设置绝对容差、相对容差和最大步长的选项，但我不明白它们与本文档中我的解决方案中的错误有何关系。

在少数情况下，解析解是可能的。例如，我在上一个问题中调试的案例。

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

\frac{d}{d r} (f \frac{d ξ}{d r}) - g ξ = 0

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

y_{0} = ξ

$y_0 = \xi$

y_{1} = f ξ^{'}

$y_1 = f \xi'$

y_{0}^{'} = \frac{y_{1}}{f}

$y_0' = \frac{y_1}{f}$

y_{1}^{'} = y_{0} g

$y_1' = y_0 g$

$f$ 和 $g$ 是磁场和压力梯度的复杂表达式。对于更复杂的情况，我使用 scipy.interpolate.InterpolateUnivariateSpline() 生成的样条曲线来描述磁场和压力梯度。 $f$ 和 $g$ 可以快速变化并且 $f$ 可以在几个位置为零，从而导致奇点。我使用 Frobenius 展开来找到接近奇点的解。

这是一个简单案例的简短示例代码 $f$ 和 $g$ , 只有一个奇点在 $r=0$ .

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