Diferencia entre revisiones de «Usuario:Carlosferber»

Revisión del 11:14 4 mar 2015

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]

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:

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);
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 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 tiende a desaparecer sino a crecer exponencialmente. Como es una función creciente no posee máximo. 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

Evolución gráfica de los infectados conforme avanzan las semanas siendo la gráfica situada arriba a la izq el modelo I(0)=9700, arriba a la derecha el modelo I(0)=10200 y la restante el modelo I(0)=30000

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

% 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

dhasjkdahasjk