Modelo para epidemias (Grupo 6-A)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Modelo para epidemias. Grupo 6-A |
| Asignatura | Ecuaciones Diferenciales |
| Curso | Curso 2014-15 |
| Autores | Álvaro Baeza Cabrero, Daniel Fojo Berlana, Gonzalo Bolea Muguruza, Pablo Carrasco del Olmo, Jorge Juan Fernández Díaz |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
1 Introducción
En este artículo se va a realizar el estudio del comportamiento temporal de un enfermedad infecciosa. Usaremos las siguientes hipótesis:
- La población será un número fijo N en el que todos serán susceptibles a la enfermedad.
- Los individuos infectados no se curarán durante el periodo de la enfermedad y serán contagiosos.
- La unidad de tiempo será la semana.
- Durante cada unidad de tiempo cada persona infectada tiene c contactos.
2 Determinación del número de contactos c
Primero en nuestro estudio vamos a determinar el número de contactos que tiene una persona infectada por unidad de tiempo de la siguiente ecuación diferencial que no proporciona el enunciado del trabajo:
- [math] I'(t) = \frac{c}{N}I(t)(N-I(t)) [/math]
c sería un número de dos decimales comprendido entre [math] 0.01 [/math] y [math] 0.99 [/math]. Supongamos que N sea igual a 500000 y que tiene las siguientes condiciones iniciales. [math] I(0)=200 [/math] y [math] I(1)=500 [/math] Para determinar c vamos a imponer una condición, que [math] |I(1)-500| [/math] tenga el valor más bajo posible, es decir, aquel valor que minimice lo posible el error absoluto. Para hallar ese valor lo resolveremos numéricamente con el siguiente programa escrito en MATLAB:
y0 = 200; % Valor inicial
y = [];
c = 1/100:0.01:99/100; % Posibles valores de c
N = 500000; % Población total
k = length(c);
C_min = c(1);
Y = N;
for i = 1:k
Y_final = Y;
Y = N*y0/(y0+(N-y0).*exp(-c(i))); % Modelo para t=1
y = [y Y];
if abs(Y-500) <= abs(Y_final-500)
C_min = c(i);
end
end
min(abs(y-500)) % Error absoluto mínimo que se comete
C_min
plot(c,abs(y-500))Ejecutando este programa nos da el valor c que minimiza el error que será igual a [math] 0.92 [/math] También hemos dibujado un gráfico para ver como evolucionaría el error para cada valor de c:
3 Estudio de la evolución de la enfermedad
Una vez hallado el valor c vamos a estudiar cómo la enfermedad va evolucionando a lo largo de las semanas, para ello utilizaremos el método de Heun y Runge-Kutta.
3.1 Método de Heun
Es una aproximación de la solución exacta de una ecuación diferencial en el que se descompone la solución en dos soluciones más simples y se halla una función media de las dos anteriores, se calcula los siguientes elementos:
- [math] K_{1}=f(t_{n},y_{n}) [/math]
- [math] K_{2}=f(t_{n}+h,y_{n}+K_{1}h) [/math]
- [math] y_{n+1}=y_{n}+\frac{h}{2}(K_{1}+K_{2}) [/math]
Siendo h la longitud de paso entre dos valores consecutivos
3.2 Método de Runge-Kutta
Sería similar al método de Heun pero se descompone en cuatro soluciones, se hace del siguiente modo:
- [math] K_{1}=f(t_{n},y_{n}) [/math]
- [math] K_{2}=f(t_{n}+\frac{h}{2},y_{n}+\frac{h}{2}K_{1}) [/math]
- [math] K_{3}=f(t_{n}+\frac{h}{2},y_{n}+\frac{h}{2}K_{2}) [/math]
- [math] K_{4}=f(t_{n}+h,y_{n}+K_{3}h) [/math]
- [math] y_{n+1}=y_{n}+\frac{h}{6}(K_{1}+2K_{2}+2K_{3}+K_{4}) [/math]
3.3 Gráficos del modelo de epidemias
Para realizar un gráfico del modelo utilizaremos los métodos numéricos explicados anteriormente, además veremos el número de infectados al final de la sexta semana y el tiempo necesario para que el número de infectados llegue a 400000. Lo haremos con el siguiente código MATLAB:
t0 = 0;h = 0.01;
r = [];
y = [];
y(1) = 200;
r(1) = 200;
i = 1;
c = 0.92;
N = 500000;
while y < 400000 && r < 400000
% Método Runge-Kutta
k1 = (c*y(i))/N*(N-y(i));
k2 = (c*(y(i)+1/2*k1*h))/N*(N-(y(i)+1/2*k1*h));
k3 = (c*(y(i)+1/2*k2*h))/N*(N-(y(i)+1/2*k2*h));
k4 = (c*(y(i)+k3*h))/N*(N-(y(i)+h*k3));
y(i+1) = y(i) + (h/6)*(k1+2*k2+2*k3+k4);
% Método de Heun
K1 = (c*r(i))/N*(N-r(i));
K2 = (c*(r(i)+K1*h))/N*(N-(r(i)+K1*h));
r(i+1) = r(i) + (h/2)*(K1+K2);
i = i + 1;
end
tN1 = h*(length(y)-1);
tN2 = h*(length(r)-1);
t1 = t0:h:tN1;
t2 = t0:h:tN2;
Semana_maxima_infeccion = [t1(end) t2(end)]
Valor_aproximado = [y(8),r(8)]
hold on
plot(t1,y)
plot(t2,r,'r')
legend('Runge Kutta','Heun','Location','best')
hold offAl final de la sexta semana el numero de infectados utilizando ambos métodos da el mismo resultado que es de 213.298.
Cuando llega a 400.000 el número de infectados, nos encontraríamos aproximadamente en la semana 10 sea cual sea el método numérico que se utilice.
Podemos ver que ambos métodos dan prácticamente el mismo gráfico, se puede apreciar un aumento continuo en el número de infectados por lo que no tiende a desaparecer ni tampoco se aprecia un momento de máximo número de infectados ya que no para de crecer.
4 Dificultad de contagio
Vamos a suponer que si durante la primera semana, el número de infectados no supera los 10.000, entonces a lo largo de las semanas los el número de individuos infectados irá disminuyendo. De lo contrario, si durante la primera semana hay más de 10.000 infectados entonces aumentará a lo largo del tiempo. Eso se debe por la dificultad que pueda tener la enfermedad de transmitirse. Por ello hay que modificar la ecuación diferencial inicial y pasar a la siguiente:
- [math] I'(t)= cI(t)(1-\frac{I(t)}{N})(\frac{I(t)}{M}-1) [/math]
Ahora con esta ecuación diferencial vamos a obtener la función I(t) utilizando el método de Heun, y lo vamos a hacer para [math] I(0) = 9700 [/math] , [math] I(0) = 10200 [/math] y [math] I(0) = 30000 [/math]. Lo haremos con el siguiente código MATLAB:
t0 = 0; h = 0.01; tN = 10; N = 500000; M = 10000; c = 0.92;
t = t0:h:tN;
n = round((tN-t0)/h);
y = zeros(1,n+1);
y0 = 9700;
y(1) = y0;
for i = 1:n
% Método de Heun
K1 = (c*y(i))/N*(N-y(i))*((y(i)/M)-1);
K2 = (c*(y(i)+K1*h))/N*(N-(y(i)+K1*h))*(((y(i)+K1*h)/M)-1);
y(i+1) = y(i) + (h/2)*(K1+K2);
end
plot(t,y)
% Para I(0) = 10200 y I(0) = 30000 sólo habría que cambiar en y0
% el resto del código sería igual
En el primer gráfico podemos ver que el número de infectados tiende a desaparecer ya que el número inicial es inferior a 10.000, por lo que el número de infectados máximo en ese caso es de 9.700 personas enfermas que es el número de contagiados que habían al principio. En cambio, en los otros dos gráficos tiende a crecer el número de individuos enfermos, por haber más de 10.000 individuos infectados con la única diferencia de que en el tercer gráfico crece más rápido que el segundo por haber más individuos infectados al principio en el tercer gráfico que en el segundo. Y el número máximo de individuos enfermos podríamos decir que es de toda la población, es decir 500.000 durante un tiempo indefinido, porque el modelo no contempla la posibilidad que algunos de sus habitantes se curen durante la epidemia.
5 Modelo de epidemias generalizado
Vamos a generalizar el modelo de epidemias planteado anterior. Además del número de infectados I(t) vamos a incluir las siguientes variables :
- S(t): Son los individuos susceptibles a infectarse.
- R(t): Serían los individuos curados, en cuarentena, o muertos. Es decir, las personas que no provocarían un aumento de infectados.
Ahora para modelizar vamos a plantearnos las siguientes ecuaciones diferenciales:
- [math] S'(t) = -rS(t)I(t) [/math]
- [math] I'(t) = rS(t)I(t) - aI(t) [/math]
- [math] R'(t) = aI(t) [/math]
- [math] S(0) = S_{0}, I(0) = I_{0}, R(0) = 0 [/math]
Ahora vamos a demostrar que la suma de las funciones planteadas es igual a la población total infectada, es decir
- [math] S(t) + I(t) + R(t) = S_{0} + I_{0} = N [/math]
Para ello vamos a sumar las ecuaciones diferenciales que se nos plantean y luego se integrarían, es decir
- [math] S'(t) + I'(t) + R'(t) = (- rS(t)I(t)) + (rS(t)I(t) - aI(t)) + (aI(t)) = 0 [/math]
- [math] S(t) + I(t) + R(t) = Constante [/math]
- [math] S(0) + I(0) + R(0) = S_{0} + I_{0} + 0 = N [/math]
5.1 Demostración analítica
Ahora vamos a demostrar que si [math] S_{0}\le\frac{a}{r} [/math] entonces [math] I(t)\le I_{0} [/math], siendo [math] \frac{a}{r}[/math] la tasa de remoción. Para ello vamos a sustituirlo en I'(t) del siguiente modo
- [math] I'(0)\le rS_{0}I_{0} - aI_{0}\le r\frac{a}{r}I_{0} - aI_{0}\le aI_{0} - aI_{0}\le 0 [/math]
En cambio si [math] S_{0}\gtfrac{a}{r} [/math] el número de individuos enfermos aumentará.
5.2 Gráficos generalizados
Ahora vamos a representar el modelo generalizado usando el método de Euler y Runge-Kutta para distintos [math] S_{0} [/math], [math] I_{0} [/math] y h. Lo hacemos con el siguiente código MATLAB:
% Datos. Irán cambiando según el caso
t0 = 0; tN = 4; S0 = 159985; I0 = 15; r = 0.0000218; a = 0.341; h = 0.1;
N = round((tN-t0)/h);
t = t0:h:tN;
S1 = zeros(1,N+1); % S' Euler
I1 = zeros(1,N+1); % I' Euler
R1 = zeros(1,N+1); % R' Euler
S2 = zeros(1,N+1); % S' Runge-Kutta
I2 = zeros(1,N+1); % I' Runge-Kutta
R2 = zeros(1,N+1); % R'Runge-Kutta
S1(1) = S0;
I1(1) = I0;
R1(1) = 0;
S2(1) = S0;
I2(1) = I0;
R2(1) = 0;
for i= 1:N
S1(i+1) = S1(i) + h*(-r*S1(i)*I1(i)); % Euler S'
I1(i+1) = I1(i) + h*(r*S1(i)*I1(i)-a*I1(i)); % Euler I'
R1(i+1) = R1(i) + h*(a*I1(i)); % Euler R'
% Runge-Kutta
k1_S = -r*S2(i)*I2(i);
k1_I = r*S2(i)*I2(i)-a*I2(i);
k1_R = a*I2(i);
k2_S = -r*(S2(i) + 1/2*k1_S*h)*(I2(i) + 1/2*k1_I*h);
k2_I = r*(S2(i) + 1/2*k1_S*h)*(I2(i) + 1/2*k1_I*h)-a*(I2(i) + 1/2*k1_I*h);
k2_R = a*(I2(i) + 1/2*k1_I*h);
k3_S = -r*(S2(i) + 1/2*k2_S*h)*(I2(i) + 1/2*k2_I*h);
k3_I = r*(S2(i) + 1/2*k2_S*h)*(I2(i) + 1/2*k2_I*h)-a*(I2(i) + 1/2*k2_I*h);
k3_R = a*(I2(i) + 1/2*k2_I*h);
k4_S = r*(S2(i) + k3_S*h)*(I2(i) + k3_I*h);
k4_I = r*(S2(i) + k3_S*h)*(I2(i) + k2_I*h)-a*(I2(i) + 1/2*k2_I*h);
k4_R = a*(I2(i) + k3_I*h);
S2(i+1) = S2(i) + (h/6)*(k1_S+2*k2_S+2*k3_S+k4_S); % Runge-Kutta S'
I2(i+1) = I2(i) + (h/6)*(k1_I+2*k2_I+2*k3_I+k4_I); % Runge-Kutta I'
R2(i+1) = R2(i) + (h/6)*(k1_R+2*k2_R+2*k3_R+k4_R); % Runge-Kutta R'
end
% Gráfico
hold on
plot(t,S1,'r','Linewidth',6)
plot(t,I1,'g','Linewidth',6)
plot(t,R1,'y','Linewidth',6)
plot(t,S2,'b','Linewidth',6)
plot(t,I2,'m','Linewidth',6)
plot(t,R2,'c','Linewidth',6)
legend('S´ Euler','I´ Euler','R´ Euler','S´ Runge-Kutta','I´ Runge-Kutta','R´ Runge-Kutta','Location','best')
hold off
En los tres casos que se nos plantean con distintos valores iniciales vemos que el número e individuos susceptibles a infectarse (S) se va reduciendo a medida que los individuos infectados (I) e individuos removidos (R) va aumentando, por lo que las variables que hemos utilizado se ajusta a los valores numéricos. Si vemos en un periodo medio de infección que se encontraría alrededor de tres semanas se puede apreciar cambios en todos los casos ya sea porque el número de susceptibles empiece a reducirse o porque dos de las tres variables se igualen en ese momento aproximadamente.
5.3 Trayectorias
Si dividimos I'(t) entre S'(t) tendríamos la siguiente ecuación diferencial.
- [math] \frac{dI}{dS} = -1 + \frac{a}{rS(t)} [/math]
- [math] I(t) + S(t) - \frac{a}{r}\log{S(t)} = k [/math]
En todas las trayectorias podemos apreciar un máximo, que al ser en función de I(t) sería el punto de mayor infectados, teniendo una forma de parábola aproximadamente en el que aumenta rápidamente el número de infectados y después del máximo, disminuye rápidamente hasta que el número de infectados tiende a cero.
5.4 Máximo número de infectados
Para terminar vamos a hallar el valor máximo de I(t) y el momento que sería el máximo número de infectados para el caso [math] S_{0} = 159985 [/math] . Daría de forma numérica 107988 infectados de forma numérica entre la tercera y cuarta semana. En cambio de forma analítica vemos que utilizando la siguiente expresión
- [math] I(máximo) = N - \frac{a}{r} + \frac{a}{r}\log\frac{a}{rS_{0}} = 500000 - \frac{0.341}{0.0000218} + \frac{0.341}{0.0000218}\log\frac{0.341}{0.0000218\cdot159985} = 447987.9965 [/math]
