[SciPy-User] Problem with ODEINT
Flavio Coelho
fccoelho at gmail.com
Thu Nov 24 10:33:57 EST 2011
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
On Thu, Nov 24, 2011 at 05:56, Lofgren, Eric <elofgren at 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
>
>
>
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
--
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20111124/72b0f17d/attachment.html>
More information about the SciPy-User
mailing list