ОДУ второго порядка dopri5 python UserWarning: больше nmax

Для ODE второго порядка (метод dopri5 в python) приведенный ниже код всегда приводит к ошибке: 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)

python scipy ode user-warning
person KeVal    schedule 04.01.2016    source источник

Ответы (1)


arrow_upward
5
arrow_downward

Положите

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))

Чтобы увидеть математическую сторону этого, обратите внимание, что начальная константа Липшица находится где-то на уровне 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.)

Резюме

Мало того, что это дифференциальное уравнение имеет чрезвычайно большую постоянную Липшица и, таким образом, плохо себя ведет или является жестким, точное решение, насколько может сказать аппроксимация уравнением Эйлера, растет так быстро, что превышает диапазон double в то время, когда численное интегрирование не работает. То есть даже такой лучший метод, как неявный интегратор, не даст лучших результатов.

person Lutz Lehmann    schedule 04.01.2016