Modelos epidemiológicos (Grupo 3A)

Trabajo realizado por estudiantes
Título	Modelos Epidemológicos (Grupo 3A)
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2015-16
Autores	Ignacio Mollá Carcaño Pablo Revuelta Aragón David González Hernández Jose María García Rodríguez Alejandro Martínez Gamonal
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

2 Resolución del problema con una sola ecuación diferencial

Primero vamamos a estudiar el modelo epidemiológico de una población en la que incialmente las 2000 personas que conforman la población se encuentran infectadas. Para hacer el seguimiento de la enfermedad contamos tanto con la constante de proporcionalidad de personas que se curan (b=0.3) y la constante de proporcionalidad del número de personas que fallecen (c=0.01). De esta forma nos quedará el siguiente PVI:

[math] \begin{cases}   \frac{dI}{dt} = -0.31 I \\I_{0}  = 2000\end{cases}  [/math]

A continuación vamos a proceder a resolverlo por el método de Euler y trapecio asignandole un tamaño de paso de h=0.1 y veremos cuanto tiempo tardará en reducirse el número de infectados a la cuarta parte.

-Método de Euler

Aproximación mediante Euler

clear all
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
y0=2000;
b=0.3;
c=0.01;

t(1)=t0;
ye(1)=y0; %Vector solución Euler

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de y(i) hasta que y(i)>500

%EULER

while ye(i)>500
    
    ye(i+1)=ye(i)+h*(-b*ye(i)-c*ye(i)); 
    t(i+1)=t(i)+h;
    i=i+1;
end

[t',ye'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,ye)
legend('Población infectada(Método de Euler)','Location','best');

-Método del trapecio

Aproximación mediante del trapecio

clear all
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
z0=2000;
b=0.3;
c=0.01;

t(1)=t0;
z(1)=z0; %Vector solución trapecio

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de y(i) hasta que y(i)>500

%TRAPECIO
while z(i)>500
     z(i+1)=(z(i)*(1-(h/2)*(b+c)))/(1+h/2*(b+c));%trapecio
     t(i+1)=t(i)+h;
     i=i+1;
end

[t',z'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,z)
legend('Población infectada(Método del trapecio)','Location','best');

Con todo esto podemos determinar que el número de infectados se reducirá a la cuarta parte en 4.5 días, esto ocurrirá o bien por el fallecimiento o por la curación de dicha enfermedad. Cabe señalar que el número de curados será mayor que el de defunciones, esto se puede ver a simple vista viendo que b>a. De cara a nuestro estudio de esta enfermedad vamos a ver como afectaría al desarrollo de la misma si introducimos 100 sujetos sanos en la población, tomando una constante de interacción entre personas infectadas y sanas de a=0.003. Procederemos igual que antes a resolverlo por el método de Euler y trapecio con un tamaño de paso de 0.1.

3 Resolución del problema completo

Vamos a resolver el problema completo para ver como evolucionarían las poblaciones en un periodo de 40 días. Primero supongamos una población infectada inicial de 20 individuos y una población en riesgo de contagio de 800. Después supondremos el caso de 40 individuos enfermos y 10000 en riesgo de contagio. Para resolverlo usaremos el método de Euler con distintas discretizaciones y así ver la influencia de estas.

3.1 Situación inicial de [math](S _{0} , I_{0} )=(800,20)[/math]

%Euler
clear all
%Datos
%Introducimos los valores de las constantes
a=0.003;
b=0.3;
c=0.01;
h=input('Introducir valores de h:'); 
t0=0;
y0=input('Introducir vector [So,Io]:');
tN=40;
% calculamos los subintervalos
t=t0:h:tN;          
N=(tN-t0)/h;
y=zeros(2,length(t));
y(:,1)=y0';
for i= 1:N;
    y(:,i+1)=y(:,i)+h*[-a*y(1,i).*y(2,i);a*y(1,i).*y(2,i)-b*y(2,i)-c*y(2,i)];
end

% grafico
plot(t,y)
% vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h
legend('Población sana','Location','best','Población enferma','Location','best');

Y en este programa, sustituyendo las h para 0.1 0.01 0.001 y 0.0001 obtendremos los siguientes resultados:

-Para h=0.1

Campo del desplazamiento u

Como se puede ver la población sana se infecta rápidamente y a partir del día 15 prácticamente toda la población ha fallecido. Nuestro programa nos indica que en este caso el momento de máximos enfermos se produce al final del tercer día (Día 2.8) y el número de enfermos en ese momento es 517.

-Para h=0.01

Campo del desplazamiento u

En el gráfico apenas se encuentra diferencia, lo cual es lógico porque es el mismo problema, pero si sobrepusiéramos las gráficas, lograríamos ver como se va ajustando a la realidad en función del aumento de h. En este caso el momento máximo de enfermos es también en el tercer día (2.7) y el número de enfermos es 506.

-Para h=0.001

Campo del desplazamiento u

Podemos decir lo mismo que antes. En este caso el máximo de enfermos se produce el tercer día (2.696) con un valor de 505 enfermos.

-Para h=0.0001

Campo del desplazamiento u

Al igual que antes, se repite el momento de máximos en el tercer día (2.695) con un valor de 505 enfermos.

Ahora vamos a resolver de nuevo el problema completo para ver como evolucionarían las poblaciones en un periodo de 40 días, pero esta vez utilizando el método de Runge-Kutta de cuarto orden.

%Runge-Kutta
clear all
%Datos
S0=input('Introducir valor inicial SO: ');
I0=input('Introducir valor inicial IO: ');
a=0.003;
b=0.3;
c=0.01;
%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;
y0=[S0;I0]; %Condiciones iniciales
y=zeros(2,N+1);
y(:,1)=y0;
for i=1:N
     K1S=-a*y(1,i)*y(2,i);
     K2S=-a*(y(1,i)+(1/2)*K1S*h)*(y(2,i)+(1/2)*K1S*h);
     K3S=-a*(y(1,i)+(1/2)*K2S*h)*(y(2,i)+(1/2)*K2S*h);
     K4S=-a*(y(1,i)+K3S*h)*(y(2,i)+K3S*h);
     y(1,i+1)=y(1,i)+(h/6)*(K1S+2*K2S+2*K3S+K4S);
     K1I=a*y(1,i)*y(2,i)-b*y(2,i)-c*y(2,i);
     K2I=a*(y(1,i)+(1/2)*K1I*h)*(y(2,i)+(1/2)*K1I*h)-b*(y(2,i)+(1/2)*K1I*h)-c*(y(2,i)+(1/2)*K1I*h);
     K3I=a*(y(1,i)+(1/2)*K2I*h)*(y(2,i)+(1/2)*K2I*h)-b*(y(2,i)+(1/2)*K2I*h)-c*(y(2,i)+(1/2)*K2I*h);
     K4I=a*(y(1,i)+K3I*h)*(y(2,i)+K3I*h)-b*(y(2,i)+K3I*h)-c*(y(2,i)+K3I*h);
     y(2,i+1)=y(2,i)+(h/6)*(K1I+2*K2I+2*K3I+K4I);
end
%Tabla de resutados
[t',y(1,:)',y(2,:)']
%Gráfico
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'k')
title('Método de Runge-Kutta');
legend('Suceptibles de contraer la enfermedad','Infectados de la enfermedad','Location','best'); 
hold off

-h=0.1

Campo del desplazamiento u