# scipy odeint: excess work done on this call and very sensitive to initial value

I am trying out odeint and received the error

'Excess work done on this call (perhaps wrong Dfun type).'.

The values returned are also super sensitive to small changes of the initial value. Some initial values don't run into this problem whereas others do.

rho = 0.1
beta = 0.01
N = 50
gamma = 0.20
theta = 0.01
delta = 0.1
alpha = 1000
S = 1000
W = 180

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

def model(L,t):

R = (L * delta * alpha * 20) / (2* theta*S*N*W)
dLdt= rho*L-beta*N+(L*delta*alpha*R*W*20)/(N*S)

return dLdt

L0 = 20

# time points
t = np.linspace(0, 500,  1000)

# solve ODE
L, infodict = odeint(model,L0,t,full_output=True)
infodict['message']


The output L goes from 20 to 100 and then spikes to 1.62e+011 before going to zero again.

Essentially, by combining all the constant factors, your ODE is $$\frac{dL}{dt} = -A + B L + C L^2$$ With the initial $$L(0)=L_0$$ large enough, the positive feed-back of the quadratic term drives the equation towards a dynamic blow-up in finite time, as you observed with the values getting very large. This is a property of the exact solution, so it is no surprise that an accurate numerical solution will replicate this.

Your equation will only give bounded solutions if the initial point is below the positive root of the quadratic polynomial on the right side.

You have

A = beta*N
B = rho
C = ((delta * alpha * 20)/(N*S))**2/(2*theta)
print([A,B,C])
print(np.roots([C,B,-A]))


giving coefficients [0.5, 0.1, 0.08] and roots of the quadratic [-3.20194102 1.95194102]. For values slightly larger than that you will still reach the end of the integration interval before the blow-up point, however the values will be getting larger.

• Thank you, I need to have some initial values in which L is positive. However, for most initial values I tried, L returns as negative in the equilibrium. For initial values that return a positive trajectory for L, the spiking issue occurs, then L becomes zero again. How can I have L follow a positive trajectory, and what aspect of this formulation causes it to spike and then return to zero again? Oct 17 '20 at 11:06
• Yes, that is correct. As the quadratic coefficient is positive, the positive root is unstable and the negative root is stable, so that all bounded solution go towards this negative equilibrium point. This can only be changed if the signs of the coefficients, esp. the quadratic one, are different. Oct 17 '20 at 11:11