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

De MateWiki
Saltar a: navegación, buscar
 
(No se muestran 75 ediciones intermedias de 3 usuarios)
Línea 1: Línea 1:
−
{| class="wikitable"
+
 
−
|-
+
{{ TrabajoED | Modelos Epidemológicos. Grupo 20-C | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | Angela Béjar Gómez, Eduardo Bonet García, Gonzalo Ubeda, Elisa Pamplona Frey, Alberto Rojas, Fernando Sancha Domínguez }}
−
| Alberto
+
−
|-
+
−
| Texto de celda
+
−
|-
+
−
| Texto de celda
+
−
|-
+
−
| Texto de celda
+
−
|-
+
−
| Texto de celda
+
−
|-
+
−
| Texto de celda
+
−
|}
+
 
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.  
 
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::  
 
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> \left \{ \begin{matrix} \frac{dS}{dt}=-aSI \\
−
<math>\frac{dI}{dt}=aSI-bI-cI</math>
+
\frac{dI}{dt}=aSI-bI-cI\end{matrix}\right. </math>
 +
Suponemos que :
 +
# 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;
 +
# 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:
 
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.
+
#"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.
+
#"b":Proporción de personas que superan la enfermedad.
−
"c":Proporción de personas que mueren a causa de la enfermedad.
+
#"c":Proporción de personas que mueren a causa de la enfermedad.
 +
== 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:
 +
[[Archivo:fotofoto2.png|400px|thumb|right|Comparación del método de Euler y del Trapecio]]
 +
  {{matlab|codigo=
 +
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.
 +
 
 +
[[Archivo:Graficazorra.png|400px|thumb|right|Comparación del método de Euler y del Trapecio]]
 +
{{matlab|codigo=
 +
 +
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.
 +
[[Archivo:Graficazorra2.png|400px|thumb|right|Comparación del método de Euler y del Trapecio]]
 +
{{matlab|codigo=
 +
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))
 +
}}
 +
== Modelo Completo==
 +
En este caso obtendremos valores de S e I en función del tiempo por el método de Euler resolviendo el sistema
 +
debido a que el sistema no es lineal aplicamos Euler para cada una de las ecuaciones, lo que nos dará el
 +
resultado gráficamente de S e I y sus evoluciones en el tiempo, además utilizaremos cuatro longitudes de paso
 +
lo cual nos reportará cuatro gráficas en cada caso. Lo resolvemos para unos valores iniciales I(0)=20 S(0)=800
 +
e I(0)=40 S(0)=10000.
 +
En el primer caso utilizamos el siguiente código MatLab:
 +
[[Archivo:graficafer1.png|400px|thumb|right|Modelo Completo]]
 +
 
 +
{{matlab|codigo=
 +
 
 +
% Euler
 +
clear all
 +
%DATOS DEL PROBLEMA
 +
t0=0;
 +
tN=40;
 +
I0=20;
 +
S0=800;
 +
h1=0.1;
 +
h2=0.01;
 +
h3=0.001;
 +
h4=0.0001;
 +
a=0.003;
 +
b=0.3;
 +
c=0.01;
 +
%Calculamos número de subintervalos
 +
N1=round((tN-t0)/h1);
 +
N2=round((tN-t0)/h2);
 +
N3=round((tN-t0)/h3);
 +
N4=round((tN-t0)/h4);
 +
% Definimos la variable independiente
 +
t1=t0:h1:tN;
 +
t2=t0:h2:tN;
 +
t3=t0:h3:tN;
 +
t4=t0:h4:tN;
 +
%Ahora vamos a guardar los valores de la solución aproximada en el vector y
 +
I=zeros(1,N1+1); 
 +
S=zeros(1,N1+1);
 +
F=zeros(1,N2+1); 
 +
R=zeros(1,N2+1);
 +
A=zeros(1,N3+1); 
 +
L=zeros(1,N3+1);
 +
G=zeros(1,N4+1); 
 +
U=zeros(1,N4+1);
 +
I(1)=I0;
 +
S(1)=S0;
 +
F(1)=I0;
 +
R(1)=S0;
 +
A(1)=I0;
 +
L(1)=S0;
 +
G(1)=I0;
 +
U(1)=S0;
 +
for i=1:N1
 +
    S(i+1)=S(i)+h1.*(-a.*S(i).*I(i)); %eulerS
 +
    I(i+1)=I(i)+h1.*(a.*S(i).*I(i)-b.*I(i)-c.*I(i)) ; %eulerI
 +
end
 +
for i=1:N2
 +
    R(i+1)=R(i)+h2.*(-a.*R(i).*F(i)); %eulerS
 +
    F(i+1)=F(i)+h2.*(a.*R(i).*F(i)-b.*F(i)-c.*F(i)) ; %eulerI
 +
end
 +
for i=1:N3
 +
    L(i+1)=L(i)+h3.*(-a.*L(i).*A(i)); %eulerS
 +
    A(i+1)=A(i)+h3.*(a.*L(i).*A(i)-b.*A(i)-c.*A(i)) ; %eulerI
 +
end
 +
for i=1:N4
 +
    U(i+1)=U(i)+h4.*(-a.*U(i).*G(i)); %eulerS
 +
    G(i+1)=G(i)+h4.*(a.*U(i).*G(i)-b.*G(i)-c.*G(i)) ; %eulerI
 +
end
 +
J1=max(I);
 +
J2=max(F);
 +
J3=max(A);
 +
J4=max(G);
 +
disp('El valor maximo en la grafica 1 es')
 +
J1
 +
  disp('El valor maximo en la grafica 2 es')
 +
J2
 +
  disp('El valor maximo en la grafica 3 es')
 +
J3
 +
  disp('El valor maximo en la grafica 4 es')
 +
J4
 +
subplot(2,2,1)
 +
hold on
 +
plot(t1,I);
 +
plot(t1,S,'r');
 +
legend('INFECTADOS','SUSCEPTIBLES')
 +
hold off
 +
subplot(2,2,2)
 +
hold on
 +
plot(t2,F);
 +
plot(t2,R,'g');
 +
legend('INFECTADOS','SUSCEPTIBLES')
 +
hold off
 +
subplot(2,2,3)
 +
hold on
 +
plot(t3,A);
 +
plot(t3,L,'y');
 +
legend('INFECTADOS','SUSCEPTIBLES')
 +
hold off
 +
subplot(2,2,4)
 +
hold on
 +
plot(t4,G);
 +
plot(t4,U,'c');
 +
legend('INFECTADOS','SUSCEPTIBLES')
 +
hold off }}
 +
 
 +
Éste programa de MatLab nos dará además el valor máximo de población infectada en cada gráfica que son:
 +
1.I=517.45
 +
2.I=506.37
 +
3.I=505.3
 +
4.I=505.19
 +
En el segundo caso utilizamos el código MatLab:
 +
[[Archivo:graficafer2.png|400px|thumb|right|Modelo Completo]]
 +
{{matlab|codigo=
 +
% Euler
 +
clear all
 +
%DATOS DEL PROBLEMA
 +
t0=0;
 +
tN=40;
 +
I0=40;
 +
S0=10000;
 +
h1=0.1;
 +
h2=0.01;
 +
h3=0.001;
 +
h4=0.0001;
 +
a=0.003;
 +
b=0.3;
 +
c=0.01;
 +
%Calculamos número de subintervalos
 +
N1=round((tN-t0)/h1); % también N=length(t)-1;
 +
N2=round((tN-t0)/h2); % también N=length(t)-1;
 +
N3=round((tN-t0)/h3); % también N=length(t)-1;
 +
N4=round((tN-t0)/h4); % también N=length(t)-1;
 +
% Definimos la variable independiente
 +
t1=t0:h1:tN;
 +
t2=t0:h2:tN;
 +
t3=t0:h3:tN;
 +
t4=t0:h4:tN;
 +
%Ahora vamos a guardar los valores de la solución aproximada en el vector y
 +
I=zeros(1,N1+1); 
 +
S=zeros(1,N1+1);
 +
F=zeros(1,N2+1); 
 +
R=zeros(1,N2+1);
 +
A=zeros(1,N3+1); 
 +
L=zeros(1,N3+1);
 +
G=zeros(1,N4+1); 
 +
U=zeros(1,N4+1);
 +
I(1)=I0;
 +
S(1)=S0;
 +
F(1)=I0;
 +
R(1)=S0;
 +
A(1)=I0;
 +
L(1)=S0;
 +
G(1)=I0;
 +
U(1)=S0;
 +
for i=1:N1
 +
    S(i+1)=S(i)+h1.*(-a.*S(i).*I(i)); %eulerS
 +
    I(i+1)=I(i)+h1.*(a.*S(i).*I(i)-b.*I(i)-c.*I(i)) ; %eulerI
 +
end
 +
for i=1:N2
 +
    R(i+1)=R(i)+h2.*(-a.*R(i).*F(i)); %eulerS
 +
    F(i+1)=F(i)+h2.*(a.*R(i).*F(i)-b.*F(i)-c.*F(i)) ; %eulerI
 +
end
 +
for i=1:N3
 +
    L(i+1)=L(i)+h3.*(-a.*L(i).*A(i)); %eulerS
 +
    A(i+1)=A(i)+h3.*(a.*L(i).*A(i)-b.*A(i)-c.*A(i)) ; %eulerI
 +
end
 +
for i=1:N4
 +
    U(i+1)=U(i)+h4.*(-a.*U(i).*G(i)); %eulerS
 +
    G(i+1)=G(i)+h4.*(a.*U(i).*G(i)-b.*G(i)-c.*G(i)) ; %eulerI
 +
end
 +
J1=max(I);
 +
J2=max(F);
 +
J3=max(A);
 +
J4=max(G);
 +
disp('El valor maximo en la grafica 1 es')
 +
J1
 +
  disp('El valor maximo en la grafica 2 es')
 +
J2
 +
  disp('El valor maximo en la grafica 3 es')
 +
J3
 +
  disp('El valor maximo en la grafica 4 es')
 +
J4
 +
subplot(2,2,1)
 +
hold on
 +
plot(t1,I);
 +
plot(t1,S,'r');
 +
legend('INFECTADOS','SUSCEPTIBLES')
 +
hold off
 +
subplot(2,2,2)
 +
hold on
 +
plot(t2,F);
 +
plot(t2,R,'g');
 +
legend('INFECTADOS','SUSCEPTIBLES')
 +
hold off
 +
subplot(2,2,3)
 +
hold on
 +
plot(t3,A);
 +
plot(t3,L,'y');
 +
legend('INFECTADOS','SUSCEPTIBLES')
 +
hold off
 +
subplot(2,2,4)
 +
hold on
 +
plot(t4,G);
 +
plot(t4,U,'c');
 +
legend('INFECTADOS','SUSCEPTIBLES')
 +
hold off }}
 +
 
 +
Éste programa de MatLab nos dará el valor máximo de población infectada en cada gráfica que son:
 +
1.I=1.26+e004
 +
2.I=9521.1
 +
3.I=9469.6
 +
4.I=9464.7
 +
Observamos en este segundo programa que las gráficas y resultados no parecen concordar con lo que sería la realidad
 +
ésto puede ser debido a que la alta diferencia entre los valores iniciales de I y S no permite elcálculo efectivo de
 +
su variación respecto del tiempo y por tanto la infección no podrá ser erradicada.
 +
 
 +
 
 +
== 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
 +
[[Archivo:Graficazorra3.png|400px|thumb|right|Comparación del método Runge-Kutta y de Euler]]
 +
  {{matlab|codigo=
 +
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');
 +
}}
 +
 
 +
== 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>. La interpretación más lógica para que la proporción de gente susceptible que se relaciona con los infectados vaya disminiyendo con el tiempo es que, debido al conocimiento de la epidemia, la gente empiece a tomar medidas para no estar en contacto con los infectados.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 infectado.
 +
 
 +
Para resolver este sistema hacemos uso del metodo de Heun, ayudandonos del siguiente código :
 +
[[Archivo:ejercicio69.png|400px|thumb|right|Método de Heun]]
 +
{{matlab|codigo=
 +
%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
 +
}}
 +
 
 +
Es destacable ver que con un valor para la a constante y con estos valores iniciales parecidos, el modelo no es válido, pero con estas consideraciones ya si nos queda un resultado que se ajusta más a la realidad.
 +
 
 +
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. Concluimos que el valor de a para el cual el máximo está más cerca de 5 es 0.0005, es decir, el menor del intervalo, y a medida que crece la I máxima se alcanza antes que t=5. A continuación se muestra el programa utilizado así como la tabla obtenida con este.
 +
[[Archivo:bonethijoputa.png|400px|thumb|right|Tabla de datos obtenidos]]
 +
 
 +
{{matlab|codigo=
 +
 
 +
clear all;
 +
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
 +
 
 +
}}
 +
 
 +
 
 +
[[Categoría:Ecuaciones Diferenciales]]
 +
 
 +
[[Categoría:ED14/15]]
 +
 
 +
[[Categoría:Trabajos 2014-15]]

Revisión actual del 20:19 3 dic 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] \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 Modelo Completo

En este caso obtendremos valores de S e I en función del tiempo por el método de Euler resolviendo el sistema debido a que el sistema no es lineal aplicamos Euler para cada una de las ecuaciones, lo que nos dará el resultado gráficamente de S e I y sus evoluciones en el tiempo, además utilizaremos cuatro longitudes de paso lo cual nos reportará cuatro gráficas en cada caso. Lo resolvemos para unos valores iniciales I(0)=20 S(0)=800 e I(0)=40 S(0)=10000. En el primer caso utilizamos el siguiente código MatLab:

Modelo Completo
% Euler 
clear all
%DATOS DEL PROBLEMA
t0=0;
tN=40;
I0=20;
S0=800;
h1=0.1;
h2=0.01;
h3=0.001;
h4=0.0001;
a=0.003;
b=0.3;
c=0.01;
%Calculamos número de subintervalos
N1=round((tN-t0)/h1); 
N2=round((tN-t0)/h2); 
N3=round((tN-t0)/h3); 
N4=round((tN-t0)/h4); 
% Definimos la variable independiente
t1=t0:h1:tN;
t2=t0:h2:tN;
t3=t0:h3:tN;
t4=t0:h4:tN;
%Ahora vamos a guardar los valores de la solución aproximada en el vector y
I=zeros(1,N1+1);   
S=zeros(1,N1+1);
F=zeros(1,N2+1);   
R=zeros(1,N2+1);
A=zeros(1,N3+1);   
L=zeros(1,N3+1);
G=zeros(1,N4+1);   
U=zeros(1,N4+1);
I(1)=I0;
S(1)=S0;
F(1)=I0;
R(1)=S0;
A(1)=I0;
L(1)=S0;
G(1)=I0;
U(1)=S0;
for i=1:N1
    S(i+1)=S(i)+h1.*(-a.*S(i).*I(i)); %eulerS
    I(i+1)=I(i)+h1.*(a.*S(i).*I(i)-b.*I(i)-c.*I(i)) ; %eulerI
 end
for i=1:N2
    R(i+1)=R(i)+h2.*(-a.*R(i).*F(i)); %eulerS
    F(i+1)=F(i)+h2.*(a.*R(i).*F(i)-b.*F(i)-c.*F(i)) ; %eulerI
 end
for i=1:N3
    L(i+1)=L(i)+h3.*(-a.*L(i).*A(i)); %eulerS
    A(i+1)=A(i)+h3.*(a.*L(i).*A(i)-b.*A(i)-c.*A(i)) ; %eulerI
 end
for i=1:N4
    U(i+1)=U(i)+h4.*(-a.*U(i).*G(i)); %eulerS
    G(i+1)=G(i)+h4.*(a.*U(i).*G(i)-b.*G(i)-c.*G(i)) ; %eulerI
 end
 J1=max(I);
 J2=max(F);
 J3=max(A);
 J4=max(G);
 disp('El valor maximo en la grafica 1 es')
 J1
  disp('El valor maximo en la grafica 2 es')
 J2
  disp('El valor maximo en la grafica 3 es')
 J3
  disp('El valor maximo en la grafica 4 es')
 J4
subplot(2,2,1)
hold on
plot(t1,I);
plot(t1,S,'r');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,2)
hold on
plot(t2,F);
plot(t2,R,'g');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,3)
hold on
plot(t3,A);
plot(t3,L,'y');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,4)
hold on
plot(t4,G);
plot(t4,U,'c');
legend('INFECTADOS','SUSCEPTIBLES')
hold off


