Modelo para epidemias (Grupo 9C, Trabajo 7)
En este trabajo vamos a estudiar el comportamiento de una enfermedad infecciosa (epidemia) en una población de N habitantes, la cual es inextensible espacialmente (no varía el número de habitantes).
Contenido
- 1 ª Parte
- 2 ª Parte
- 2.1 Demostrar que S(t) + I(t) + R(t) = S0 + I0 = N
- 2.2 Demostración analítica de la epidemia en función de S0 y a/r
- 2.3 Resolucióin numérica del sistema SIR
- 2.4 Interpretación de la epidemia
- 2.5 Trayectorias de los casos del apartado 2.3
- 2.6 Número máximo de infectados en el primer caso del apartado 3
1 ª Parte
En la primera parte del trabajo vamos a emplear un modelo sencillo que implica una ecuación diferencial de primer orden. Como datos tenemos el número total de habitantes (N=500.000 habitantes), el número inicial de infectados (I(0)=200 casos), el número de infectados tras la primera semana (I(1)=500 casos).
1.1 Número de contactos 'aproximados'
Para calcular los contactos 'aproximados', hemos empleado el método de Euler. Sabiendo que C tiene valores entre 0.01 y 0.99, hemos creado un bucle que ha pasado por cada uno de esos valores, guardando en un vector la condición de |I(primera semana) - 500| y dibujándolo en comparación con todas las C posibles para ver gráficamente, además de numéricamente, la variación de la condición. El problema que tenemos que resolver es:
Y el código empleado para estos cálculos es el siguiente:
clear all, clc, clf
% Parte 1 Apartado 1----------------------------------------------------------------------------------------------
y0=200;
h=0.01;
t=0:0.01:1; % Vector en el que se estudia la epidemia
n=length(t); % Número de puntos
c=0.01:0.01:0.99; % Vector de posibles C
m=length(c); % Número de posibles C
Y=zeros(m,n); % Matriz donde guardamos, por filas, las soluciones de los infectados según C
y(1)=y0;
% Método de Euler-------------------------------------------------------------------------------------------------
for k=1:m
Y(k,1)=y0;
for i=1:n-1
y(i+1)=y(i)+h*(c(k)/500000*(500000-y(i))*y(i));
Y(k,i+1)=y(i+1);
end
end
I=abs(Y(:,n)-500*(ones(1,m)')); % Condición de abs(I(primera semana)-500)
minimo=min(I); % Mínimo de la condición
[u,v]=find(I==minimo); % Posición del mínimo
C=u/100 % C obtenida
plot(c,I) % Gráfica comparativa de C con la condiciónY la gráfica obtenida de la comparación de la condición en función de cada valor posible de C es:
1.2 Evolución de los infectados
Utilizando el método de Heun y Runge Kutta de orden 4, con un paso de h=0.01, la gráfica de evolución de afectados frente al tiempo (I(t),t) en un tiempo t=20 semanas, que hemos escogido para ver la curva completa, y el programa son los siguientes:
clc, clear all, clf
% Parte 1 Apartado 2-----------------------------------------------
t0=0;
tN=input('Introduce las semanas de estudio: ');
y0=200;
h=0.01;
N=round(tN-t0)/h;
% Vectores---------------------------------------------------------
t=t0:h:tN; % Tiempo
y=zeros(1,N+1); % Infectados para Runge Kutta
y(1)=y0;
z=zeros(1,N+1); % Infectados para Heun
z(1)=y0;
c=0.92; % Contactos
% Bucle------------------------------------------------------------
for i=1:N
K1RK=c/500000*(500000-y(i))*y(i);
K2RK=c/500000*(500000-(y(i)+h/2*K1RK))*(y(i)+h/2*K1RK);
K3RK=c/500000*(500000-(y(i)+h/2*K2RK))*(y(i)+h/2*K2RK);
K4RK=c/500000*(500000-(y(i)+h*K3RK))*(y(i)+h*K3RK);
y(i+1)=y(i)+h/6*(K1RK+2*K2RK+2*K3RK+K4RK);
K1H=c/500000*(500000-z(i))*z(i);
K2H=c/500000*(500000-(z(i)+h*K1H))*(z(i)+h*K1H);
z(i+1)=z(i)+h/2*(K1H+K2H);
end
% Gráficas de RK4 y Heun-------------------------------------------
figure(1)
plot(t,y)
legend('Runge Kutta orden 4')
figure(2)
plot(t,z,'r')
legend('Heun')
1.3 Infectados al terminar la sexta semana
Tras la sexta semana, los infectados se muestran con el siguiente programa, siendo: I(6)=4.5032e+004.
clc, clear all, clf
% Parte 1 Apartado 3-------------------------------------
t0=0;
tN=6;
y0=200;
h=0.01;
t=t0:h:tN; % Vector tiempo
N=length(t);
y=zeros(1,N); % Vector de infectados
y(1)=y0;
c=0.92; % Contactos
% Bucle de Euler-----------------------------------------
for i=1:N-1
K1=c/500000*(500000-y(i))*y(i);
K2=c/500000*(500000-(y(i)+h/2*K1))*(y(i)+h/2*K1);
K3=c/500000*(500000-(y(i)+h/2*K2))*(y(i)+h/2*K2);
K4=c/500000*(500000-(y(i)+h*K3))*(y(i)+h*K3);
y(i+1)=y(i)+h/6*(K1+2*K2+2*K3+K4);
end
InfectadosSemana6=y(600) % El I(6) es el punto 6/h=600
1.4 Semana a partir de la cuál hay más de 400 000 infectados
Estudiando las 20 primeras semanas de la epidemia, la semana a partir de la cual los infectados superan la cifra de 400 000 es la semana 10 (con mayor exactitud, a partir de las 10.3 primeras semanas). El código empleado para su cálculo ha sido el siguiente:
clc, clear all, clf
% Parte 1 Apartado 4--------------------------------------------
t0=0;
tN=input('Introduce las semanas de estudio: ');
y0=200;
h=input('Introduce el paso:');
t=t0:h:tN; % Vector tiempo
N=length(t);
y=zeros(1,N); % Vector de infectados
y(1)=y0;
c=0.92; % Contactos
% Bucle de Euler------------------------------------------------
for i=1:N-1
K1=c/500000*(500000-y(i))*y(i);
K2=c/500000*(500000-(y(i)+h/2*K1))*(y(i)+h/2*K1);
K3=c/500000*(500000-(y(i)+h/2*K2))*(y(i)+h/2*K2);
K4=c/500000*(500000-(y(i)+h*K3))*(y(i)+h*K3);
y(i+1)=y(i)+h/6*(K1+2*K2+2*K3+K4);
if y(i)<=400000
Sem=i+1;
end
end
plot(t,y)
Semana=Sem*0.01 % La posición de la semana será (i+1)*h
1.5 Infectados con dos umbrales
Manteniendo la población en 500 000 habitantes, utilizando el método de Heun, una longitud de paso h = 0.01 y un umbral M = 10 000, el código y las gráficas de evolución de infectados para los infectados iniciales 9 700, 10 200 y 30 000 habitantes son las siguientes:
clear all, clc, clf
% Parte 1 Apartado 5--------------------------------------------------------
t0=0;
tN=input('Introduce el numero de semanas: '); % Semanas
y0=input('Inserte valor de Infectados iniciales: '); % Infectados iniciales
h=0.01; % Paso
P=500000; % Población
M=10000; % Umbral
t=t0:h:tN; % Tiempo
N=length(t);
c=0.92; % Contactos
y=zeros(1,N);
y(1)=y0;
% Método de Heun------------------------------------------------------------
for i=1:N-1
K1=c*y(i)*(1-(y(i)/P))*(y(i)/M-1);
K2=c*(y(i)+K1*h)*(1-1/P*(y(i)+K1*h))*(1/M*(y(i)+K1*h)-1);
y(i+1)=y(i)+h/2*(K1+K2);
end
% Gráfica-------------------------------------------------------------------
plot(t,y)
Se puede ver que tenderán a desaparecer los infectados únicamente en el primero de los casos para un número inicial de infectados de 9700, esto es debido a que no se supera el umbral de 10000 (M) infectados. Cuando se supera este umbral M el número de infectados tiende a N=500000 en el infinito.
2 ª Parte
Vamos a considerar ahora un método más general, empleando un sistema de ecuaciones que es el siguiente:
Llamado método SIR, donde S(t) es el número de personas susceptibles de infectarse, I(t) es el número de infectados y R(t) es el número de recuperados e inmunes (al curarse de la enfermedad no se pueden volver a contagiar) y de habitantes aislados, y S(0), I(0) y R(0) son los datos iniciales de esos grupos.
2.1 Demostrar que S(t) + I(t) + R(t) = S0 + I0 = N
Si sumamos S' + I' + R' definidas en el sistema, tenemos que S' + I' + R' = (-rSI) + (rSI - aI) + (aI) = 0, lo que indica que la suma de S + I + R es constante y además S(t) + I(t) + R(t) = S(0) + I(0) + R(0) = S0 + I0 = N, es decir, la población total.
2.2 Demostración analítica de la epidemia en función de S0 y a/r
Demostración analítica, caso (1)y (2):
(1) En el primero de los caso si el número de susceptibles iniciales es menor o igual que la tasa de remoción relativa, no habrá epidemia, el número de infectados no aumenta y tienden a desaparecer en poco tiempo.
(2) En el segundo de los casos si el número de susceptibles iniciales es mayor que la tasa de remoción relativa, entonces hay epidemia en el sentido que los infectados empiezan a aumentar.
2.3 Resolucióin numérica del sistema SIR
Para resolver numéricamente el sistema utilizamos el método de Euler y Runge-Kutta de orden 4 para los distintos valores iniciales de [S0,I0] ([159 985,15],[140 000,20 000],[159 999,1]) y para los distintos intervalos de tiempo ([0,4] y [0,20]), con r=0.0000218 y a=0.341. El código empleado para su cálculo ha sido el siguiente:
clear all;clf;clc
% Introducimos las condiciones iniciales
t0=0;
tN=input('introduce la semana: ');
h=input('introduce el paso: ');
t=t0:h:tN;
S0=input('introduce los susceptibles iniciales: ');
R0=0;
I0=input('introduce el numero de infectados iniciales: ');
% Variables de RK4
I=zeros(1,length(t));
R=zeros(1,length(t));
S=zeros(1,length(t));
I(1)=I0;
S(1)=S0;
R(1)=R0;
% Variables de Euler
Ie=zeros(1,length(t));
Re=zeros(1,length(t));
Se=zeros(1,length(t));
Ie(1)=I0;
Se(1)=S0;
Re(1)=R0;
%datos
a=0.341;
r=0.0000218;
% Runge-Kutta de orden 4
for n=1:length(t)-1
%primer término
K1s=-r*S(n)*I(n);
K1i=r*S(n)*I(n)-a*I(n);
K1r=a*I(n);
%segundo término
K2s=-r*(S(n)+h/2*K1s)*(I(n)+h/2*K1i);
K2i=r*(S(n)+h/2*K1s)*(I(n)+h/2*K1i)-a*(I(n)+h/2*K1i);
K2r=a*(I(n)+h/2*K1i);
%tercer término
K3s=-r*(S(n)+h/2*K2s)*(I(n)+h/2*K2i);
K3i=r*(S(n)+h/2*K2s)*(I(n)+h/2*K2i)-a*(I(n)+h/2*K2i);
K3r=a*(I(n)+h/2*K2i);
%cuarto término
K4s=-r*(S(n)+h/2*K3s)*(I(n)+h/2*K3i);
K4i=r*(S(n)+h/2*K3s)*(I(n)+h/2*K3i)-a*(I(n)+h/2*K3i);
K4r=a*(I(n)+h/2*K3i);
S(n+1)=S(n)+h/6*(K1s+2*K2s+2*K3s+K4s);
I(n+1)=I(n)+h/6*(K1i+2*K2i+2*K3i+K4i);
R(n+1)=R(n)+h/6*(K1r+2*K2r+2*K3r+K4r);
end
% Euler
for i=1:length(t)-1
Se(i+1)=Se(i)+h*(Se(i)*(-0.0000218)*Ie(i)); %Euler primera ecuación
Ie(i+1)=Ie(i)+h*(0.0000218*Se(i)*Ie(i)-0.341*Ie(i)); %Euler segunda ecuación
Re(i+1)=Re(i)+h*(0.341*Ie(i)); %Euler tercera ecuación
end
hold on
plot(t,S,'g','linewidth',2,t,I,'b','linewidth',2,t,R,'r','linewidth',2)
plot(t,Se,'m','linewidth',2,t,Ie,'c','linewidth',2,t,Re,'k','linewidth',2)
legend('S con RK4','I con RK4', 'R con RK4','S con Euler','I con Euler', 'R con Euler','location','best')
hold off
CONDICIONES 1: S0=159 985, I0=15.
CONDICIONES 2: S0=140 000, I0=20 000.
CONDICIONES 3: S0=159 999, I0=1.
2.4 Interpretación de la epidemia
2.5 Trayectorias de los casos del apartado 2.3
Para los casos estudiados en
clear all, clf, clc
% Introducimos las condiciones iniciales----------------------------------------------------
t0=0;
tN=input('Introduce la semana: ');
h=input('Introduce el paso: ');
t=t0:h:tN;
% Primer caso---------------------------------
S10=159985;
I10=15;
I1=zeros(1,length(t));
S1=zeros(1,length(t));
I1(1)=I10;
S1(1)=S10;
% Segundo caso--------------------------------
S20=140000;
I20=20000;
I2=zeros(1,length(t));
S2=zeros(1,length(t));
I2(1)=I20;
S2(1)=S20;
% Tercer caso---------------------------------
S30=159999;
I30=1;
I3=zeros(1,length(t));
S3=zeros(1,length(t));
I3(1)=I30;
S3(1)=S30;
%datos
a=0.341;
r=0.0000218;
% Runge-Kutta de orden 4 Primer caso--------------------------------------------------------
for n=1:length(t)-1
%primer término del runge-kutta del sistema
K1s=-r*S1(n)*I1(n);
K1i=r*S1(n)*I1(n)-a*I1(n);
%segundo término del runge-kutta del sistema
K2s=-r*(S1(n)+h/2*K1s)*(I1(n)+h/2*K1i);
K2i=r*(S1(n)+h/2*K1s)*(I1(n)+h/2*K1i)-a*(I1(n)+h/2*K1i);
%tercer término del runge-kutta del sistema
K3s=-r*(S1(n)+h/2*K2s)*(I1(n)+h/2*K2i);
K3i=r*(S1(n)+h/2*K2s)*(I1(n)+h/2*K2i)-a*(I1(n)+h/2*K2i);
%cuarto término del runge-kutta del sistema
K4s=-r*(S1(n)+h/2*K3s)*(I1(n)+h/2*K3i);
K4i=r*(S1(n)+h/2*K3s)*(I1(n)+h/2*K3i)-a*(I1(n)+h/2*K3i);
S1(n+1)=S1(n)+h/6*(K1s+2*K2s+2*K3s+K4s);
I1(n+1)=I1(n)+h/6*(K1i+2*K2i+2*K3i+K4i);
end
% Runge-Kutta de orden 4 Segundo caso-------------------------------------------------------
for n=1:length(t)-1
%primer término del runge-kutta del sistema
K1s=-r*S2(n)*I2(n);
K1i=r*S2(n)*I2(n)-a*I2(n);
%segundo término del runge-kutta del sistema
K2s=-r*(S2(n)+h/2*K1s)*(I2(n)+h/2*K1i);
K2i=r*(S2(n)+h/2*K1s)*(I2(n)+h/2*K1i)-a*(I2(n)+h/2*K1i);
%tercer término del runge-kutta del sistema
K3s=-r*(S2(n)+h/2*K2s)*(I2(n)+h/2*K2i);
K3i=r*(S2(n)+h/2*K2s)*(I2(n)+h/2*K2i)-a*(I2(n)+h/2*K2i);
%cuarto término del runge-kutta del sistema
K4s=-r*(S2(n)+h/2*K3s)*(I2(n)+h/2*K3i);
K4i=r*(S2(n)+h/2*K3s)*(I2(n)+h/2*K3i)-a*(I2(n)+h/2*K3i);
S2(n+1)=S2(n)+h/6*(K1s+2*K2s+2*K3s+K4s);
I2(n+1)=I2(n)+h/6*(K1i+2*K2i+2*K3i+K4i);
end
% Runge-Kutta de orden 4 Tercer caso--------------------------------------------------------
for n=1:length(t)-1
%primer término del runge-kutta del sistema
K1s=-r*S3(n)*I3(n);
K1i=r*S3(n)*I3(n)-a*I3(n);
%segundo término del runge-kutta del sistema
K2s=-r*(S3(n)+h/2*K1s)*(I3(n)+h/2*K1i);
K2i=r*(S3(n)+h/2*K1s)*(I3(n)+h/2*K1i)-a*(I3(n)+h/2*K1i);
%tercer término del runge-kutta del sistema
K3s=-r*(S3(n)+h/2*K2s)*(I3(n)+h/2*K2i);
K3i=r*(S3(n)+h/2*K2s)*(I3(n)+h/2*K2i)-a*(I3(n)+h/2*K2i);
%cuarto término del runge-kutta del sistema
K4s=-r*(S3(n)+h/2*K3s)*(I3(n)+h/2*K3i);
K4i=r*(S3(n)+h/2*K3s)*(I3(n)+h/2*K3i)-a*(I3(n)+h/2*K3i);
S3(n+1)=S3(n)+h/6*(K1s+2*K2s+2*K3s+K4s);
I3(n+1)=I3(n)+h/6*(K1i+2*K2i+2*K3i+K4i);
end
N=S10+I10;
h1=10000;
J=0:h1:N;
D=zeros(1,length(J));
P=N*ones(1,length(J));
for i=1:length(J)
D(i)=P(i)-J(i);
end
RO=a/r;
%Para k=0
L1=RO*log(S1)-S1;
L2=RO*log(S2)-S2;
L3=RO*log(S3)-S3;
% Trayectorias
figure(1)
clf
hold on
plot(t,L1,'g','linewidth',2)
plot(t,L2,'linewidth',2)
plot(t,L3,'r','linewidth',2)
legend('Caso 1','Caso 2','Caso 3','location','best')
hold off
% S-I
figure(2)
clf
hold on
plot(S1,I1,'g','linewidth',2)
plot(S2,I2,'linewidth',2)
plot(S3,I3,'r','linewidth',2)
plot(D,J,'k')
legend('Caso 1','Caso 2','Caso 3','location','best')
hold off
.png)
.png)
Las trayectorias obtenidas en el plano de fases (I,S), están necesariamente acotadas por el triángulo de vértices (S,I)=(0,0); (S,I)=(0,N); (S,I)=(N,0). Esto ocurre porque las variables son positivas y dado que el sistema tiene tamaño fijo.
2.6 Número máximo de infectados en el primer caso del apartado 3
Para el primer caso del apartado 2.3 (S0=159 985, I0=15, r=0.0000218, a=0.341), si calculamos analíticamente el número de infectados máximos, tanto para 4 semanas de estudio como para 20, resulta ser:
I(máx) = N - rho + rho * ln (rho / S0) = S0 - a/r + a/r * ln (a / (S0*r)) = 159 985 - 15 642'2 + 15 642'2 * ln (15 642'2 / 159 985) = 107 988'0008 infectados
Y para calcularlo numéricamente, hemos empleado el programa usado en el apartado 2.3 añadiendo solamente:
format long e
Imax=max(I)
Y nos da un valor de Imax para los pasos h = 0.1, 0.01 y 0.001 siguientes:
Imax = 1.08520960052860e+005 = 108 521 infectados para h = 0.1
Imax = 1.08041178514538e+005 = 108 041 infectados para h = 0.01
Imax = 1.07993269666941e+005 = 107 993 infectados para h = 0.001
Como se puede comprobar, cuanto menor es el paso usado menor es el error del método numérico frente al analítico.
