[SciPy-User] SciPy ODEINT Problem

[Lofgren, Eric]

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?

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)