Usuario:Carlosferber

Trabajo realizado por estudiantes
Título	Modelo para epidemias. Grupo 23
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2014-15
Autores	Álvaro Roales Blanco, Carlos Fernández Bermejo, Adrián Gómez Apiñániz, Marta Torra Escánez, Ricardo Mazón Cabrera
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

El objeto de estudio de nuestro trabajo es el modelo de expansión de una epidemia a lo largo de un tiempo. Usaremos las siguientes hipótesis:

1. La población (constante) se indicará como N en el que todos serán susceptibles a la enfermedad.
2. La duración de la enfermedad es larga de modo que los individuos no se recuperarán durante el periodo de estudio.
3. La unidad de tiempo escogida será la semana.
4. Durante cada unidad de tiempo cada persona infectada tiene c contactos.

2 Determinación del número de contactos c

El primer paso a realizar es determinar el número de contactos que tiene una persona infectada en la unidad de tiempo elegida, es decir, la semana. Emplearemos la siguiente ecuación diferencial que nos muestra el enunciado del trabajo:

[math]I'(t) = \frac{c}{N}I(t)(N-I(t))[/math]

Asignando a I(t) el número de personas infectadas.

c es un número de dos decimales comprendido entre [math] 0.01 [/math] y [math] 0.99 [/math]. Siendo N igual a 500000 con las siguientes condiciones iniciales: [math] I(0)=200 [/math] y [math] I(1)=500 [/math] Para determinar c imponemos la siguiente condición, que [math] |I(1)-500| [/math] tenga el menor valor posible, es decir, aquel valor que minimice lo posible el error absoluto. Para hallar ese valor lo hallaremos numéricamente con el siguiente programa escrito en MATLAB, basado en el método de Euler para la resolución de ecuaciones diferenciales:

y0 = 200; % Número inicial de infectados
y = []; 
c = 1/100:0.01:99/100; % Posibles valores de c
N = 500000; % Población total
k = length(c);
Cminimo = c(1);
for i = 1:k
   Yfin = N;
   Y = N*y0/(y0+(N-y0).*exp(-c(i))); % Modelo para t=1
   y = [y Y];
   if abs(Y-500) <= abs(Yfin-500)
       Cminimo = c(i);
       end
  end
min(abs(y-500)) % Error mínimo absoluto cometido  
C_min
plot(c,abs(y-500),'r')

Al ejecutar el programa anteriormente escrito obtenemos el valor c minimizando así el error que será igual a [math] 0.92 [/math] Por otra parte hemos dibujado un gráfico para ver la variación del error para cada valor de c:

Evolución del error según los diferentes valores de c

3 Estudio de la evolución de la enfermedad

A continuación haremos el estudio de la evolución de la epidemia a lo largo de las semanas mediante los dos siguientes métodos.

Primero explicaremos el Método de Heun:

Consiste en la aproximación a la pendiente mediante la aplicación de dos derivadas del intervalo, una en el punto inicial y otra en punto final. La aproximación mejorada de la pendiente será el promedio de las dos derivadas, 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

En segundo lugar expondremos el 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]

Mediante los métodos expuestos anteriormente podemos determinar el gráfico de evolución de los infectados conforme avanzan las semanas. Para ello aplicaremos el siguiente código MATLAB:

Evolución de los infectados conforme avanzan las semanas mediante el método de Heun y el método de Runge-Kutta

t0 = 0;h = 0.01;
r = [];
y = [];
y(1) = 200;
r(1) = 200;
i = 1;
c = 0.92;
N = 500000;
for i=1:N;
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
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,'k')
plot(t2,r,'g')
legend('Runge Kutta','Heun','Location','best')
hold off

Como podemos observar en el gráfico adjunto, los infectados no tienden a desaparecer sino a crecer exponencialmente. Como es una función creciente no posee máximo, sin embargo como el número de habitantes está acotados el máximo se produce cuando todos los habitantes están infectados, es decir 500.000. Este fenómeno se produce pasada la semana número 10.

Al finalizar la sexta semana el número aproximado de infectados es 110.000.

A partir de la décima semana, el número de infectados es mayor a 400.000.

4 Susceptibilidad al contagio

t0 = 0;% Tiempo inicial
h = 0.01;% Paso del método numérico 
tN = 10; % Última semana
N = 500000; % Tamaño de población
M = 10000; % Umbral de contagio
c = 0.92; % Susceptibilidad al contagio
t = t0:h:tN; % Valores del tiempo
n = round((tN-t0)/h); % Número de valores del tiempo(instantes de tiempo)
y = zeros(1,n+1); % Valores de personas contagiadas
y1(1) = 9700;
y2(1)=10200;
y3(1)=30000;
for i = 1:n
% Método de Heun
K11 = (c*y1(i))/N*(N-y1(i))*((y1(i)/M)-1);
K21 = (c*(y1(i)+K11*h))/N*(N-(y1(i)+K11*h))*(((y1(i)+K11*h)/M)-1);
y1(i+1) = y1(i) + (h/2)*(K11+K21);

