Modelos Epidemológicos Grupo 4-C
| 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:
- "a":Proporción de la interacción entre personas susceptibles a ser infectados y personas ya infectadas.
- "b":Proporción de personas que superan la enfermedad.
- "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:
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.
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.
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.
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.
%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 offTEXTO
%PROBLEMA 7, TRABAJO 4
% 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; %n˙mero 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
