#bloc1.py
#integracio de les equacions del model SIR
from numpy import *
import pylab as p
from scipy import integrate
# Definim els parametres
# Infeccio
r=0.01
# Recuperacio
b=0.4
# Condicions inicials
S0=990
I0=10
R0=0
def SIR(X, t=0):
""" Retorna les equacions. El t=0 es perque el sistema es autonom """
return array([-r*X[0]*X[1],r*X[0]*X[1] -b*X[1],b*X[1] ])
# Passem a la integracio
t = linspace(0, 10, 1000) # temps integracio
X0 = array([S0,I0,R0])
X = integrate.odeint(SIR, X0, t)
# Passem a representar-ho:
S, I, R = X.T
#Aixo seria el mateix que si fessim:
#presa=X[:,0]
#predador=X[:,1]
f1 = p.figure() #iniciem una figura
p.plot(t, S, 'g-', label='Susceptibles x (x(0)=%g)' %(S[0]))
p.plot(t, I, 'r-', label='Infectats y (y(0)=%g)' %(I[0]))
p.plot(t, R, 'y-', label='Recuperats y (y(0)=%g)' %(R[0]))
p.grid() # Ens mostra una graella
p.legend(loc='best') # Ens mostra una llegenda al millor lloc
p.xlabel('$t$')
p.ylabel('poblacio')
p.title("r=%g, b=%g, $R_0$=%g" % (r,b,r*X0[0]/b))
p.show()6 Models Epidèmics
6.1 El model SIR
El model epidemiològic SIR \[\begin{array}{rcl} S' & = & -rSI \\ I' & = & rSI - b I \\ R' & = & bI \end{array}\] el podem simular integrant les equacions diferencials.
Si volem veure les solucions en 3d, podem usar el paquet mlab del mayavi. Per això cal haver obert el Ipython amb la comanda -wthread
# bloc2.py
# ho fem en 3d
from enthought.mayavi.tools import mlab
from numpy import *
from scipy import integrate
# Definim els parametres
# Infeccio
r=0.01
# Recuperacio
b=0.4
# Condicions inicials
S0=990
I0=10
R0=0
def SIR(X, t=0):
""" Retorna les equacions. El t=0 es perque el sistema es autonom """
return array([-r*X[0]*X[1],r*X[0]*X[1] -b*X[1],b*X[1] ])
# Passem a la integracio
t = linspace(0, 10, 1000) # temps integracio
X0 = array([S0,I0,R0])
X = integrate.odeint(SIR, X0, t)
# Passem a representar-ho:
S, I, R = X.T
mlab.plot3d(S,I,R)Si volem veure un model més espectacular podem fer-ho amb l’atractor de Rössler:
# bloc3.py
# ho fem en 3d
# En un model mes interessant
from enthought.mayavi.tools import mlab
from numpy import *
from scipy import integrate
# Definim els parametres
# Condicions inicials
S0=0
I0=0
R0=0
a=0.432
b=2.0
c=4.0
def Rossler(X, t=0):
""" Retorna les equacions. El t=0 es perque el sistema es autonom """
return array([-X[1]-X[2],X[0]+a*X[1],b+X[2]*(X[0]-c)])
# Passem a la integracio
t = linspace(0, 1000, 10000) # temps integracio
X0 = array([S0,I0,R0])
X = integrate.odeint(Rossler, X0, t)
# Passem a representar-ho:
S, I, R = X.T
mlab.plot3d(S,I,R)6.2 Exercici: Ajustament de la plaga de Bombay
El fitxer bombay.txt inclou les dades de mortalitat setmanal de l’epidèmia de pesta de bombay del 1905-1906 que varen usar Kermack-McKendrick (1927) per formular el seu model SIR. En primer lloc ho representem:
# bloc3.py
# ho fem en 3d
# En un model mes interessant
from enthought.mayavi.tools import mlab
from numpy import *
from scipy import integrate
# Definim els parametres
# Condicions inicials
S0=0
I0=0
R0=0
a=0.432
b=2.0
c=4.0
def Rossler(X, t=0):
""" Retorna les equacions. El t=0 es perque el sistema es autonom """
return array([-X[1]-X[2],X[0]+a*X[1],b+X[2]*(X[0]-c)])
# Passem a la integracio
t = linspace(0, 1000, 10000) # temps integracio
X0 = array([S0,I0,R0])
X = integrate.odeint(Rossler, X0, t)
# Passem a representar-ho:
S, I, R = X.T
mlab.plot3d(S,I,R)Segons s’explica en el llibre de Murray(pàg. 324), Kermack i McKendrick varen aproximar la funció \[\frac{dR}{dt}= a \mathrm{sech }(b t + c)\] on els paràmetres \(a\), \(b\) i \(c\) tenen relació amb els paràmetres originals. Trobeu l’ajustament que varen aconseguir. Podeu provar com a valors inicials \[a \approx 890, \quad b= 0.2 \quad c= -3.4\] Usant la referència del llibre de Murray, identifiqueu el valor de \(R_0\).