[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.
>
> Thanks in advance,
>
> 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
More information about the SciPy-User
mailing list