I'm implementing a very simple Susceptible-Infected-Recovered model with a steady population for an idle side project - normally a pretty trivial task. But I'm running into solver errors using either PysCeS or SciPy, both of which use lsoda as their underlying solver. This only happens for particular values of a parameter, and I'm stumped as to why. The code I'm using is as follows:

import numpy as np
from pylab import *
import scipy.integrate as spi

#Parameter Values
S0 = 99.
I0 = 1.
R0 = 0.
N0 = S0 + I0 + R0
PopIn= (S0, I0, R0, N0)
beta= 0.50     
mu = 1/25550.
t_end = 15000.
t_start = 1.
t_step = 1.
t_interval = np.arange(t_start, t_end, t_step)

#Solving the differential equation. Solves over t for initial conditions PopIn
def eq_system(PopIn,t):
    '''Defining SIR System of Equations'''
    #Creating an array of equations
    Eqs= np.zeros((3))
    Eqs[0]= -beta * (PopIn[0]*PopIn[1]/PopIn[3]) - mu*PopIn[0] + mu*PopIn[3]
    Eqs[1]= (beta * (PopIn[0]*PopIn[1]/PopIn[3]) - gamma*PopIn[1] - mu*PopIn[1])
    Eqs[2]= gamma*PopIn[1] - mu*PopIn[2]
    return Eqs

SIR = spi.odeint(eq_system, PopIn, t_interval)

This produces an error:

 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.7732042715460E+04
Excess work done on this call (perhaps wrong Dfun type).
Run with full_output = 1 to get quantitative information.

However, changing "mu" (which is currently 1/70 years, a pretty common value for life expectancy) to either 26550 or 22550 produces no such error, though some values between 25550 and 22550 seem to.

I'm a little puzzled as to why that combination in particular is producing errors. Anyone have some insight?


2 Answers 2


It seems that your system blows up slightly above t=7732. You can check this by setting t_end=7732 and plotting the result using -

plot( SIR[:,0] ); show()

Your differential equation might have a problem.

  • $\begingroup$ There's clearly something wrong, but having stared at it, I'm not seeing it. $\endgroup$
    – Fomite
    Apr 4, 2013 at 21:34
  • $\begingroup$ The initial condition is a vector of length 4 whereas the function returns a vector of length 3. Could that be the problem? $\endgroup$
    – Rajeev
    Apr 5, 2013 at 4:29

I am guessing that the derivative of one of the functions is already machine precision 0 as the error also indicates. Your entire ODE system evolves within tend < 150 and is in steady state during the next following 14850 steps. At some point the derivative between two consecutive time steps is evaluate to 0.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.