Second-order ODE dopri5 python UserWarning: larger nmax

Question

For a second-order ODE (dopri5 method in python) the code below always results in an error: 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)). I've changed the parameters but nothing seems to help. Even setting nsteps=100000 doesn't work. Is there any other way to solve this instead of just increasingnsteps?

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

Put a

print t,y

as first line in fun to see that your solution is rapidly exploding while employing very small step sizes. The last lines of that log are

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

---

To see the mathematical side of it observe that the initial Lipschitz constant is somewhere at L=1e+18.

• Useful step sizes for numerical integration have to observe L*dt < 10, probably a smaller upper bound, to stay inside the region of A-stability for explicit methods.

• The magnification rate from local to global error is (exp(L*T)-1), where T is the length of the integration interval. This means that meaningful results can only be expected for, optimistically, L*T < 50, which gives T<5e-17.

Under those theoretical restrictions, the dopri5 integrator proves as quite robust in practice since it bails out only at T=2.5e-7.

---

Perturbation to the Euler form

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

gives initial growth in the range of

(t/t0) ^ 3e6

and since the largest double is around 10^300 the numerical range gets exceeded at around

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

which is the closest to where the numerical integration bails out and thus probably the true reason. (Better bounds give 10^(308/2.6e6)=1.00027280498.)

---

Summary

Not only has this differential equation an exceedingly large Lipschitz constant and is thus ill-behaved or stiff, also the exact solution, as far as the approximation by an Euler equation can tell, grows so fast to exceed the double range at the time where the numerical integration breaks down. That is, even a better method like an implicit integrator will not give better results.