二阶 ODE dopri5 python UserWarning：更大的 nmax

Question

对于二阶 ODE（python 中的 dopri5 方法），下面的代码总是会导致错误：C:\Users\MY\Anaconda3\lib\site-packages\scipy\integrate\_ode.py:1019: UserWarning: dopri5: larger nmax is needed self.messages.get(idid, 'Unexpected idid=%s' % idid))。我已经更改了参数，但似乎没有任何帮助。即使设置 nsteps=100000 也不起作用。除了增加nsteps?

，还有其他方法可以解决这个问题吗？

from scipy.integrate import ode
import numpy as np

def fun(t, y):
    return np.array([y[1], -3/t*y[1] + 7/(t**6)*y[0]])

yinit = np.array([0.01, 0.2])

dt = 0.01
t_stop = 2

solver = ode(fun).set_integrator('dopri5', nsteps=100000).set_initial_value(yinit)
solver.t = 0.001
t_RK4_sci = [0]
x_RK4_sci = [yinit]
while solver.successful() and solver.t < t_stop:
    solver.integrate(solver.t+dt, step=True)
    t_RK4_sci.append(solver.t)
    x_RK4_sci.append(solver.y)
t_RK4_sci = np.array(t_RK4_sci)
x_RK4_sci = np.array(x_RK4_sci)

Answer 1

放一个

print t,y

作为 fun 中的第一行，可以看到您的解决方案在使用非常小的步长时迅速爆炸。该日志的最后几行是

0.00100025397168 [  2.57204893e+289   6.79981963e+298]
0.00100025397168 [  2.57204893e+289   6.79981963e+298]
0.00100025397168 [  2.57204893e+289   6.79981964e+298]
0.00100025397168 [  2.57204893e+289   6.79981964e+298]
0.00100025397168 [  2.57204897e+289   6.79981974e+298]
0.00100025397168 [  2.57204899e+289               inf]
0.00100025397168 [ inf  nan]
0.00100025397168 [ nan  nan]
0.00100025397168 [ nan  nan]
0.00100025397168 [ nan  nan]
0.00100025397168 [  2.57204894e+289   6.79981966e+298]
0.00100025397168 [  2.57204894e+289               inf]
0.00100025397168 [ inf  nan]
0.00100025397168 [ nan  nan]
0.00100025397168 [ nan  nan]
0.00100025397168 [ nan  nan]
/usr/lib/python2.7/dist-packages/scipy/integrate
/_ode.py:1018: UserWarning: dopri5: step size becomes too small
  self.messages.get(idid, 'Unexpected idid=%s' % idid))

要查看它的数学方面，请观察初始 Lipschitz 常数位于 L=1e+18 的某处。

数值积分的有用步长必须遵守 L*dt < 10，可能是一个更小的上限，以保持在 A-stability 区域内用于显式方法。
局部到全局误差的放大倍数为(exp(L*T)-1)，其中T为积分区间的长度。这意味着只能乐观地预期 L*T < 50 会产生有意义的结果，这会给出 T<5e-17.

在这些理论限制下，dopri5 积分器在实践中证明非常稳健，因为它仅在 T=2.5e-7 时退出。

欧拉形式的微扰

t²·y'' + 3t·y' - 7/t0^4·y = 0

在

范围内给出初始增长

(t/t0) ^ 3e6

并且由于最大的双倍数在 10^300 左右，因此数值范围在

左右超出

t/t0 = 10 ^ 1e-4 = 1.00023028502  or t=0.00100023028502

这是最接近数值积分退出的位置，因此可能是真正的原因。（更好的界限给出 10^(308/2.6e6)=1.00027280498。）

总结

这个微分方程不仅有一个非常大的 Lipschitz 常数，因此是 ill-behaved 或刚性的，而且精确解，就欧拉方程的近似值而言，增长得如此之快以超过double 数值积分失效时的范围。也就是说，即使是像隐式积分器这样更好的方法也不会给出更好的结果。

二阶 ODE dopri5 python UserWarning：更大的 nmax

Second-order ODE dopri5 python UserWarning: larger nmax

python

scipy

user-warning

ode

总结