Modelo para epidemias. Grupo 6-C
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Modelo para epidemias. Grupo 6-C |
| Asignatura | Ecuaciones Diferenciales |
| Curso | Curso 2015-16 |
| Autores |
Manuel Jesús García Vega |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
1 Introducción
En este artículo vamos a centrar nuestro estudio en el desarrollo de modelos que nos permitan hacer una estimación del comportamiento temporal de una enfermedad infecciosa sin extensión espacial.
Suponemos las siguientes hipótesis:
1. La población es un número fijo N y cada miembro de la población es susceptible a la enfermedad.
2. La duración de la enfermedad es larga, de manera que no se cura durante el periodo de estudio.
3. Todos los individuos infectados son contagiosos y circulan libremente entre la población.
4. Durante cada unidad de tiempo cada persona infectada tiene c contactos y cada contacto con una persona no infectada redunda en la transmisión de la enfermedad.
2 Cálculo del número de contactos aproximados 'c'
Para realizar el cálculo del número de contactos ‘aproximados’ c que tiene una persona por unidad de tiempo emplearemos un método numérico basado en la siguiente ecuación logística:
donde:
I(t): Número de personas contagidas al cabo de un tiempo de t semanas.
c: Número de contactos que tiene una persona pasadas t semanas.
N: Población total en el lugar de estudio.
Además sabemos que los servicios de salud pública registran la difusión de una epidemia de gripe de duración particularmente larga en una ciudad de 500 000 personas. Al inicio de la primera semana de registro se habían contabilizado 200 casos; durante la primera semana aparecieron 300 nuevos casos.
Sabiendo que 'c' tiene valores entre 0.01 y 0.99, hemos creado un bucle que pasa por cada uno de esos valores, guardando en un vector la condición de |I(primera semana) - 500|, se ha representado gráficamente en comparación con todos los valores posibles de 'c' y también hemos obtenido la solución de forma numérica de cómo varía según la condición.
Para obtener el número de contactos utilizaremos el siguiente método númerico, el cual está basado en el método de Euler.
clear, close all
%Datos del problema
y=[];
y0=200;
c=linspace(1/100,99/100,99);
N=500000; %Población total;
x=linspace(0,300,99);
n=length(c); %numero de puntos
c0=c(1);
Y=N
%Esquema númérico
for i=1:n
Yfinal=Y;
Y=N*y0/(y0+(N-y0).*exp(-c(i))); % en t=1
y=[y;Y]
if abs(Y-500)<=abs(Yfinal-500)
c0=c(i)
end
end
min(abs(y-500)) %Error absoluto mínimo
plot(c,abs(y-500))
Al ejecutar el programa obtenemos el siguiente gráfico que representa la evolución del error que se produce para cada valor de 'c' . Y se obtiene el valor de 'c' que minimiza el error el cual es igual a 0,92.
3 Evolución de los Infectados por la Epidemia
Para determinar la evolución de los infectados por la epidemia conforme avanza el tiempo, emplearemos el método de Hern y el método de Runge Kutta de orden 4.
Método de Heun:
Método de Runge-Kutta:
Mediante los métodos expuestos, el de Heun y Runge-Kutta, podemos determinar la evolución que experimentan los infectados a lo largo de las semanas. Dicha Evolución la vamos a representar gráficamente mediante el siguiente código empleado en MATLAB:
t0=0;
h=0.01;
z=[];
y=[];
y0=200;
z0=200;
y(1)=y0;
z(1)=z0;
c =0.92;
N=500000;
for i=1:N;
while y<400000 & z<400000
% Método de Heun
K1H=(c*z(i))/N*(N-z(i));
K2H=(c*(z(i)+K1H*h))/N*(N-(z(i)+K1H*h));
z(i+1)=z(i)+(h/2)*(K1H+K2H);
%Método de Runge Kutta
k1=(c*y(i))/N*(N-y(i));
k2=(c*(y(i)+1/2*k1*h))/N*(N-(y(i)+1/2*k1*h));
k3=(c*(y(i)+1/2*k2*h))/N*(N-(y(i)+1/2*k2*h));
k4=(c*(y(i)+k3*h))/N*(N-(y(i)+h*k3));
y(i+1)=y(i)+(h/6)*(k1+2*k2+2*k3+k4);
i=i+1;
end
end
%Discretización variable independiente
tN1=h*(length(y)-1);
tN2=h*(length(z)-1);
ty=t0:h:tN1;
tz=t0:h:tN2;
Semana_maxima_infeccion=[ty(end) tz(end)]
ValorAproximado=[y(6/h),z(6/h)]
hold on
plot(ty,y,'b')
plot(tz,z,'r')
xlabel('Tiempo en semananas');
ylabel('Número de infectados');
legend('Runge Kutta','Heun','Location','best')
hold off
Tal y cómo se puede apreciar en los gráficos adjuntos, el número de infectados es siempre creciente y lo hace de forma exponencial, de manera que el número de infectados seguirá aumentando hasta llegar a contagiar a toda la población, que en nuestro caso es de 500.000 habitantes. Al tratarse de una función exponencial que está acotada se comprueba la veracidad del máximo anteriormente dicho (500.000 habitantes).
También se aprecia que hay muy poca diferencia entre las gráficas obtenidas por Heun y Runge-Kutta, de hecho, para ver de forma más clara la diferencia entre éstas, hemos tenido que aumentar hasta el orden de la diezmilésima en el eje de las X y teniendo en cuenta que la unidad es la semana.
Además de las gráficas, el código nos devuelve el número de infectados pasadas 6 semanas y el tiempo que pasa hasta llegar a 400.000 personas infectadas.
-Utilizando el método de Heun:
A las 6 semanas hay 45.029 infectados. Se llega a los 400.000 infectados pasadas 10 semanas y un día y medio aproximadamente.
-Utilizando el método de Runge-Kutta:
A las 6 semanas hay 45.032 infectados. Se llega a los 400.000 infectados pasadas 10 semanas y un día y medio aproximadamente.
4 2a Parte: Modelo SIR
Ahora vamos a usar un modelo más general para estudiar la enfermedad, que divide la población total N en susceptibles de ser infectados(S), infectados(I), y removidos, que agrupa tanto a los curados como a los muertos (R). Así, el modelo, que tiene tres incógnitas, S(t), I(t) y R(t), y llamando r a la tasa de infección y a a la de remoción, lo podemos escribir:
4.1 Demostración de que S\left ( t \right )+I\left (t\right )+R\left ( t\right )=N
Sumando la variación de los susceptibles, los infectados y los removidos, vemos que esta es igual a 0, lo que nos quiere decir que la población total N es constante: \frac{\partial S}{\partial t}+\frac{\partial I}{\partial t}+\frac{\partial R}{\partial t}= S’\left ( t \right )+I’\left (t\right )+R’\left ( t\right )\left ( -rSI \right )+\left ( rSI-aI \right )+\left ( aI\right )=0, Además, si sumamos sus valores iniciales vemos que estos son igual a la población total, y al ser ésta constante, vemos que se demuestra el enunciado inicial S\left ( 0 \right )+I\left (0\right )+R\left ( 0\right )=N
4.2 Demostraciones analíticas con la tasa de remoción relativa
Llamamos tasa de remoción relativa a: \rho =\frac{r}{a} Vamos a demostrar que si dicha tasa es mayor o igual que el número de personas susceptibles iniciales, S\left ( 0 \right ), no hay epidemia, ya que el número de infectados no aumenta: También vamos a demostrar que, si dicha tasa es menor que la población susceptible inicial, habrá epidemia, pues aumentarán el número de infectados:
4.3 Resolución numérica
Ahora vamos a resolver numéricamente el problema de epidemias, usando los métodos de Euler, Euler implícito y Runge-Kutta de orden 4. Los aplicaremos a 3 casos, cada uno con diferentes datos, y probando con diferentes anchuras de paso: Mostramos el código para el primer caso con h=0.1
t0=0;%Intervalo de tiempo y anchura de paso
;tf=4;
h= 0.1
S0=159985; I0=15; % datos de la población
r=0.0000218; a=0.341; R0=0;
t=t0:h:tf;
k=((tf-t0)/h);
Se=[]; % método de euler
Ie=[];
Re=[];
Si=[]; % euler implícito
Ii=[];
Ri=[];
Sr=[]; % runge-kutta
Ir=[];
Rr=[];
Se(1)=S0;
Ie(1)=I0;
Re(1)=R0;
Si(1)=S0;
Ii(1)=I0;
Ri(1)=R0;
Sr(1)=S0;
Ir(1)=I0;
Rr(1)=R0;
% esquema numérico euler
for i=1:k
Ie(i+1)=Ie(i)+h*(r*Se(i)*Ie(i)-a*Ie(i));
Se(i+1)=Se(i)+h*(-r*Se(i)*Ie(i));
Re(i+1)=Re(i)+h*(a*Ie(i));
end
% esquema numérico euler implícito (es necesario despejar la I, la S y la R de (j+1) para definir el código
for j=1:k
Ii(j+1)=Ii(j)/(1-(h*r*Si(j))+a*h);
Si(j+1)=Si(j)/(1+r*h*Ii(j));
Ri(j+1)=Ri(j)+h*(a*Ii(j));
end
%esquema numerico runge-kutta
for n=1:k
k1S = -r*Sr(n)*Ir(n);
k1I = r*Sr(n)*Ir(n)-a*Ir(n);
k1R = a*Ir(n);
k2S = -r*(Sr(n) + 1/2*k1S*h)*(Ir(n) + 1/2*k1I*h);
k2I = r*(Sr(n) + 1/2*k1S*h)*(Ir(n) + 1/2*k1I*h)-a*(Ir(n) + 1/2*k1I*h);
k2R = a*(Ir(n) + 1/2*k1I*h);
k3S = -r*(Sr(n) + 1/2*k2S*h)*(Ir(n) + 1/2*k2I*h);
k3I = r*(Sr(n) + 1/2*k2S*h)*(Ir(n) + 1/2*k2I*h)-a*(Ir(n) + 1/2*k2I*h);
k3R = a*(Ir(n) + 1/2*k2I*h);
k4S = r*(Sr(n) + k3S*h)*(Ir(n) + k3I*h);
k4I = r*(Sr(n) + k3S*h)*(Ir(n) + k2I*h)-a*(Ir(n) + 1/2*k2I*h);
k4R = a*(Ir(n) + k3I*h);
Sr(n+1) = Sr(n) + (h/6)*(k1S+2*k2S+2*k3S+k4S);
Ir(n+1) = Ir(n) + (h/6)*(k1I+2*k2I+2*k3I+k4I);
Rr(n+1) = Rr(n) + (h/6)*(k1R+2*k2R+2*k3R+k4R);
end
hold on
plot(t,Ia,'--g',t,Sa,'--m',t,Ra,'--b')
plot(t,Ib,':g',t,Sb,':m',t,Rb,':b')
plot(t,Ic,'-g',t,Sc,'-m',t,Rc,'-b')
legend('I Euler','S Euler','R Euler','I Euler implícito','S Euler implícito','R Euler implícito','I Runge-Kutta','S Runge-Kutta','R Runge-Kutta','Location','best')
hold offA continuacion mostramos las gráficas obtenidas, para todas r=0,0000218 a=0,341:
| Archivo:Graficoo2.8.jpg Caso 2:S(0)=140000, I(0)=20000, h=0,1, 20 semanas |
Como se puede observar en las diferentes gráficas, cuanto menor es el tamaño de paso mas aproximado es el método, siéndolo el que más el de Runge-Kutta que además es el único de orden 4.
