Diferencia entre revisiones de «Modelo para epidemias (Grupo 6-A)»

Revisión actual del 20:16 10 mar 2015

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:

1. La población será un número fijo N en el que todos serán susceptibles a la enfermedad.
2. Los individuos infectados no se curarán durante el periodo de la enfermedad y serán contagiosos.
3. La unidad de tiempo será la semana.
4. 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 nos 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:

Error absoluto en I(1) para cada 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:

1. [math] K_{1}=f(t_{n},y_{n}) [/math]
2. [math] K_{2}=f(t_{n}+h,y_{n}+K_{1}h) [/math]
3. [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:

1. [math] K_{1}=f(t_{n},y_{n}) [/math]
2. [math] K_{2}=f(t_{n}+\frac{h}{2},y_{n}+\frac{h}{2}K_{1}) [/math]
3. [math] K_{3}=f(t_{n}+\frac{h}{2},y_{n}+\frac{h}{2}K_{2}) [/math]
4. [math] K_{4}=f(t_{n}+h,y_{n}+K_{3}h) [/math]
5. [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:

Modelo de epidemias con el método de Heun y Runge-Kutta

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 off

Al 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, habrán pasado aproximadamente 10 semanas desde el comienzo de la infección 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

modelo para I(0) = 9700

modelo para I(0) = 10200

modelo para I(0) = 30000

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 :

1. S(t): Son los individuos susceptibles a infectarse.
2. 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}\gt\frac{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

Para [math] S_{0} = 159985 [/math] e [math] I_{0} = 15 [/math]

Para [math] S_{0} = 140000 [/math] e [math] I_{0} = 20000 [/math]

Para [math] S_{0} = 159999 [/math] e [math] I_{0} = 1 [/math]

Para [math] S_{0} = 159985 [/math] e [math] I_{0} = 15 [/math]

Para [math] S_{0} = 140000 [/math] e [math] I_{0} = 20000 [/math]

Para [math] S_{0} = 159999 [/math] e [math] I_{0} = 1 [/math]

Para [math] S_{0} = 159985 [/math] e [math] I_{0} = 15 [/math]

Para [math] S_{0} = 140000 [/math] e [math] I_{0} = 20000 [/math]

Para [math] S_{0} = 159999 [/math] e [math] I_{0} = 1 [/math]

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]

Utilizando el método de Euler

Utilizando el método de Runge-Kutta

El intervalo (S,I) en el que se encuentran las trayectorias para [math] S_{0} =159985 [/math] e [math] I_{0} = 15 [/math].

El intervalo (S,I) en el que se encuentran las trayectorias para [math] S_{0} =140000 [/math] e [math] I_{0} = 20000 [/math].

El intervalo (S,I) en el que se encuentran las trayectorias para [math] S_{0} =159999 [/math] e [math] I_{0} = 1 [/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 447988 infectados de forma numérica casi 4 semanas después de que la infección comenzara. 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]