Éste programa de MatLab nos dará además el valor máximo de población infectada en cada gráfica que son: 1.I=517.45 2.I=506.37 3.I=505.3 4.I=505.19 En el segundo caso utilizamos el código MatLab:

Modelo Completo
% Euler 
clear all
%DATOS DEL PROBLEMA
t0=0;
tN=40;
I0=40;
S0=10000;
h1=0.1;
h2=0.01;
h3=0.001;
h4=0.0001;
a=0.003;
b=0.3;
c=0.01;
%Calculamos número de subintervalos
N1=round((tN-t0)/h1); % también N=length(t)-1;
N2=round((tN-t0)/h2); % también N=length(t)-1;
N3=round((tN-t0)/h3); % también N=length(t)-1;
N4=round((tN-t0)/h4); % también N=length(t)-1;
% Definimos la variable independiente
t1=t0:h1:tN;
t2=t0:h2:tN;
t3=t0:h3:tN;
t4=t0:h4:tN;
%Ahora vamos a guardar los valores de la solución aproximada en el vector y
I=zeros(1,N1+1);   
S=zeros(1,N1+1);
F=zeros(1,N2+1);   
R=zeros(1,N2+1);
A=zeros(1,N3+1);   
L=zeros(1,N3+1);
G=zeros(1,N4+1);   
U=zeros(1,N4+1);
I(1)=I0;
S(1)=S0;
F(1)=I0;
R(1)=S0;
A(1)=I0;
L(1)=S0;
G(1)=I0;
U(1)=S0;
for i=1:N1
    S(i+1)=S(i)+h1.*(-a.*S(i).*I(i)); %eulerS
    I(i+1)=I(i)+h1.*(a.*S(i).*I(i)-b.*I(i)-c.*I(i)) ; %eulerI
 end
