python - 二阶 ODE dopri5 python UserWarning : larger nmax

Question

对于二阶 ODE(python 中的 dopri5 方法)，下面的代码总是会导致错误:C:\Users\MY\Anaconda3\lib\site-packages\scipy\integrate\_ode.py: 1019:用户警告:dopri5:需要更大的 nmax
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 稳定性区域内。

• 从局部到全局误差的放大率为(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 常数，因此表现不佳或僵硬，而且就欧拉方程的近似值而言，精确解的增长速度如此之快，超过了 double。数值积分失效时的范围。也就是说，即使是像隐式积分器这样更好的方法也不会给出更好的结果。