scipy odeint result depends on input time array

Question

I'm trying to solve the following one-dimensional ODE (2nd order) using the scipy odeint function:

z'' = 2*F*cos(vt-kz)*(z*cos(vt-k*z) - k*(1+z^2)*sin(vt-kz))/(1+z^2)^2

And my python script is

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def func(var,t,F,k,v):
    z,Vz = var
    tmp = v*t-k*z
    cos_tmp = np.cos(tmp)
    denom = 1+z**2

    return [Vz,
            2*F*cos_tmp*(z*cos_tmp - k*denom*np.sin(tmp))/(denom*denom)]

I understand that odeint uses an adaptive time step and it returns the solution evaluated at the time points given by the input t array. However, I observed that the solution diverges when the time array is not fine enough which doesn't make sense to me if odeint chooses the internal time step automatically.

k= 12184.
F = -0.000125597154176
v = 0.3
y0 = [-1.0,0.]

t1 = np.linspace(0,10,10000)
t2 = np.linspace(0,10,100)
t3 = np.linspace(0,10,10)

result1 = odeint(func,y0,t1,args=(F,k,v))
result2 = odeint(func,y0,t2,args=(F,k,v))
result3 = odeint(func,y0,t3,args=(F,k,v))

fig,ax = plt.subplots(1,2,figsize=(14,5))

ax[0].plot(t1,result1[:,0],label='t1')
ax[0].plot(t2,result2[:,0],'r.',label='t2')
ax[1].plot(t3,result3[:,0],'r.',label='t3')

ax[0].legend(frameon=False,loc='upper left',numpoints=1)
ax[1].legend(frameon=False,loc='upper left',numpoints=1)

plt.show()

An image of the solution is here:

Warren Weckesser · Accepted Answer · 2017-07-11T12:25:46.487

When you ran the code, did you notice the warning that is printed? It says

 lsoda--  at current t (=r1), mxstep (=i1) steps   
       taken on this call before reaching tout     
      in above message,  i1 =       500
      in above message,  r1 =  0.3622415764602D+00
[...]/scipy/integrate/odepack.py:218: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
  warnings.warn(warning_msg, ODEintWarning)

It is a bit cryptic, but what it means is that odeint hit the limit of 500 internal time steps before reaching a requested time.

You can increase that maximum by setting the mxstep argument. I just tried your code using

result3 = odeint(func,y0,t3,args=(F,k,v), mxstep=2000)

and it worked.

scipy odeint result depends on input time array

1 Answers1