Modelos Epidemológicos Grupo 4-C

De MateWiki
Revisión del 00:49 5 mar 2015 de Rojasrivero (Discusión | contribuciones) (→‎Método de Heun)

Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Modelos Epidemológicos. Grupo 20-C
Asignatura Ecuaciones Diferenciales
Curso Curso 2014-15
Autores Angela Béjar Gómez, Eduardo Bonet García, Gonzalo Ubeda, Elisa Pamplona Frey, Alberto Rojas, Fernando Sancha Domínguez
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

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] 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. Según se interpreta en la gráfica(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 codigo 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 sera el tiempo que tarde en erradicarse.

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.

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 requiere el uso de un código Matlab.

Método Runge-Kutta
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]

Para resolver este sistema hacemos uso del metodo de Heun, ayudandonos del siguiente codigo:


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
Tabla de datos obtenidos
%sistemas no lineales, no puedo usar matrices 
clear ; clc;
t0=0; 
tN=10;
h=0.0001;
N=round((tN-t0)/h); %n˙mero 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
Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Modelos_Epidemológicos_Grupo_4-C&oldid=25672»