Modelo de Epidemia (G-23)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Modelo para epidemias. Grupo 23 |
| Asignatura | Ecuaciones Diferenciales |
| Curso | Curso 2014-15 |
| Autores | Álvaro Roales Blanco, Carlos Fernández Bermejo, Adrián Gómez Apiñániz, Marta Torra Escánez, Ricardo Mazón Cabrera |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
1 Introducción
El objeto de estudio de nuestro trabajo es el modelo de expansión de una epidemia a lo largo de un tiempo. Usaremos las siguientes hipótesis:
- La población (constante) se indicará como N en el que todos serán susceptibles a la enfermedad.
- La duración de la enfermedad es larga de modo que los individuos no se recuperarán durante el periodo de estudio.
- La unidad de tiempo escogida será la semana.
- Durante cada unidad de tiempo cada persona infectada tiene c contactos.
2 Determinación del número de contactos c
El primer paso a realizar es determinar el número de contactos que tiene una persona infectada en la unidad de tiempo elegida, es decir, la semana. Emplearemos la siguiente ecuación diferencial que nos muestra el enunciado del trabajo:
- [math]I'(t) = \frac{c}{N}I(t)(N-I(t))[/math]
c es un número de dos decimales comprendido entre [math] 0.01 [/math] y [math] 0.99 [/math]. Siendo N igual a 500000 con las siguientes condiciones iniciales: [math] I(0)=200 [/math] y [math] I(1)=500 [/math] Para determinar c imponemos la siguiente condición, que [math] |I(1)-500| [/math] tenga el menor valor posible, es decir, aquel valor que minimice lo posible el error absoluto. Para hallar ese valor lo hallaremos numéricamente con el siguiente programa escrito en MATLAB, basado en el método de Euler para la resolución de ecuaciones diferenciales:
y0 = 200; % Número inicial de infectados
y = [];
c = 1/100:0.01:99/100; % Posibles valores de c
N = 500000; % Población total
k = length(c);
Cminimo = c(1);
for i = 1:k
Yfin = N;
Y = N*y0/(y0+(N-y0).*exp(-c(i))); % Modelo para t=1
y = [y Y];
if abs(Y-500) <= abs(Yfin-500)
Cminimo = c(i);
end
end
min(abs(y-500)) % Error mínimo absoluto cometido
C_min
plot(c,abs(y-500),'r')Al ejecutar el programa anteriormente escrito obtenemos el valor c minimizando así el error que será igual a [math] 0.92 [/math] Por otra parte hemos dibujado un gráfico para ver la variación del error para cada valor de c:
3 Estudio de la evolución de la enfermedad
A continuación haremos el estudio de la evolución de la epidemia a lo largo de las semanas mediante los dos siguientes métodos.
Primero explicaremos el Método de Heun:
Consiste en la aproximación a la pendiente mediante la aplicación de dos derivadas del intervalo, una en el punto inicial y otra en punto final. La aproximación mejorada de la pendiente será el promedio de las dos derivadas, se calcula los siguientes elementos:
- [math] K_{1}=f(t_{n},y_{n}) [/math]
- [math] K_{2}=f(t_{n}+h,y_{n}+K_{1}h) [/math]
- [math] y_{n+1}=y_{n}+\frac{h}{2}(K_{1}+K_{2}) [/math]
Siendo h la longitud de paso entre dos valores consecutivos
En segundo lugar expondremos el Método de Runge-Kutta:
Sería similar al método de Heun pero se descompone en cuatro soluciones, se hace del siguiente modo:
- [math] K_{1}=f(t_{n},y_{n}) [/math]
- [math] K_{2}=f(t_{n}+\frac{h}{2},y_{n}+\frac{h}{2}K_{1}) [/math]
- [math] K_{3}=f(t_{n}+\frac{h}{2},y_{n}+\frac{h}{2}K_{2}) [/math]
- [math] K_{4}=f(t_{n}+h,y_{n}+K_{3}h) [/math]
- [math] y_{n+1}=y_{n}+\frac{h}{6}(K_{1}+2K_{2}+2K_{3}+K_{4}) [/math]
Mediante los métodos expuestos anteriormente podemos determinar el gráfico de evolución de los infectados conforme avanzan las semanas. Para ello aplicaremos el siguiente código MATLAB:
t0 = 0;h = 0.01;
r = [];
y = [];
y(1) = 200;
r(1) = 200;
i = 1;
c = 0.92;
N = 500000;
for i=1:N;
while y < 400000 & r < 400000
% Método 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);
% Método de Heun
K1 = (c*r(i))/N*(N-r(i));
K2 = (c*(r(i)+K1*h))/N*(N-(r(i)+K1*h));
r(i+1) = r(i) + (h/2)*(K1+K2);
i = i + 1;
end
end
tN1 = h*(length(y)-1);
tN2 = h*(length(r)-1);
t1 = t0:h:tN1;
t2 = t0:h:tN2;
Semana_maxima_infeccion = [t1(end) t2(end)]
Valor_aproximado = [y(8),r(8)]
hold on
plot(t1,y,'k')
plot(t2,r,'g')
legend('Runge Kutta','Heun','Location','best')
hold off
Como podemos observar en el gráfico adjunto, los infectados no tienden a desaparecer sino a crecer exponencialmente. Como es una función creciente no posee máximo, sin embargo como el número de habitantes está acotados el máximo se produce cuando todos los habitantes están infectados, es decir 500.000. Este fenómeno se produce pasada la semana número 10.
Al finalizar la sexta semana el número aproximado de infectados es 110.000.
A partir de la décima semana, el número de infectados es mayor a 400.000.
4 Susceptibilidad al contagio
t0 = 0;% Tiempo inicial
h = 0.01;% Paso del método numérico
tN = 10; % Última semana
N = 500000; % Tamaño de población
M = 10000; % Umbral de contagio
c = 0.92; % Susceptibilidad al contagio
t = t0:h:tN; % Valores del tiempo
n = round((tN-t0)/h); % Número de valores del tiempo(instantes de tiempo)
y = zeros(1,n+1); % Valores de personas contagiadas
y1(1) = 9700;
y2(1)=10200;
y3(1)=30000;
for i = 1:n
% Método de Heun
K11 = (c*y1(i))/N*(N-y1(i))*((y1(i)/M)-1);
K21 = (c*(y1(i)+K11*h))/N*(N-(y1(i)+K11*h))*(((y1(i)+K11*h)/M)-1);
y1(i+1) = y1(i) + (h/2)*(K11+K21);
K12 = (c*y2(i))/N*(N-y2(i))*((y2(i)/M)-1);
K22 = (c*(y2(i)+K12*h))/N*(N-(y2(i)+K12*h))*(((y2(i)+K12*h)/M)-1);
y2(i+1) = y2(i) + (h/2)*(K12+K22);
K13 = (c*y3(i))/N*(N-y3(i))*((y3(i)/M)-1);
K23 = (c*(y3(i)+K13*h))/N*(N-(y3(i)+K13*h))*(((y3(i)+K13*h)/M)-1);
y3(i+1) = y3(i) + (h/2)*(K13+K23);
end
figure(1)
axis equal
plot(t,y1)
figure(2)
axis equal
plot(t,y2)
figure(3)
axis equal
plot(t,y3)
Como podemos observar en la gráfica, los infectados desaparecen a la décima semana siendo el máximo el valor alcanzado en la primera semana debido a que la función es estrictamente decreciente.
En estas dos últimas gráficas la solución es creciente por lo que el número de infectados nunca es cero, alcanzando el máximo en 500.000 (valor de la población)
La diferencia entre estos dos últimos modelos es que el máximo se alcanza en la cuarta semana para el segundo modelo y antes de concluir la primera semana para el último modelo. Esto es coherente debido a que cuanto mayor sea el número de infectados iniciales antes se contagiará al resto ya que hay menos gente sana.
5 Modelo SIR
Definimos a continuación los siguientes parámetros, los cuales fijarán el número y las variaciones de infectados.
-S: Número de personas susceptibles a la enfermedad.
-I: Número de infectados
-R: Número de personas en cuarentena, muertos o que han superado o son inmunes a la enfermedad.
Definimos S(t), I(t), R(t), cantidad de miembros de cada clase al transcurrir un tiempo t.
Posibles suposiciones acerca de los procesos de incubación y transmisión (posible cambiar) considerando el tiempo de infección despreciable.
-rSI, r>0: variacíon del número de infectados y susceptibles. Se define también r como Tasa de infección.
-aI, a>0: existe una relación proporcional entre el número de infectados y la tasa de miembros infectados que pasan a la clase R.
-a: se define como la tasa de remoción.
Partimos de la suposición inicial considerandio S(0)=So, I(0)=Io y R(0)=0; siendo N=So + Io.
El número de la N (población total) no varía puesto que es constante por lo que podemos afirmar que la suma de S(t),I(t),R(t) se corresponde con N y esto con el número inicial de población susceptible e infectada, puesto que R(0)=0
Si[math]S_{0}\le\frac{a}{r}[/math] entonces [math]I(t)\le I_{0}[/math], siendo [math]\frac{a}{r}[/math] la tasa de remoción. Sustituyendo en I'(t) se demuestra que
- [math]I'(0)\le rS_{0}I_{0} - aI_{0}\le r\frac{a}{r}I_{0} - aI_{0}\le aI_{0} - aI_{0}\le 0[/math]
Análogamente si [math]S_{0}\gtfrac{a}{r}[/math] el número de individuos enfermos aumentará
- [math]I'(0)\geq rS_{0}I_{0} - aI_{0}\geq r\frac{a}{r}I_{0} - aI_{0}\geq aI_{0} - aI_{0}\geq 0[/math] .
5.1 Gráficas e interpretación con distintos valores
El siguiente código Matlab, que utiliza el método Euler y Runge-Kutta, se ejecuta para distintos valores de So, Io y h. Se encuentran ordenados con el mismo So e Io pero variando la h, tomando los valores 0.1, 0.01, 0.001.
% Datos. Irán cambiando según el caso
t0 = 0;
a = 0.341;
tN = 3.5;
S0 = 140000;
I0 = 20000;
r = 0.0000218;
h = 0.1;
N = round((tN-t0)/h);
t = t0:h:tN;
S1 = zeros(1,N+1);
I1 = zeros(1,N+1);
R1 = zeros(1,N+1);
S2 = zeros(1,N+1);
I2 = zeros(1,N+1);
R2 = zeros(1,N+1);
S1(1) = S0;
I1(1) = I0;
R1(1) = 0;
S2(1) = S0;
I2(1) = I0;
R2(1) = 0;
for i= 1:N
S1(i+1) = S1(i) + h*(-r*S1(i)*I1(i));
I1(i+1) = I1(i) + h*(r*S1(i)*I1(i)-a*I1(i));
R1(i+1) = R1(i) + h*(a*I1(i));
k1S = -r*S2(i)*I2(i);
k1I = r*S2(i)*I2(i)-a*I2(i);
k1R = a*I2(i);
k2S = -r*(S2(i) + 1/2*k1S*h)*(I2(i) + 1/2*k1I*h);
k2I = r*(S2(i) + 1/2*k1S*h)*(I2(i) + 1/2*k1I*h)-a*(I2(i) + 1/2*k1I*h);
k2R = a*(I2(i) + 1/2*k1I*h);
k3S = -r*(S2(i) + 1/2*k2S*h)*(I2(i) + 1/2*k2I*h);
k3I = r*(S2(i) + 1/2*k2S*h)*(I2(i) + 1/2*k2I*h)-a*(I2(i) + 1/2*k2I*h);
k3R = a*(I2(i) + 1/2*k2I*h);
k4S = r*(S2(i) + k3S*h)*(I2(i) + k3I*h);
k4I = r*(S2(i) + k3S*h)*(I2(i) + k2I*h)-a*(I2(i) + 1/2*k2I*h);
k4R = a*(I2(i) + k3I*h);
S2(i+1) = S2(i) + (h/6)*(k1S+2*k2S+2*k3S+k4S);
I2(i+1) = I2(i) + (h/6)*(k1I+2*k2I+2*k3I+k4I);
R2(i+1) = R2(i) + (h/6)*(k1R+2*k2R+2*k3R+k4R);
end
E1=S1-I1;
v1=find(E1<0);
t1=v1(1)*h % Momento en el que se igualan los susceptibles y los infectados mediante euler
E2=S2-I2;
v2=find(E2<0);
t2=v2(1)*h % Momento en el que se igualan los susceptibles y los infectados mediante Runge-Kutta
figure(1)
hold on
plot(t,S1,'b','Linewidth',2)
plot(t,I1,'r','Linewidth',2)
plot(t,R1,'g','Linewidth',2)
legend('S´ Euler','I´ Euler','R´ Euler','Location','best')
hold off
figure(2)
hold on
plot(t,S2,'b','Linewidth',2)
plot(t,I2,'r','Linewidth',2)
plot(t,R2,'g','Linewidth',2)
legend('S´ Runge-Kutta','I´ Runge-Kutta','R´ Runge-Kutta','Location','best')
hold off
Según observamos en las gráficas, el modelo que mejor se ajusta a los parámetros dados es el que tiene los valores S0=159985, I0=15 siendo los tiempos en los que se igualan los infectados y los susceptibles, que coincide con el comienzo de un crecimiento exponencial de los removidos, t1=3,5 semanas según el método de Euler y t2=3,1 semanas para el método de Runge-Kutta para un paso de h=0.1, t1=3,03 semanas para Euler y t2=2,99 semanas para Runge-Kutta para un paso de h=0.01 y por último para h=0.001 t1=2,985 según el método de Euler, t2=2,972 para el método de Runge-Kutta observando así una mayor exactitud con un paso menor (h=0.001) debido a que el parámetro dado es 1/a (2.932).
Cabe destacar el siguiente punto de corte entre las curvas de removidos e infectados que nos indica que se iguala el número de removidos e infectados apreciando que el número de susceptibles es prácticamente nulo, marcando el final del contagio y la recuperación o muerte de los infectados.
5.2 Trayectorias de los infectados
El primer caso que hemos realizado ha sido hallar la constante k para los diferentes valores iniciales de susceptibles e infectados. Posteriormente resolvemos la ecuación diferencial
- [math]I(t)= -S(t) + \frac{a}{r}\log{S(t)} + k[/math]
que dependerá unicamente de S por lo que podemos dibujar los siguientes gráficos de I frente a S.
h=100;
a=0.341;
r=0.0000218;
S0=159985;
S=[158990:-h:0];
I0=15;
p=a/r;
k=I0-p*log(S(1))+S(1);
N=length(S);
I=zeros(1,N);
for i=1:N
I(i)=p*log(S(i))+k-S(i);
end
plot(S,I)
title('Solución analítica')
xlabel('Individuos susceptibles de contagio')
ylabel('Individuos contagiados')
5.3 Obtención numérica del número máximo de infectados
Derivando la ecuación diferencial del apartado anterior obtenemos lo siguiente:
[math]\[ I(máximo) = N - \frac{a}{r} + \frac{a}{r}\log\frac{a}{rS_{0}}\ltmath\gt Para realizar el último apartado hemos creado el siguiente código MATLAB que nos calculará el número máximo de infectados {{matlab|codigo= N=500000; a=0.341; r=0.0000218; z=a/r; s0=159985; I=N-z+z*log(z/s0) fprintf('El número máximo de infectados es %f',I) }} Dándonos en la ventana de comandos "El número máximo de infectados es 447987.99652" obteniendo el mismo resultado mediante el procedimiento analítico, por lo que la máxima cantidad de infectados es 447988 siendo en la semana 3.8420, hallada utilizando el código de Matlab del Modelo SIR mediante el método de Runge-Kutta con los valores So=159985, Io=15 y un paso h=0,001 y añadiendo el siguiente código hemos hallado la semana en la que se alcanza el número máximo de infectados: {{matlab|codigo= valormaximoinfec=max(I2); find(I2==w); % Encontramos la posición en la que se encuentra el número máximo de infectados semanamaxinfec=t(3843)%nos da la semana en la que el número de infectados es máximo }} [[Categoría:Ecuaciones Diferentiales]] [[Categoría:ED14/15]] [[Categoría:Trabajos 2014-15]][/math]
