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 en el triángulo (0,0)-(0,N)-(N,0)
- 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)
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
2.3 Resolucióin numérica del sistema SIR
Para resolver numéricamente el sistema utilizamos el método de Euler para distintos valores iniciales de [S0,I0] ([159 985,15],[140 000,20 000],[159 999,1]) y para 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:
clc, clf, clear all
t0=0;
S0=input('Introduce los susceptibles de ser infectados inicialmente: ');
I0=input('Introduce el numero de infectados inicialmente: ');
R0=0;
tN=input('Introduce el numero de semanas de estudio: ');
h=input('Introduce el paso: ');
N=round((tN-t0)/h); %numero de subintervalos
% Definimos la variable independiente---------------------------------------------------------------------------
t=t0:h:tN;
% Definimos las variables dependientes--------------------------------------------------------------------------
S=zeros(1,N+1); % Susceptibles de ser infectados (primera incógnita)
I=zeros(1,N+1); % Infectados (segunda incógnita)
R=zeros(1,N+1); % Removidos (tercera incógnita)
% Datos iniciales-----------------------------------------------------------------------------------------------
S(1)=S0;
I(1)=I0;
R(1)=R0;
% Bucle Euler---------------------------------------------------------------------------------------------------
for i=1:N
S(i+1)=S(i)+h*(S(i)*(-0.0000218)*I(i)); %Euler primera ecuación
I(i+1)=I(i)+h*(0.0000218*S(i)*I(i)-0.341*I(i)); %Euler segunda ecuación
R(i+1)=R(i)+h*(0.341*I(i)); %Euler tercera ecuación
end
% Gráficas de S, I y R------------------------------------------------------------------------------------------
hold on
plot(t,S,'Linewidth',3)
plot(t,I,'g',Linewidth',3)
plot(t,R,'r',Linewidth',3)
legend('Variable S','Variable I','Variable R','Location','best')
hold off
Para resolver numéricamente el sistema utilizamos el método de Runge-Kutta para distintos valores iniciales de [S0,I0] ([159 985,15],[140 000,20 000],[159 999,1]) y para 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
% 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: ');
I=zeros(1,length(t));
R=zeros(1,length(t));
S=zeros(1,length(t));
I(1)=I0;
S(1)=S0;
R(1)=R0;
%datos
a=0.341;
r=0.0000218;
% Runge-Kutta de orden 4
for n=1:length(t)-1
%primer término del runge-kutta del sistema
K1s=-r*S(n)*I(n);
K1i=r*S(n)*I(n)-a*I(n);
K1r=a*I(n);
%segundo término del runge-kutta del sistema
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 del runge-kutta del sistema
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 del runge-kutta del sistema
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
hold on
plot(t,S,'g','Linewidth',3)
hold on
plot(t,I,'Linewidth',3)
hold on
plot(t,R,'r','Linewidth',3)
2.4 Interpretación de la epidemia
2.5 Trayectorias en el triángulo (0,0)-(0,N)-(N,0)
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.
