[SciPy-User] Problem with ODEINT
Lofgren, Eric
elofgren@email.unc....
Thu Nov 24 01:56:36 CST 2011
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
More information about the SciPy-User
mailing list