Diferencia entre revisiones de «Modelos epidemiológicos (Grupo 3A)»

Revisión del 21:22 1 may 2016

1 Introducción y planteamiento del problema

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.

[math]\Rightarrow[/math] 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');

[math]\Rightarrow[/math] Método del trapecio

Aproximación mediante el método 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 z(i) hasta que z(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. Primero vamos a definir nuestro nuevo PVI:

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

[math]\Rightarrow[/math] Método de Euler

Aproximación mediante Euler para S=100

%Euler
clear all
%Datos
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%a: Es la constante de proporcionalidad entre los infectados por cada
%interacción
%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;
S=100;
a=0.003;
b=0.3;
c=0.01;

t(1)=t0;
y(1)=y0;

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
while y(i)>500
    
    y(i+1)=y(i)+h*(a*S*y(i)-b*y(i)-c*y(i));
    t(i+1)=t(i)+h;
    i=i+1;
end

[t',y'];

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

[math]\Rightarrow[/math] Método del trapecio

Aproximación mediante el método del trapecio para S=100

%Trapecio
clear all
%Datos
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%a: Es la constante de proporcionalidad entre los infectados por cada
%interacción
%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;
S=100;
a=0.003;
b=0.3;
c=0.01;

t(1)=t0;
z(1)=z0;

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

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

%TRAPECIO
  while z(i)>500
      z(i+1)=(z(i)*(1-(h/2)*(b+c-S*a)))/(1+h/2*(b+c-S*a));
      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 para S=100(Método del trapecio)','Location','best');

Ahora con estos nuevos datos vemos que para que se reduzca a la cuarta parte el número de infectados tendrán que pasar 138.7 días, algo lógico al fijarnos en la ED, ambas son ecuaciones decrecientes pero para el caso de S=100 la pendiente será mucho menor. Si continuamos haciendo pruebas vemos para valores superiores de s, por ejemplo, S=200 el matlab no nos da ninguna respuesta. Esto ocurre porque al aumentar S la ecuación pasa a ser creciente, por lo que lógicamente si intentas buscar un valor del número de infectados menor al inicial no lo podrás encontrar.

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 y más tarde el método de Runge-Kutta de orden cuatro con distintas discretizaciones para ver la influencia de estas.

[math]\Rightarrow[/math] Método de Euler

%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');

[math]\Rightarrow[/math] Método de Runge-Kutta

%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;

%Vector solucion
y0=[S0;I0];
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')
legend('Poblacion sana','Poblacion enferma','Location','best'); 
hold off

% 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

La dificultad que hay en el uso del método del trapecio o cualquier otro método implícito está en que las incógnitas de nuestro sistema [math](S,I)[/math] quedan implícitas en ambas ecuaciones y necesitan que se despejen. Al comenzar a realizar cálculos nos damos cuenta que el despeje es imposible, ya que las dos variables incógnita aparecen multiplicándose entre sí, quedando una variable dependiendo de la otra siempre.

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

En estos programas, sustituyendo las [math]h[/math] para 0.1, 0.01, 0.001 y 0.0001 obtendremos los siguientes resultados:

-Para [math]h=0.1[/math]

Método de Euler

Método de Runge-Kutta

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.

• Método de Euler: nuestro programa nos indica que en este caso el momento en el que los enfermos son máximos se produce al tercer día (Día 2.8) y el número de enfermos en ese momento es 517.
• Método de Runge-Kutta: el momento en el que los enfermos son máximos se produce al tercer día (Día 2.6) y el número de enfermos en ese momento es 650.

-Para [math]h=0.01[/math]

Método de Euler

Método de Runge-Kutta

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.

• Método de Euler: en este caso el momento en el que los enfermos son máximos es también en el tercer día (2.7) y el número de enfermos es 506.
• Método de Runge-Kutta: el momento en el que los enfermos son máximos se produce al tercer día (Día 2.69) y el número de enfermos en ese momento es 517.

-Para [math]h=0.001[/math]

Método de Euler

Método de Runge-Kutta

• Método de Euler: en este caso el momento en el que los enfermos son máximos se produce el tercer día (2.696) con un valor de 505 enfermos.
• Método de Runge-Kutta: el momento en el que los enfermos son máximos se produce el tercer día (2.694) con un valor de 506 enfermos.

-Para [math]h=0.0001[/math]

Método de Euler

Método de Runge-Kutta

• Método de Euler: al igual que antes, se repite el momento en el que los enfermos son máximos en el tercer día (2.695) con un valor de 505 enfermos.
• Método de Runge-Kutta: se repite el momento en el que los enfermos son máximos en el tercer día (2.6948) con un valor de 505 enfermos.

3.2 Situación inicial de [math](S _{0} , I_{0} )=(10000,40)[/math]

En estos programas, sustituyendo las [math]h[/math] para 0.1, 0.01, 0.001 y 0.0001 obtendremos los siguientes resultados:

-Para [math]h=0.1[/math]

[[Archivo:|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

• Método de Euler:
• Método de Runge-Kutta:

-Para [math]h=0.01[/math] [[Archivo:|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

• Método de Euler:
• Método de Runge-Kutta:

-Para [math]h=0.001[/math] [[Archivo:|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

• Método de Euler:
• Método de Runge-Kutta:

-Para [math]h=0.0001[/math] [[Archivo:|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

• Método de Euler:
• Método de Runge-Kutta:

4 Resolución del problema completo suponiendo [math] a [/math] una función variable con el tiempo

Para el siguiente caso supondremos que el coeficiente [math] a [/math] (tasa de contagio en población sana) es una función dada por :

[math] a(t)= \frac{0.003}{(1+t)} [/math] Suponiendo además las condiciones iniciales [math] (S_{0},I_{0})=(1640,40) [/math] , vamos a dibujar las gráficas de infectados y sanos por el método de Heun, tomando un tamaño de paso [math] h=0.01 [/math]

Población infectada-Población sana para a=a(t)

%Heun. Coeficiente variable a(t)

%Datos
S0=input('Introducir valor inicial SO: ');
I0=input('Introducir valor inicial IO: ');
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;

%Vector solucion
y0=[S0;I0];
y=zeros(2,N+1);
y(:,1)=y0;
for i=1:N
     K1S=-(0.003/(1+t(i)))*y(1,i)*y(2,i);
     K2S=-(0.003/(1+t(i)+h))*(y(1,i)+K1S*h)*(y(2,i)+K1S*h);
     y(1,i+1)=y(1,i)+(h/2)*(K1S+K2S);
     K1I=(0.003/(1+t(i)))*y(1,i)*y(2,i)-b*y(2,i)-c*y(2,i);
     K2I=(0.003/(1+t(i)+h))*(y(1,i)+K1I*h)*(y(2,i)+K1I*h)-b*(y(2,i)+K1I*h)-c*(y(2,i)+K1I*h);
     y(2,i+1)=y(2,i)+(h/2)*(K1I+K2I);
end

%Tabla de resutados
[t',y(1,:)',y(2,:)']
%Graficos
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'k')
title('Método de Heun');
legend('Poblacion sana','Población infectada','Location','best'); 
hold off

% 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

Como se puede ver la población sana se infecta rápidamente y a partir del día 15 casi toda la población ha fallecido. Al segundo día [math](t=2,15)[/math] observamos que el número de infectados es de 994 individuos [math](I=994,3)[/math], siendo este el máximo que se alcanzará.

A medida que va pasando el tiempo el coeficiente [math] a(t) [/math] (tasa de contagio en población sana) va a ir decreciendo poco a poco hasta aproximarse a cero según podemos ver en el siguiente gráfico: [math] \lim_{t\to \infty } a(t)=0[/math]

a(t) para t=[0,40]

%tasa de contagio en poblacion sana a(t)

%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector a
a=zeros(1,N+1);
for i=1:N+1
    a(i)=0.003/(1+t(i));
end

%Grafica
plot(t,a)

Poblaciones según la tasa de contagio a(t)

Ambas poblaciones parten del mismo inicio (para [math]t=0[/math] las poblaciones son iguales). Podemos observar que la población sana (susceptible de contraer la enfermedad) en el caso en el que la tasa de contagio en población sana [math](a)[/math] es constante [math](a=0.003)[/math] se ve reducida antes que en el caso en el que depende del tiempo [math](a=a(t))[/math].

[math]a(t)[/math] disminuye a la vez que transcurre el tiempo, por lo que en general las pendientes tanto de la población sana como de la infectada serán menores para este caso que para el caso en el que [math]a=0.003[/math]

5 Calibración del parámetro [math]a[/math] en base a una experiencia previa