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

Revisión del 13:01 30 abr 2016

1 Introducción y planteamiento del problema

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

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 /I_{0}=20 y S_{0}=800/)

%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