Diferencia entre revisiones de «Modelos Epidemológicos Grupo 4-C»

Revisión del 12:03 5 mar 2015

En el desarrollo de una epidemia se distingue dos tipos de individuos: Los que ya han contraido de la enfermedad I, y los que son susceptibles de contraerla por encontrarse en zona de riesgo S. Consideramos las variables: t tiempo, S(t) población de individuos susceptibles a contraer la enfermedad, I(t) población de individuos infectados; y el sistema:: [math] \left \{ \begin{matrix} \frac{dS}{dt}=-aSI \\ \frac{dI}{dt}=aSI-bI-cI\end{matrix}\right. [/math] Suponemos que :

1. La población de personas infectadas se altera por el fallecimiento o la cura de las mismas. En ambos casos, la tasa de cambio depende del numero de personas infectadas;
2. La tasa de individuos que pasan de ser susceptibles a contraer la enfermedad a estar infectados es proporcional a la interacción entre el numero de individuos en ambas clases

Los parámetros a, b y c como:

1. "a":Proporción de la interacción entre personas susceptibles a ser infectados y personas ya infectadas.
2. "b":Proporción de personas que superan la enfermedad.
3. "c":Proporción de personas que mueren a causa de la enfermedad.

1 Método de Euler y Trapecio

Consideramos la ecuación: [math]\frac{dI}{dt}=aSI-bI-cI[/math] Con S=0, es decir, [math]\frac{dI}{dt}=-bI-cI[/math]. Conocidos los valores b=0.3y c=0.01, nos queda el problema de valor inicial siguiente: [math]\left \{ \begin{matrix}\frac{dI}{dt}=-bI-cI; \\ I_{0}=2000 \end{matrix}\right.[/math]: [math]b=0.3 \\ c=0.01[/math] Aplicamos el método de Euler para resolver la ecuación diferencial utilizamos un intervalo de tiempo de 0 a 50 porque a los 41 días el mínimo de infectados se reduce a cero. La gráfica muestra que el numero de infectados se reduce exponencialmente con el tiempo, en concreto,(ayudándonos de un bucle while de matlab) se observa que el número de infectados se reduce a la cuarta parte en 4,5 días. El código matlab utilizado es:

Comparación del método de Euler y del Trapecio

t0=0;
tN=50; 
% consideramos el intervalo de tiempo [0,50] ya que el numero 
%de infectados será cero a los 41 días
I0=2000;
S=0
h=0.1;
N=round((tN-t0)/h); 
t=t0:h:tN;
I=zeros(1,N+1); %euler 
I(1)=I0;
Z=zeros(1,N+1); % Trapecio
Z(1)=I0;
G=zeros(1,N+1); 
% Vector para cáculo del tiempo que tarda en reducirse I(i)a la cuarta parte
G(1)=I0;
a=0.003;
b=0.3;
c=0.01;
for i=1:N
 I(i+1)=I(i)+h*(a*S.*I(i)-b.*I(i)-c.*I(i)); %euler
 Z(i+1)=(Z(i).*(1+h/2*(-b-c)))/(1+h/2*(b+c)); %trapecio
end
i=1;
tg0=0;
tg(1)=tg0;
while G(i)>(I0/4)
 G(i+1)=G(i)+h*(a*S.*G(i)-b.*G(i)-c.*G(i));
%bucle que calcula el tiempo que tarda en reducirse I(i) a la cuarta parte
 tg(i+1)=tg(i)+h;
 i=i+1;
end
%sacamos tabla de resultados
[t',I',Z'] 
hold on
plot(t,I)
plot(t,Z,'r')
hold off
legend('Euler','Trapecio','Location','best');
disp('Tiempo final:')
disp(tg(end))

Consideramos ahora S=100 constante a lo largo del tiempo y a=0,003, hacemos uso de un código de Matlab, como el del apartado anterior para resolver esta ecuación y se interpreta en la gráfica que cuanto mayor sea la población susceptible a contraer la enfermedad, mayor será el tiempo que tarde en erradicarse. Calculamos a su vez cuanto tarda en reducirse el numero de infectados a la cuarta parte y resulta 138.6000 días.

Comparación del método de Euler y del Trapecio

t0=0;
tN=600; % consideramos el intervalo de tiempo [0,600]
I0=2000;
S=100
h=0.1;
N=round((tN-t0)/h); 
t=t0:h:tN;
I=zeros(1,N+1); %euler 
I(1)=I0;
Z=zeros(1,N+1); % Trapecio
Z(1)=I0;
G=zeros(1,N+1); 
% Vector para cáculo del tiempo que tarda en reducirse I(i)a la cuarta parte
G(1)=I0;
a=0.003;
b=0.3;
c=0.01;
for i=1:N
 I(i+1)=I(i)+h*(a*S.*I(i)-b.*I(i)-c.*I(i)); %euler
 Z(i+1)=(Z(i).*(1+h/2*(a*S-b-c)))/(1-h/2*(a*S-b-c)); %trapecio
end
i=1;
tg0=0;
tg(1)=tg0;
while G(i)>500
 G(i+1)=G(i)+h*(a*S.*G(i)-b.*G(i)-c.*G(i));
%bucle que calcula el tiempo que tarda en reducirse I(i)a la cuarta parte
 tg(i+1)=tg(i)+h;
 i=i+1;
end
%sacamos tabla de resultados
[t',I',Z'] 
hold on
plot(t,I)
plot(t,Z,'r')
hold off
legend('Euler','Trapecio','Location','best');
disp('Tiempo final:')
disp(tg(end))

Repetimos el proceso para S=200 constante a lo largo del tiempo y vemos en la gráfica que la ecuación ya no es valida para este valor de S debido a que sobrepasa el límite y la enfermedad no podrá ser erradicada, ya que la ecuación crece exponencialmente con el tiempo, es decir que el numero de infectados no para de crecer. Por tanto, no tiene sentido calcular el tiempo que tarda en reducirse el numero de infectados a la cuarta parte.

Comparación del método de Euler y del Trapecio

t0=0;
tN=600; % consideramos el intervalo de tiempo [0,600]
I0=2000;
S=200
h=0.1;
N=round((tN-t0)/h); 
t=t0:h:tN;
I=zeros(1,N+1); %euler 
I(1)=I0;
Z=zeros(1,N+1); % Trapecio
Z(1)=I0;
G=zeros(1,N+1); 
% Vector para cáculo del tiempo que tarda en reducirse I(i)
%a la cuarta parte
G(1)=I0;
a=0.003;
b=0.3;
c=0.01;
for i=1:N
 I(i+1)=I(i)+h*(a*S.*I(i)-b.*I(i)-c.*I(i)); %euler
 Z(i+1)=(Z(i).*(1+h/2*(a*S-b-c)))/(1-h/2*(a*S-b-c)); %trapecio
end
i=1;
tg0=0;
tg(1)=tg0;
while G(i)>500
 G(i+1)=G(i)+h*(a*S.*G(i)-b.*G(i)-c.*G(i));
%bucle que calcula el tiempo que tarda en reducirse I(i) 
%a la cuarta parte
 tg(i+1)=tg(i)+h;
 i=i+1;
end
%sacamos tabla de resultados
[t',I',Z'] 
hold on
plot(t,I)
plot(t,Z,'r')
hold off
legend('Euler','Trapecio','Location','best');
disp('Tiempo final:')
disp(tg(end))

2 Método Runge-Kutta

Haciendo uso de este método de cuarto orden para resolver la ecuación [math]\frac{dI}{dt}=-bI-cI[/math] Al comparar el resultado con el método de Euler observamos que los resultados para ambos métodos son muy próximos. No es aconsejable usar un metodo implicito como el trapezoidal debido a que es demasiado complejo despejar la ecuación analíticamente Por tanto se usa normalmente un método como el de Runge-Kutta. Según vemos en la gráfica la ecuación crece exponencialmente con el tiempo, esto es, que el numero de infectados no para de crecer por lo que la enfermedad no podrá ser erradicada

Comparación del método Runge-Kutta y de Euler

t0=0;
tN=25; 
% consideramos el intervalo de tiempo [0,50]
%ya que el numero de infectados será cero a los 41 días
I0=2000;
S=0
h=0.1;
N=round((tN-t0)/h); 
t=t0:h:tN;
I=zeros(1,N+1); %euler 
I(1)=I0;
a=0.003;
b=0.3;
c=0.01;
r=zeros(1,N+1); %RK4
r(1)=I0;
for i=1:N
 I(i+1)=I(i)+h*(a*S.*I(i)-b.*I(i)-c.*I(i)); %euler
 % RK4
 K1=a*S.*r(i)-b.*r(i)-c.*r(i);
 K2=-b.*(r(i)+(1/2).*K1.*h)-c.*(r(i)+(1/2).*K1.*h);
 K3=-b.*(r(i)+(1/2).*K2.*h)-c.*(r(i)+(1/2).*K2.*h);
 K4=-b.*(r(i)+K3*h)-c.*(r(i)+K3*h);
 r(i+1)= r(i)+h/6*(K1+2*K2+2*K3+K4) 
end

%sacamos tabla de resultados
[t',I',r'] 
hold on
plot(t,I)
plot(t,r,'r')
hold off
legend('Euler','Runge Kutta','Location','best');

3 Método de Heun

Suponemos ahora que [math]a\left ( t \right )=\frac{0.003}{1+t}[/math],lo que significa que la interaccion entre infectados y susceptibles disminuye a medida que aumenta el tiempo. Tambien tomamos [math] I_{O}=40[/math] y [math] S_{O}=1640[/math]. Con lo que el sistema nos quedaría asi:

[math] \left \{ \begin{matrix}\frac{dS}{dt}= \frac{-0.003}{1+t}\cdot S\cdot I \\ \frac{dI}{dt}= \frac{0.003}{1+t}\cdot S\cdot I-0.3\cdot I-0.01\cdot I\end{matrix}\right.[/math] En la gráfica vemos que el número de infectados alcanza un máximo entorno a t=2, y luego decrece exponencialmente con el tiempo, sin embargo el número de susceptibles empieza en 1640 y decrece exponencialmente con el tiempo más rápido que el número de infectados

Para resolver este sistema hacemos uso del metodo de Heun, ayudandonos del siguiente código :

Método de Heun

%sistemas no lineales, no puedo usar matrices 
clear ; clc;
t0=0; 
tN=10;
h=0.0001;
N=round((tN-t0)/h); %numero de subintervalos

t=t0:h:tN;
t(1)=t0;
I=zeros(1,N+1); %INFECTADOS
S=zeros(1,N+1); %SUBCEPTIBLES 
I(1)=40;
S(1)=1640;
a=zeros(1,N+1);
a(1)=0.003/(1+t0);

for i=1:N
    a(i+1)=0.003/(1+t(i));
    K1S=-a(i).*S(i).*I(i); 
    K2S=-a(i)*S(i).*I(i)+h*K1S;
    S(i+1)=S(i)+(h/2).*(K1S+K2S);
    K1I=a(i).*S(i)*I(i)-0.3.*I(i)-0.01*I(i);
    K2I=a(i).*S(i)*I(i)-0.31*I(i)+K1I*h;
    I(i+1)=I(i)+(h/2)*(K1I+K2I);
end
[S,I,t]';
hold on 
plot(t,S);
plot(t,I,'r');
legend('SUSCEPTIBLES','INFECTADOS');
hold off

Basándonos en los resultados obtenidos, resolvemos el sistema para diferentes valores de la constante a. Con el programa de matlab desarrollado a continuación obtenemos una tabla con todos los valores de tiempo en los que se alcanza el número máximo de infectados para cada valor de a. A continuación se muestra el programa utilizado así como la tabla obtenida con este.

Tabla de datos obtenidos

% nos dicen que Imax se alcanza en t=5 en una experiencia 
% Resolver el sistema para a(0.0005,0.002) y ver cual se ajusta a la exp
%usar un paso de 10^-4 y usar bucle para el valor max 
%nos quedamos con el valor de a que da un maximo mas cerca de t=5 
clear; clc; 
h=0.0001;
t0=0; 
tN=5;
N=(tN-t0)/h; %numero de subintervalos
% Definimos la variable independiente
t=t0:h:tN;
t(1)=t0;
I=zeros(1,N+1); %INFECTADOS
S=zeros(1,N+1); %SUBCEPTIBLES 
I(1)=40;
S(1)=1640;
c=[]
for a=0.0005:h:0.002
    for i=1:N
        K1S=a*S(i)*I(i);
        K2S=a*(S(i)+K1S*h)*(I(i)+h*K1S);
        K1I=a*S(i)*I(i)-0.3*I(i)-0.01*I(i);
        K2I=a*(S(i)+h*K1S)*(I(i)+K1I*h)-0.3*(I(i)+K1I*h)-(0.01*(I(i)+K1I*h));
        I(i+1)=I(i)+(h/2)*(K1I+K2I);
        S(i+1)=S(i)+(h/2)*(K1S+K2S);
        [v,w]=max(I);
    end
    c=[c;[a,t(w)]];
end