# [SciPy-User] Problem with ODEINT

Warren Weckesser warren.weckesser@enthought....
Thu Nov 24 10:07:11 CST 2011

```On Thu, Nov 24, 2011 at 1:56 AM, Lofgren, Eric <elofgren@email.unc.edu>wrote:

> 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.
>
>
> 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.

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20111124/b3a262d2/attachment.html
```