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

Question

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.

Answer 1

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.