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).

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ón

Y la gráfica obtenida de la comparación de la condición en función de cada valor posible de C es:

Error absoluto en I(1) para cada c

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;
% Vector tiempo
t=t0:h:tN;
% Vector de infectados para Runge Kutta
y=zeros(1,N+1);
y(1)=y0;
% Vector de infectados para Heun
z=zeros(1,N+1);
z(1)=y0;
% Contactos
c=0.92;
% Bucle de Euler
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')

Gráfica de la variación de I en 20 semanas según Heun

Gráfica de la variación de I en 20 semanas según RK4

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;
% Vector tiempo
t=t0:h:tN;
N=length(t);
% Vector de infectados
y=zeros(1,N);
y(1)=y0;
c=0.92;
% 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:');
% Vector tiempo
t=t0:h:tN;
N=length(t);
% Vector de infectados
y=zeros(1,N);
y(1)=y0;
c=0.92;
% 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;
N=length(t);
c=0.92;
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)

Gráfica para 9 700 infectados iniciales

Gráfica para 10 200 infectados iniciales

Gráfica para 30 000 infectados iniciales

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

Utilizamos el método de Euler para distintos valores S0, I0 y para distintos intervalos [0,4] y [0,20] El código empleado para su cálculo ha sido el siguiente:

clc
clf
clear all
t0=0;
S0=input('Introduce los susceptiles de ser infectados inicialmente: ');
I0=input('Introduce el número 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)
plot(t,I,'g')
plot(t,R,'r')
legend('Variable S','Variable I','Variable R','Location','best')
hold off