for i=1:N2
    R(i+1)=R(i)+h2.*(-a.*R(i).*F(i)); %eulerS
    F(i+1)=F(i)+h2.*(a.*R(i).*F(i)-b.*F(i)-c.*F(i)) ; %eulerI
 end
for i=1:N3
    L(i+1)=L(i)+h3.*(-a.*L(i).*A(i)); %eulerS
    A(i+1)=A(i)+h3.*(a.*L(i).*A(i)-b.*A(i)-c.*A(i)) ; %eulerI
 end
for i=1:N4
    U(i+1)=U(i)+h4.*(-a.*U(i).*G(i)); %eulerS
    G(i+1)=G(i)+h4.*(a.*U(i).*G(i)-b.*G(i)-c.*G(i)) ; %eulerI
 end
 J1=max(I);
 J2=max(F);
 J3=max(A);
 J4=max(G);
 disp('El valor maximo en la grafica 1 es')
 J1
  disp('El valor maximo en la grafica 2 es')
 J2
  disp('El valor maximo en la grafica 3 es')
 J3
  disp('El valor maximo en la grafica 4 es')
 J4
subplot(2,2,1)
hold on
plot(t1,I);
plot(t1,S,'r');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,2)
hold on
plot(t2,F);
plot(t2,R,'g');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,3)
hold on
plot(t3,A);
plot(t3,L,'y');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,4)
hold on
plot(t4,G);
plot(t4,U,'c');
legend('INFECTADOS','SUSCEPTIBLES')
hold off


