二阶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稳定性区域内。

• 从局部误差到整体误差的放大率是(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常数，因此表现不佳或僵硬，而且精确的解决方案（通过Euler方程的近似值可以看出）增长得如此之快，超过double了数值时的范围。整合失败。也就是说，即使是更好的方法（如隐式积分器）也不会给出更好的结果。