[SciPy-User] SciPy ODEINT Problem

Ralf Gommers ralf.gommers@gmail....
Mon Oct 8 13:54:54 CDT 2012


On Mon, Oct 8, 2012 at 8:33 AM, Lofgren, Eric <elofgren@email.unc.edu>wrote:

> I've been working on a set of ordinary differential equations for an
> epidemic model. Normally odeint works swimmingly for this kind of thing,
> though I'll admit this one is a touch more complex than most of the ones
> I've thrown at it. I'm getting the following error message when I try to
> run the code below:
>
>  lsoda--  at t (=r1) and step size h (=r2), the
>        corrector convergence failed repeatedly
>        or with abs(h) = hmin
>       In above,  R1 =  0.7016749763132E+04   R2 =  0.1954514552051E-06
> Repeated convergence failures (perhaps bad Jacobian or tolerances).
>
> Having used the atol and rtol options to relax the tolerances by quite a
> bit, the error then converts to:
>
> RuntimeWarning: overflow encountered in double_scalars
>
> This feels to me like a coding error at that point, rather than some issue
> with the solver itself, but I've not managed to find anything in
> particular. Any ideas?
>

Have you looked at the returned solution? If you plot it:

  import matplotlib.pyplot as plt
  plt.plot(model)
  plt.ylim([-10, 20])
  plt.show()

you'll see that a few of the curves are stable and a few other ones run off
to infinity quickly. Reducing the time step doesn't change that. So your
model may have an error in it.

Ralf



> Thanks,
>
> Eric
>
> # Imports
> import numpy as np
> from pylab import *
> import scipy.integrate as spi
>
> # Initial Population States - Model is in Individuals
> Us0 = 2.0
> H0 = 1.0
> Up0 = 4.0
> Cp0 = 4.0
> Ua0 = 4.0
> Ca0 = 4.0
> D0 = 0.0
> PopIn = (Us0, H0, Up0, Cp0, Ua0, Ca0, D0)
>
> # Model parameters
> # Time is currently in MINUTES
> N = Us0 + H0 + Up0 + Cp0 + Ua0 + Ca0 + D0
> M = Up0 + Cp0 + Ua0 + Ca0 + D0
> n_contacts = (3.0)*(1/20.0)
> p_contacts = nurse_contacts/N
> rho_p = p_contacts
> rho_d = p_contacts
> rho_a = p_contacts
> sigma_p = 0.05
> sigma_d = 0.25
> sigma_a = 0.05
> omega = 14400
> alpha = 0.25
> psi_p = 0.10
> psi_a =  0.10
> mu_p = 0.0
> mu_a = 0.0
> theta_p = 1.0/10080.0
> theta_a = 1.0/10080.0
> nu_cp = 0.07
> nu_ca = 0.07
> nu_d = 0.0
> nu_up = 0.43
> nu_ua = 0.43
> kappa = 0.10
> tau = 3.50
> iota =  1/20.0 * 0.60 * 1.00
>
> zeta =  0.2785 * (1.0/17280.0) # Pr(death) and time until death
> gamma = (1.0-0.2785) * (1.0/12902.4) # Pr(discharge) and time until
> discharge
> theta = theta_p + theta_a
>
>
> # ODE Fit and Graphing Parameters
> t_end = 144000
> t_start = 1
> t_step = 0.1
> t_interval = np.arange(t_start, t_end, t_step)
>
> # The actual model running part
>
> def eq_system(PopIn,t):
>     #Creating an array of equations
>     Eqs= np.zeros((7))
>     Eqs[0] = ((iota*PopIn[1]) - (rho_p*sigma_p*PopIn[3]*(PopIn[0]/N)) -
>         (rho_d*sigma_d*PopIn[6]*(PopIn[0]/N)) -
> (rho_a*sigma_a*PopIn[5]*(PopIn[0]/N)))
>     Eqs[1] = ((rho_p*sigma_p*PopIn[3]*(PopIn[0]/N)) +
> (rho_d*sigma_d*PopIn[6]*(PopIn[0]/N))
>         + (rho_a*sigma_a*PopIn[5]*(PopIn[0]/N)) - (iota*PopIn[1]))
>     Eqs[2] = (((1/omega)*PopIn[4]) - alpha*PopIn[2] -
> (rho_p*psi_p*PopIn[2]*(PopIn[1]/N))
>         - (mu_p*sigma_p*PopIn[2]*(PopIn[3]/N)) -
> (mu_a*sigma_a*PopIn[2]*(PopIn[5]/N))
>         - theta_p*PopIn[2] +
> nu_up*((theta*M)+(zeta*PopIn[6])+(gamma*PopIn[6])))
>     Eqs[3] = (alpha*PopIn[2] - ((1/omega)*PopIn[4]) -
> (rho_a*psi_a*PopIn[4]*(PopIn[1]/N))
>         - (mu_p*sigma_p*PopIn[4]*(PopIn[3]/N)) -
> (mu_a*sigma_a*PopIn[4]*(PopIn[5]/N))
>         - theta_a*PopIn[4] +
> nu_ua*((theta*M)+(zeta*PopIn[6])+(gamma*PopIn[6])))
>     Eqs[4] = (((1/omega)*PopIn[5]) + (rho_p*psi_p*PopIn[2]*(PopIn[1]/N))
>         + (mu_p*sigma_p*PopIn[4]*(PopIn[3]/N)) +
> (mu_a*sigma_a*PopIn[2]*(PopIn[5]/N))
>         - alpha*PopIn[3] - kappa*PopIn[3] - theta_p*PopIn[3]
>         + nu_cp*((theta*M)+(zeta*PopIn[6])+(gamma*PopIn[6])))
>     Eqs[5] = (alpha*PopIn[3] + (rho_a*psi_a*PopIn[4]*(PopIn[1]/N))
>         + (mu_p*sigma_p*PopIn[4]*(PopIn[3]/N)) +
> (mu_a*sigma_a*PopIn[4]*(PopIn[5]/N)) - ((1/omega)*PopIn[5])
>         - kappa*tau*PopIn[5] - theta_a*PopIn[5] +
> nu_ca*((theta*M)+(zeta*PopIn[6])+(gamma*PopIn[6])))
>     Eqs[6] = (kappa*PopIn[3] + kappa*tau*PopIn[5] - gamma*PopIn[6]
>               - zeta*PopIn[6] +
> nu_d*((theta*M)+(zeta*PopIn[6])+(gamma*PopIn[6])))
>     return Eqs
>
> # Model Solver
> model = spi.odeint(eq_system, PopIn, t_interval)
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20121008/04126682/attachment.html 


More information about the SciPy-User mailing list