2

对于二阶 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)
4

1 回答 1

5

放一个

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了数值集成崩溃。也就是说，即使是像隐式积分器这样更好的方法也不会给出更好的结果。

于 2016-01-04T14:45:33.997 回答