Modelos epidemiológicos Grupo 6C
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Modelos epidemiológicos. Grupo 6-C |
| Asignatura | Ecuaciones Diferenciales |
| Curso | Curso 2014-15 |
| Autores |
Manuel Morales López 1175 David Toledo Menéndez 1228 Sergio Rodríguez Torcal 994 Jose María Rodríguez Vicente 1213 Lourdes Sánchez-Ocaña Merino 1248 Jorge Villa Lobo 1237 |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
1 Enunciado
En el desarrollo de una epidemia se distinguen dos tipos de individuos: los que ya han contraido la enfermedad o infectados I, y los que son susceptibles de contraerla por encontrarse en zona de riesgo S. Supongamos que se dan las siguientes hip´otesis: 1. La poblaci´on de personas infectadas se altera por el fallecimiento o la cura de las mismas. En ambos casos, la tasa de cambio depende del n´umero 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´on entre el n´umero de individuos en ambas clases. Consideramos las variables: t tiempo, S(t) poblaci´on de individuos susceptibles a contraer la enfermedad, I(t) poblaci´on de individuos infectados; y el sistema: dS dt = −aSI dI dt = aSI − bI − cI donde a, b, c son parametros.
2 Interpretación de parámetros
En el problema: el coeficiente "a" es la tasa de infectados por contagio, "b" la de muertos y "c" la de curados.
3 Método de Euler y Trapecio con S=0
%Apartado 2 trabajo1 2015
clear all
t0=0;
tN=20;
y0=2000;
h=0.1;
t=t0:h:tN;
N=(tN-t0)/h;
y=zeros(1,N+1); %euler
y(1)=y0;
z=zeros(1,N+1); %trapecio
z(1)=y0;
w=zeros(1,N+1);
w(1)=y0;
x=zeros(1,N+1);
x(1)=y0;
for i=1:N
y(i+1)=y(i)-h*(0.3+0.001)*y(i);
z(i+1)=(1/(1-h/2*(-0.31)))*(z(i))+h/2*(-0.31*z(i));
end
[t',y',z'];
i=1;
while abs(y0-w(i))<3/4*y0
w(i+1)=(1/(1-h/2*(-0.31)))*(w(i))+h/2*(-0.31*w(i)); %trapecio, evitamos la variable auxiliar yy
t(i+1)=t(i)+h;
i=i+1;
end
disp('Tiempo final con trapecio:')
disp(t(i))
%gráfico
i=1;
while abs(y0-x(i))<3/4*y0
x(i+1)=x(i)-h*(0.3+0.001)*x(i); %euler, evitamos la variable auxiliar yy
t(i+1)=t(i)+h;
i=i+1;
end
disp('Tiempo final con Euler:')
disp(t(i))
hold on
plot(t,y,'b')
plot(t,z,'r')
plot(t,w,'g')
plot(t,x,'k')
legend('Euler','Trapecio','Tiempo que tarda por Trapecio','Tiempo que tarda por Euler','Location','best');
hold offEsta gráfica nos muestra que en t=4.60 se alcanza la condición final, que el número de infectados se reduzca a 500 mientras que con el trapecio el tiempo final es t=4.50
4 Método de Euler y Trapecio con S=100
Se puede interpretar como que a partir de un valor limite de "S" entre 100 y 200, el numero de infectados se mantiene constante en el tiempo. De igual forma, si "S" es menor que este valor limite, el numero de infectados desciende en el tiempo, y si es mayor que el valor limite "S", asciende.
%Apartado 3 trabajo 6-C 2015
clear all
t0=0; tN=40; %Intervalo de tiempo tomado
S0=input('introduce valor de población susceptible: ');
I0=2000; %Valores iniciales de infectados(I) y susceptibles(S)
h=0.1; %Determinación del paso
t=t0:h:tN; %Desarrollo del tiempo desde t0 hasta tN tomando intervalos de paso h
a=0.003; %valores de los parametros
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0; %Asignacion del valor incial para la primera componente
I(1)=I0;
%Resolucion del sistema de forma matricial
for n=1:N
A=[S0;I(n)]+h*[-a*S0*I(n);a*S0*I(n)-(b+c)*I(n)];
S(n+1)=S0;%Asignacion de los distintos valores de S como la primera compnente de la matriz A
I(n+1)=A(2); %Asignacion de los distintos valores de I como la primera componente de la matriz A
end
%Interpretación gráfica
hold on
figure(1)
plot(t,S,'r')
plot(t,I,'b')
hold off
Podemos llegar a la conclusión de que el número de personas infectadas varía notablemente cuando el número de susceptibles se mantiene constante a lo largo del tiempo. Esto se debe a que cuanto mayor sea el número de susceptibles, mayor será el número de infectados, es decir, con 200 susceptibles, mantiene un ritmo constante hasta que llega un momento en el cual el número de infectados se dispara, debido a que la tasa de personas susceptibles es mayor que las personas curadas o fallecidas. En conclusión, debe de haber un valor de personas susceptibles para que la función de infectados sea lineal, es decir, sea constante en el tiempo.
5 Modelo completo
%Apartado 4 trabajo 6-C 2015
clear all
t0=0; tN=40; %Intervalo de tiempo tomado
S0=input('introduce valor de población susceptible: ');
I0=input('introduce valor de población infectada inicial: '); %Valores iniciales de infectados(I) y susceptibles(S)
h=input('introduce valor del paso de tiempo: '); %Determinación del paso
t=t0:h:tN; %Desarrollo del tiempo desde t0 hasta tN tomando intervalos de paso h
a=0.003; %valores de los parametros
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0; %Asignacion del valor incial para la primera componente
I(1)=I0;
%Resolucion del sistema de forma matricial
for n=1:N
A=[S(n);I(n)]+h*[-a*S(n)*I(n);a*S(n)*I(n)-(b+c)*I(n)];
S(n+1)=A(1);%Asignacion de los distintos valores de S como la primera compnente de la matriz A
I(n+1)=A(2); %Asignacion de los distintos valores de I como la primera componente de la matriz A
end
%Interpretación gráfica
hold on
figure(1)
plot(t,S,'r')
plot(t,I,'b')
hold off
maximo=max(I)
posicion=find(I==maximo)
tiempo=t(posicion)
Como podemos apreciar en las gráficas el número máximo de enfermos esperados con 800 susceptibles y 20 infectados inicialmente, es aproximadamente 500, mientras que para 10000 susceptibles y 40 infectados inicialmente, observamos un aumento de personas enfermas, aproximadamente de 9500.
6 Comparación Runge-Kutta con Euler
%Apartado 5 trabajo 6-C 2015
clear all
t0=0; tN=30; %Intervalo de tiempo tomado
y1=input('introduce valor de población susceptible inicial: ');
y2=input('introduce valor de población infectada inicial: ');
y0=[y1;y2]; %Valores de población iniciales (Susceptibles e Infectados)
h=input('introduce valor del paso del tiempo: ');
N=(tN-t0)/h;
y=y0;
S(1)=y(1); %Asignación del valor incial para la primera componente de S
I(1)=y(2); %Asignación del valor incial para la primera componente de I
a=0.003; % Valores de los parámetros
b=0.3; c=0.01; %Resolución empleando el método Runge Kutta
for n=1:N
k1=[-a*y(1)*y(2);a*y(1)*y(2)-(b+c)*y(2)];
k2=[-a*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2));a*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2)-(b+c)*(y(2)+1/2*h*k1(2)))];
k3=[-a*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2));a*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2)-(b+c)*(y(2)+1/2*h*k2(2)))];
k4=[-a*(y(1)+h*k3(1))*(y(2)+h*k3(2));a*(y(1)+h*k3(1))*(y(2)+h*k3(2)-(b+c)*(y(2)+h*k3(2)))];
y=y+h/6*(k1+2*k2+2*k3+k4);
S(n+1)=y(1);
I(n+1)=y(2);
end
t=t0:h:tN;
figure(1)
hold on
plot(t,S,'b')
plot(t,I,'r')
hold off
for n=1:N
A=[S(n);I(n)]+h*[-a*S(n)*I(n);a*S(n)*I(n)-(b+c)*I(n)];
S(n+1)=A(1);%Asignacion de los distintos valores de S como la primera compnente de la matriz A
I(n+1)=A(2); %Asignacion de los distintos valores de I como la primera componente de la matriz A
end
%Interpretación gráfica
figure(2)
hold on
plot(t,S,'b')
plot(t,I,'r')
hold off
La principal dificultad a la hora de emplear un método implícito radica en que estos métodos, al contrario que los explícitos no emplean la información del punto anterior para calcular el siguiente y requieren por tanto algoritmos mas complejos. Sin embargo presentan mejores aproximaciones que los métodos explícitos.
7 Tasa de infectados por contagio dependiente del tiempo. Método de Heun
clear all
%Datos
t0 = 0;
tN = 40;
h = 0.1;
N = round((tN-t0)/h);
t = t0:h:tN;
S = zeros(1,N+1);
I = zeros(1,N+1);
S(1) = 1600;
I(1) = 40;
%Valores de los parámetros
a = zeros(1,N+1);
a(1) = 0.003/(1+t0);
b = 0.3;
c = 0.01;
%Calcular los vectores infectados y susceptibles.
for i = 1:N
a(i+1) = 0.003/(1+t(i));
K1 = -a(i)*S(i)*I(i);
K2 = -a(i)*S(i)*I(i) + K1*h;
S(i+1) = S(i) + h/2*(K1+K2);
Z1 = (a(i)*S(i)*I(i)) -I(i)*(b+c);
Z2 = (a(i)*S(i)*I(i)) -I(i)*(b+c) + Z1*h;
I(i+1) = I(i) + h/2*(Z1+Z2);
end
%Gráficas
hold on
plot (t,S,'-')
plot (t,I,'-r')
hold off
Con el paso del tiempo el valor de la variable "a" disminuye considerablemente, puesto que en este apartado la "a" varia segun el tiempo, por esta razon, se ve amortiguada la caida de gente susceptible en la funcion.
8 Calibración del parámetro a
clear all
%Datos
t0 = 0;
tN = 8;
h = 0.0001;
N = (tN-t0)/h;
t = t0:h:tN;
S = zeros(1,N+1);
I = zeros(1,N+1);
S(1) = 1600;
I(1) = 40;
%Valores de los parámetros
a = 0.0005:h:0.0020;
b = 0.3;
c = 0.01;
%Definimos los vectores M, T y D, donde guardaremos respectivamente
%los máximos de infectados, el tiempo en el que se producen y la
%diferencia entre ese tiempo y los 5 días
M=zeros(1,length(a));
T=zeros(1,length(a));
D=zeros(1,length(a));
%Calcular los vectores infectados y susceptibles
for j=1:length(a)
A=a(j);
for i = 1:N
K1 = -A*S(i)*I(i);
K2 = -A*S(i)*I(i) + K1*h;
S(i+1) = S(i) + h/2*(K1+K2);
Z1 = (A*S(i)*I(i)) -I(i)*(b+c);
Z2 = (A*S(i)*I(i)) -I(i)*(b+c) + Z1*h;
I(i+1) = I(i) + h/2*(Z1+Z2);
end
M(j)=max(I)
posicion=find(I==M(j));
tiempo=t(posicion);
T(j)=tiempo
D(j)=abs(T(j)-5);
end
[M']
[D']
[T']
T1=min(D);
posicion=find(D==T1);
%soluciones obtenidas
A_buscado=a(posicion)
Maximo_en_ese_caso=M(posicion)
Cuando calibramos el coeficiente de a con una experiencia anterior en la que el máximo de personas infectadas se alcanza a los 5 días obtenemos el resultado de que a=0.0008 es el valor que mas se ajusta a esa experiencia, tomando como paso h=10^-4.
