[SciPy-User] Problem with ODEINT

[Lofgren, Eric]

I'm attempting to solve a fairly simple compartmental SIR model (http://en.wikipedia.org/wiki/SIR_Model) using SciPy and particularly odeint . I've solved things like this before using this code, and they've always worked, so I'm frankly quite puzzled. The code is as follows:

---

#Python implementation of continuous SIR model

#Import Necessary Modules
import numpy as np
import scipy.integrate as spi

#Solving the differential equation. Solves over t for initial conditions PopIn

def eq_system(startState,t):
'''Defining SIR System of Equations'''
beta = startState[3]
gamma = startState[4]
#Creating an array of equations
Eqs= np.zeros((3))
Eqs[0]= -beta * startState[0]*startState[1]
Eqs[1]= beta * startState[0]*startState[1] - gamma*startState[1]
Eqs[2]= gamma*startState[1]
return Eqs

def model_solve(t):
'''Stores all model parameters, runs ODE solver and tracks results'''
#Setting up how long the simulation will run
t_start = 1
t_end = t
t_step = 0.02
t_interval = np.arange(t_start,t_end, t_step)
n_steps = (t_end-t_start)/t_step
#Setting up initial population state
#n_params is the number of parameters (beta and gamma in this case)
S0 = 0.99999
I0 = 0.00001
R0 = 0.0
beta = 0.50
gamma = 1/10.
n_params = 2
startPop = (S0, I0, R0, beta, gamma)
#Create an array the size of the ODE solver output with the parameter values
params = np.zeros((n_steps,n_params))
params[:,0] = beta
params[:,1] = gamma
timer = np.arange(n_steps).reshape(n_steps,1)
SIR = spi.odeint(eq_system, startPop, t_interval)
#Glue together ODE model output and parameter values in one big array
output = hstack((timer,SIR,params))
return output

def MonteCarlo_SIR(runs):
holder = [ ]
for i in range(runs):
holder.append(model_solve(100))
results = np.hstack(holder)
return results

testing = MonteCarlo_SIR(10)

print testing
print testing.shape

---

Ignore for a moment that there's absolutely no reason for a Monte Carlo simulation of a system made up of nothing but constants. And the poorly documented code - this is essentially a running musing right now, and wasn't intended to see the light of day until I hit this problem.

When I run this code, sometimes it just works, and I get a nice tidy array of results. But sometimes, errors like the following crop up:

lsoda-- warning..internal t (=r1) and h (=r2) are
such that in the machine, t + h = t on the next step
(h = step size). solver will continue anyway
In above, R1 = 0.1000000000000E+01 R2 = 0.0000000000000E+00
intdy-- t (=r1) illegal
In above message, R1 = 0.1020000000000E+01
t not in interval tcur - hu (= r1) to tcur (=r2)
In above, R1 = 0.1000000000000E+01 R2 = 0.1000000000000E+01
intdy-- t (=r1) illegal
In above message, R1 = 0.1040000000000E+01
t not in interval tcur - hu (= r1) to tcur (=r2)
In above, R1 = 0.1000000000000E+01 R2 = 0.1000000000000E+01
lsoda-- trouble from intdy. itask = i1, tout = r1
In above message, I1 = 1
In above message, R1 = 0.1040000000000E+01
Illegal input detected (internal error).
Run with full_output = 1 to get quantitative information.

lsoda-- at t (=r1) and step size h (=r2), the
corrector convergence failed repeatedly
or with abs(h) = hmin
In above, R1 = 0.6556352055692E+01 R2 = 0.2169078563087E-05
Repeated convergence failures (perhaps bad Jacobian or tolerances).
Run with full_output = 1 to get quantitative information.

lsoda-- warning..internal t (=r1) and h (=r2) are
such that in the machine, t + h = t on the next step
(h = step size). solver will continue anyway
In above, R1 = 0.1000000000000E+01 R2 = 0.5798211925922E-81
lsoda-- warning..internal t (=r1) and h (=r2) are
such that in the machine, t + h = t on the next step
(h = step size). solver will continue anyway
In above, R1 = 0.1000000000000E+01 R2 = 0.8614947121323E-85
...
lsoda-- warning..internal t (=r1) and h (=r2) are
such that in the machine, t + h = t on the next step
(h = step size). solver will continue anyway
In above, R1 = 0.1000000000000E+01 R2 = 0.1722989424265E-83
lsoda-- above warning has been issued i1 times.
it will not be issued again for this problem
In above message, I1 = 10

I believe those three are are a full tour of the error messages that are cropping up. This definitely occurs a minority of the time, but it's common enough that 4 out of 10 runs of the code above produces at least one error message like that. It seems that the problem is the step size getting small enough that it's beyond the precision of the machine to deal with, but passing an argument of something like hmin = 0.0000001 or something that should be well within range doesn't seem to help.

Any idea what's going on? The end goal of this is a somewhat more complex set of equations, and the parameters (i.e. beta and gamma) to be drawn from a distribution, but I'm somewhat concerned about the trouble the solver is seeming to have with what should be a pretty straightforward case. Any suggestions on how to fix this? Or is there another type of ODE solver I should be using? I confess my knowledge of how these things work is...limited.