K12 = (c*y2(i))/N*(N-y2(i))*((y2(i)/M)-1);
K22 = (c*(y2(i)+K12*h))/N*(N-(y2(i)+K12*h))*(((y2(i)+K12*h)/M)-1);
y2(i+1) = y2(i) + (h/2)*(K12+K22);

K13 = (c*y3(i))/N*(N-y3(i))*((y3(i)/M)-1);
K23 = (c*(y3(i)+K13*h))/N*(N-(y3(i)+K13*h))*(((y3(i)+K13*h)/M)-1);
y3(i+1) = y3(i) + (h/2)*(K13+K23);
end

figure(1)
axis equal
plot(t,y1)
figure(2)
axis equal
plot(t,y2)
figure(3)
axis equal
plot(t,y3)

4.1 Modelo SIR

En primer lugar definimos las variables:

-S: Número de personas susceptibles a la enfermedad.

-I: Número de infectados

-R: Número de personas en cuarentena, muertos o que han superado o son inmunes a la enfermedad.

Definimos S(t), I(t), R(t), cantidad de miembros de cada clase al transcurrir un tiempo t.

4.2 df

El siguiente código Matlab, que utiliza el método Euler y Runge-Kutta, se ejecuta para distintos valores de So, Io y h. Se encuentran ordenados con el mismo So e Io pero variando la h, tomando los valores 0.1, 0.01, 0.001.

% 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);
I1 = zeros(1,N+1); 
R1 = zeros(1,N+1); 
S2 = zeros(1,N+1); 
I2 = zeros(1,N+1); 
R2 = zeros(1,N+1); 
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));
 I1(i+1) = I1(i) + h*(r*S1(i)*I1(i)-a*I1(i));
 R1(i+1) = R1(i) + h*(a*I1(i)); 
 
 k1S = -r*S2(i)*I2(i); 
 k1I = r*S2(i)*I2(i)-a*I2(i);
 k1R = a*I2(i);
 k2S = -r*(S2(i) + 1/2*k1S*h)*(I2(i) + 1/2*k1I*h);
 k2I = r*(S2(i) + 1/2*k1S*h)*(I2(i) + 1/2*k1I*h)-a*(I2(i) + 1/2*k1I*h);
 k2R = a*(I2(i) + 1/2*k1I*h);
 k3S = -r*(S2(i) + 1/2*k2S*h)*(I2(i) + 1/2*k2I*h);
 k3I = r*(S2(i) + 1/2*k2S*h)*(I2(i) + 1/2*k2I*h)-a*(I2(i) + 1/2*k2I*h);
 k3R = a*(I2(i) + 1/2*k2I*h);
 k4S = r*(S2(i) + k3S*h)*(I2(i) + k3_I*h);
 k4I = r*(S2(i) + k3S*h)*(I2(i) + k2_I*h)-a*(I2(i) + 1/2*k2I*h);
 k4R = a*(I2(i) + k3I*h);
S2(i+1) = S2(i) + (h/6)*(k1S+2*k2S+2*k3S+k4S); 
I2(i+1) = I2(i) + (h/6)*(k1I+2*k2I+2*k3I+k4I); 
R2(i+1) = R2(i) + (h/6)*(k1R+2*k2R+2*k3R+k4R); 
 end


hold on
plot(t,S1,'b','Linewidth',10)
plot(t,I1,'k','Linewidth',10)
plot(t,R1,'g','Linewidth',10)
plot(t,S2,'r','Linewidth',10)
plot(t,I2,'b','Linewidth',10)
plot(t,R2,'m','Linewidth',10)
legend('S´ Euler','I´ Euler','R´ Euler','S´ Runge-Kutta','I´ Runge-Kutta','R´ Runge-Kutta','Location','best')
hold off

Gráfica con So= 159985; Io= 15; h= 0.1

Gráfica con So= 159985; Io= 15; h= 0.01

Gráfica con So= 159985; Io= 15; h=0.001

Gráfica con So= 140000; Io= 20000; h=0.1

Gráfica con So= 140000; Io= 20000; h=0.01

Gráfica con So= 140000; Io= 20000; h=0.001

Gráfica con So= 159999; Io= 1; h=0.1

Gráfica con So= 159999; Io= 1; h=0.01

Gráfica con So= 159999; Io= 1; h=0.001

En los gráficos se observa que los distintos valores iniciales de S (número de susceptibles a contraer la enfermedad), se ven reducidos a medida que aumenta el valor de I (número de infectados) , por lo tanto estas datos se corresponden con los valores númericos del enunciado. Tambien cabe destacar que a medida tomamos un intervalo h más pequeño, la exactitud es mayor por lo que el error cometido es menor.