s3Eva_2020PAOI_T2 Modelo epidemiológico no letal

El ejercicio representa un sistema de ecuaciones diferenciales ordinarias, que serán resueltas usando Runge-Kutta de 2do Orden.

De compararse con la curva de contagios de Covid-19 se tienen diferencias en la población recuperada, pues el modelo se considera no letal por lo que no se contabiliza el número de fallecidos.

El módelo es el más básico y permite cambiar por ejemplo la tasa de infección, y se ve los cambios en la curva de infectados. Se puede observar lo que se indicaba como objetivo de «aplanar la curva» al disminuir la población expuesta mediante cambiar la tasa de infección al exponer más o menos población al contagio por iteacción entre «suceptibles» e «infectados.


Desarrollo analítico

Las fórmulas para el algoritmo se identifican como:

binfecta = 1.4
grecupera = 1/4
# Ecuaciones
fS = lambda t,S,I,R : -binfecta*S*I
gI = lambda t,S,I,R : binfecta*S*I - grecupera*I
qR = lambda t,S,I,R : grecupera*I

que luego se usan en cada iteración que se registra en la tabla empezando con las condiciones iniciales

itera Si Ii Ri
0 1 0.001 0
1 0.997797107 0.002809143 0.00039375
2 0.9916392093856501 0.007862023500353846 0.0014987671139961277
3 Tarea

itera= 1

K1S = h * fS(ti,Si,Ii,Ri) 
    = 1(-1.4*1*0.001) = -0.0014
K1I = h * gI(ti,Si,Ii,Ri) 
    = 1(1.4*1*0.001 - (1/4)*0.001) = 0.00115
K1R = h * qR(ti,Si,Ii,Ri) 
    = 1((1/4)*0.001)  = 0.00025

K2x = h * fS(ti+h, Si + K1S, Ii+K1I, Ri +K1R) 
    = 1(-1.4*(1-0.0014)(0.001+0.00115)) = -0.003005786
K2y = h * gI(ti+h, Si + K1S, Ii+K1I, Ri +K1R) 
    = 1(1.4*(1-0.0014)*(0.001+0.00115)
        -(1/4)*(0.001+0.00115)) 
    = 0.002468286
K2z = h * qR(ti+h, Si + K1S, Ii+K1I, Ri +K1R) 
    = 1( (1/4)*(0.001+0.00115)) = 0.0005375

Si = Si + (1/2)*(K1S+K2S) 
   = 1 + (-0.0014 -0.003005786)/2 = 0.997797107
Ii = Ii + (1/2)*(K1I+K2I) 
   = 0.001 + (0.00115+0.002468286)/2 = 0.002809143
Ri = Ri + (1/2)*(K1R+K2R) 
   = 0 + (0.00025 + 0.0005375)/2 = 0.00039375
ti = ti + h = 1 + 1 = 2

itera = 2

K1S = 1(-1.4*0.997797107* 0.002809143)
    = -0.003924136661969021
K1I = 1(1.4*0.997797107* 0.002809143 
        - (1/4)*0.002809143)
    = 0.003221850911969021
K1R = 1((1/4)*0.002809143)
    = 0.00070228575

K2S = 1(-1.4*(0.997797107-0.003924136661969021)*
        *(0.002809143+0.003221850911969021))
    = -0.008391658566730926
K2I = 1(1.4*(0.997797107-0.003924136661969021)*
        *(0.002809143+0.003221850911969021) 
- (1/4)*(0.002809143+0.003221850911969021))
    = 0.006883910088738671
K2R = 1( (1/4)*(0.002809143+0.003221850911969021))
    = 0.0015077484779922553

Si = 0.997797107 + 
     + (-0.003924136661969021 -0.008391658566730926)/2
   = 0.9916392093856501
Ii =  0.002809143 + 
      + (0.003221850911969021+0.006883910088738671)/2
   = 0.007862023500353846
Ri = 0.00039375 + 
     + (0.00070228575 + 0.0015077484779922553)/2
   = 0.0014987671139961277
ti = ti + h = 2 + 1 = 3

itera 3 – TAREA


Instrucciones en Python

# 3Eva_2020PAOI_T2 Modelo epidemiológico no letal
# Sistemas EDO con Runge Kutta de 2do Orden
import numpy as np

def rungekutta2_fgq(f,g,q,t0,x0,y0,z0,h,muestras):
    tamano = muestras +1
    tabla = np.zeros(shape=(tamano,4),dtype=float)
    tabla[0] = [t0,x0,y0,z0]
    ti = t0
    xi = x0
    yi = y0
    zi = z0
    for i in range(1,tamano,1):
        K1x = h * f(ti,xi,yi,zi)
        K1y = h * g(ti,xi,yi,zi)
        K1z = h * q(ti,xi,yi,zi)
        
        K2x = h * f(ti+h, xi + K1x, yi+K1y, zi +K1z)
        K2y = h * g(ti+h, xi + K1x, yi+K1y, zi +K1z)
        K2z = h * q(ti+h, xi + K1x, yi+K1y, zi +K1z)

        xi = xi + (1/2)*(K1x+K2x)
        yi = yi + (1/2)*(K1y+K2y)
        zi = zi + (1/2)*(K1z+K2z)
        ti = ti + h
        
        tabla[i] = [ti,xi,yi,zi]
    tabla = np.array(tabla)
    return(tabla)

# Programa
# Parámetros de las ecuaciones

binfecta = 1.4
grecupera = 1/4
# Ecuaciones
fS = lambda t,S,I,R : -binfecta*S*I
gI = lambda t,S,I,R : binfecta*S*I - grecupera*I
qR = lambda t,S,I,R : grecupera*I
# Condiciones iniciales
t0 = 0
S0 = 1.0
I0 = 0.001
R0 = 0.0

# parámetros del algoritmo
h = 1.0
muestras = 51

# PROCEDIMIENTO
tabla = rungekutta2_fgq(fS,gI,qR,t0,S0,I0,R0,h,muestras)
ti = tabla[:,0]
Si = tabla[:,1]
Ii = tabla[:,2]
Ri = tabla[:,3]
# SALIDA
np.set_printoptions(precision=6)
print(' [ ti, Si, Ii, Ri]')
print(tabla)

# Grafica tiempos vs población
import matplotlib.pyplot as plt
plt.plot(ti,Si, label='S')
plt.plot(ti,Ii, label='I')
plt.plot(ti,Ri, label='R')
plt.title('Modelo SIR')
plt.xlabel('t tiempo')
plt.ylabel('población')
plt.legend()
plt.grid()
plt.show()