Thanks in advance,

Eric

_______________________________________________
SciPy-User mailing list
SciPy...@scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user

[Flavio Coelho]

Hi,

I have seen this before with odeint. It is probably related to the parameterization of your model, which is leading very large derivatives. This is causing the numerical solver to try to reduce the step size (h) to try to resolve the dynamics. h is apparently hitting its default lower limit (something in the vicinity of 1E-85).  

You can try to set a lower hmin, when calling odeint, but it would be best to go over your parameter values and check if they are reasonable.

Disclaimer: Ihave not run your script nor checked your parameters, just took a lok at the error messages.

good luck,

Flávio

--
Flávio Codeço Coelho
================
+55(21) 3799-5567
Professor
Escola de Matemática Aplicada 
Fundação Getúlio Vargas
Rio de Janeiro - RJ
Brasil

[Warren Weckesser]

Eric,

You have given odeint an initial condition of length 5, but the function
that defines your system is returning a vector of only length 3.  Don't do
that.

Instead of including the parameters beta and gamma in the initial
condition, you can make them explicit arguments.  For example:
 
 def eq_system(state, t, beta, gamma):

   '''Defining SIR System of Equations'''

   Eqs = np.zeros((3))
   Eqs[0] = -beta * state[0]*state[1]
   Eqs[1] = beta * state[0]*state[1] - gamma*state[1]
   Eqs[2] = gamma*state[1]
   return Eqs


Then call odeint like this:

   startPop = (S0, I0, R0)
   SIR = spi.odeint(eq_system, startPop, t_interval, args=(beta, gamma))


When I make these changes to your script, I don't get any errors, and
it runs much faster.

There are additional examples of using parameters with odeint in the
SciPy Cookbook:
    http://www.scipy.org/Cookbook/CoupledSpringMassSystem
    http://www.scipy.org/Cookbook/KdV
   
   
If you really must include the parameters in the state vector, then you
must also define differential equations for them in your system
(i.e. d(beta)/dt = 0, etc):

def eq_system(state,t):

   '''Defining SIR System of Equations'''

   beta = state[3]
   gamma = state[4]

   #Creating an array of equations

   Eqs = np.zeros(5)
   Eqs[0] = -beta * state[0]*state[1]
   Eqs[1] = beta * state[0]*state[1] - gamma*state[1]
   Eqs[2] = gamma*state[1]
   Eqs[3] = 0
   Eqs[4] = 0
   return Eqs


Warren

[Lofgren, Eric]

I checked the parameters, and they seem to both be reasonable values for the system in question. What's especially confusing to be is if I change the code to something like:

SIR = spi.odeint(eq_system, startPop, t_interval, hmin=1e-20)

where I've set a smaller step-size, I still get the occasional error like this:

lsoda-- warning..internal t (=r1) and h (=r2) are
such that in the machine, t + h = t on the next step
(h = step size). solver will continue anyway

In above, R1 = 0.1000000000000E+01 R2 = 0.5798216400482E-81
lsoda-- at t(=r1) and step size h(=r2), the error
test failed repeatedly or with abs(h) = hmin
In above, R1 = 0.1000000000000E+01 R2 = 0.5798216400482E-81
Repeated error test failures (internal error).

Run with full_output = 1 to get quantitative information.

Which suggests the solver is blazing right past the minimum step size I just set.

Eric

On Nov 24, 2011, at 10:56 AM, <scipy-use...@scipy.org>
<scipy-use...@scipy.org> wrote:

> Date: Thu, 24 Nov 2011 13:33:57 -0200
> From: Flavio Coelho <fcco...@gmail.com>
> Subject: Re: [SciPy-User] Problem with ODEINT
> To: SciPy Users List <scipy...@scipy.org>
> Message-ID:
> <CAAbyzbQP0JV3fyKsH3XiisZAhr-tBTik1Z8HEwi=_VJWg...@mail.gmail.com>
> Content-Type: text/plain; charset="iso-8859-1"

>
> Hi,
>
> I have seen this before with odeint. It is probably related to the
> parameterization of your model, which is leading very large derivatives.
> This is causing the numerical solver to try to reduce the step size (h) to
> try to resolve the dynamics. h is apparently hitting its default lower
> limit (something in the vicinity of 1E-85).
>
> You can try to set a lower hmin, when calling odeint, but it would be best
> to go over your parameter values and check if they are reasonable.
>
> Disclaimer: Ihave not run your script nor checked your parameters, just
> took a lok at the error messages.
>
> good luck,
>

> Fl?vio

[Lofgren, Eric]

> Eric,
> You have given odeint an initial condition of length 5, but the function
> that defines your system is returning a vector of only length 3. Don't do
> that.

> ...

> Warren

Warren-

This does indeed seem to solve the problem, I haven't hit any errors in 1,000 or so runs, and it does indeed make the code run considerably faster. Thank you for the advice and help.

Eric