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 gamma=1/10. 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= -beta * (PopIn*PopIn/PopIn) - mu*PopIn + mu*PopIn Eqs= (beta * (PopIn*PopIn/PopIn) - gamma*PopIn - mu*PopIn) Eqs= gamma*PopIn - mu*PopIn 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?