Éste programa de MatLab nos dará el valor máximo de población infectada en cada gráfica que son: 1.I=1.26+e004 2.I=9521.1 3.I=9469.6 4.I=9464.7 Observamos en este segundo programa que las gráficas y resultados no parecen concordar con lo que sería la realidad ésto puede ser debido a que la alta diferencia entre los valores iniciales de I y S no permite elcálculo efectivo de su variación respecto del tiempo y por tanto la infección no podrá ser erradicada.


3 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');


4 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]. La interpretación más lógica para que la proporción de gente susceptible que se relaciona con los infectados vaya disminiyendo con el tiempo es que, debido al conocimiento de la epidemia, la gente empiece a tomar medidas para no estar en contacto con los infectados.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 infectado.

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


Es destacable ver que con un valor para la a constante y con estos valores iniciales parecidos, el modelo no es válido, pero con estas consideraciones ya si nos queda un resultado que se ajusta más a la realidad.

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. Concluimos que el valor de a para el cual el máximo está más cerca de 5 es 0.0005, es decir, el menor del intervalo, y a medida que crece la I máxima se alcanza antes que t=5. A continuación se muestra el programa utilizado así como la tabla obtenida con este.

Tabla de datos obtenidos
clear all;
 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
Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Modelos_Epidemológicos_Grupo_4-C&oldid=32987»