Diferencia entre revisiones de «Modelos Epidemológicos Grupo 4-C»
(→Método Runge-Kutta) |
(→Método de Euler y Trapecio) |
||
| Línea 123: | Línea 123: | ||
Z(1)=I0; | Z(1)=I0; | ||
G=zeros(1,N+1); | G=zeros(1,N+1); | ||
| − | % Vector para cáculo del tiempo que tarda en reducirse I(i) | + | % Vector para cáculo del tiempo que tarda en reducirse I(i) |
| + | %a la cuarta parte | ||
G(1)=I0; | G(1)=I0; | ||
a=0.003; | a=0.003; | ||
| Línea 137: | Línea 138: | ||
while G(i)>500 | while G(i)>500 | ||
G(i+1)=G(i)+h*(a*S.*G(i)-b.*G(i)-c.*G(i)); | 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) | + | %bucle que calcula el tiempo que tarda en reducirse I(i) |
| + | %a la cuarta parte | ||
tg(i+1)=tg(i)+h; | tg(i+1)=tg(i)+h; | ||
i=i+1; | i=i+1; | ||
Revisión del 13:58 4 mar 2015
| 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]\frac{dS}{dt}=-aSI[/math]: [math]\frac{dI}{dt}=aSI-bI-cI[/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]\frac{dI}{dt}=-bI-cI;[/math]: [math]b=0.3 \\ c=0.01 \\ I_{0}=2000[/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]\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[/math]
Para resolver este sistema hacemos uso del metodo de Heun, ayudandonos del siguiente codigo:
Basandonos en los resultados obtenidos, resolvemos el sistema para diferentes valores de la constante a.
Vemos que el que mas se ajusta más a la experiencia es el...
