Diferencia entre revisiones de «Modelos Epidemiológicos Grupo 3A»
(→Resolución del modelo completo mediante Runge-Kutta de orden 4) |
(→Resolución del modelo completo mediante Runge-Kutta de orden 4) |
||
| Línea 275: | Línea 275: | ||
}} | }} | ||
A continuación, comparamos las gráficas obtenidas con (S<sub>0</sub>, I<sub>0</sub>)=(800, 20) entre Runge-Kutta y Euler: | A continuación, comparamos las gráficas obtenidas con (S<sub>0</sub>, I<sub>0</sub>)=(800, 20) entre Runge-Kutta y Euler: | ||
| − | |||
[[File:Rungekutta3A15.jpg|500px|thumb|left|Aproximación del sistema mediante RK-4]] | [[File:Rungekutta3A15.jpg|500px|thumb|left|Aproximación del sistema mediante RK-4]] | ||
[[File:Euler3A15.jpg|500px|thumb|right|Aproximación del sistema mediante Euler]] | [[File:Euler3A15.jpg|500px|thumb|right|Aproximación del sistema mediante Euler]] | ||
| − | |||
Al trabajar con los valores (S<sub>0</sub>, I<sub>0</sub>)=(10000, 40) las aproximaciones mediante Euler y RK4 muestran gráficas incongruentes con la evolución de personas susceptibles e infectadas observada hasta ahora, en las aproximaciones anteriores hemos obtenido gráficas de tipo exponencial con un máximo de infectados y un descenso gradual de ambas poblaciones cuando el tiempo toma valores altos, sin embargo, en este caso mediante Euler hemos obtenido gráficas con una sucesión de máximos y mínimos en intervalos de tiempo muy pequeños( en cuestión de 1 día las cifras de infectados oscilan entre pocos cientos y más de 10000 alternativamente) y además las personas susceptibles adquieren valores negativos(Físicamente imposible) lo cual nos indica que '''el método de Euler es inadecuado''' para la aproximación requerida, seguramente debido a la inestabilidad que presenta en la resolución de sistemas.Con el método de Runge-Kutta se obtiene de la misma manera un resultado gráfico desconcertante con respecto a las gráficas expuestas anteriormente ofreciendo resultados que nos son equiparables a una evolución epidemiológica real. Otro posible motivo del error en las gráficas es la incompatibilidad de la pareja de valores iniciales con la aproximación de una epidemia expuesta en este trabajo, es decir, para valores muy altos de S(t) la modelización de una epidemia no puede ser aproximada por un sistema de ecuaciones diferenciales como el proporcionado para este trabajo(Obtenemos valores que no concuerdan con la realidad). | Al trabajar con los valores (S<sub>0</sub>, I<sub>0</sub>)=(10000, 40) las aproximaciones mediante Euler y RK4 muestran gráficas incongruentes con la evolución de personas susceptibles e infectadas observada hasta ahora, en las aproximaciones anteriores hemos obtenido gráficas de tipo exponencial con un máximo de infectados y un descenso gradual de ambas poblaciones cuando el tiempo toma valores altos, sin embargo, en este caso mediante Euler hemos obtenido gráficas con una sucesión de máximos y mínimos en intervalos de tiempo muy pequeños( en cuestión de 1 día las cifras de infectados oscilan entre pocos cientos y más de 10000 alternativamente) y además las personas susceptibles adquieren valores negativos(Físicamente imposible) lo cual nos indica que '''el método de Euler es inadecuado''' para la aproximación requerida, seguramente debido a la inestabilidad que presenta en la resolución de sistemas.Con el método de Runge-Kutta se obtiene de la misma manera un resultado gráfico desconcertante con respecto a las gráficas expuestas anteriormente ofreciendo resultados que nos son equiparables a una evolución epidemiológica real. Otro posible motivo del error en las gráficas es la incompatibilidad de la pareja de valores iniciales con la aproximación de una epidemia expuesta en este trabajo, es decir, para valores muy altos de S(t) la modelización de una epidemia no puede ser aproximada por un sistema de ecuaciones diferenciales como el proporcionado para este trabajo(Obtenemos valores que no concuerdan con la realidad). | ||
A modo de aclaración se adjuntan dichas gráficas a continuación: | A modo de aclaración se adjuntan dichas gráficas a continuación: | ||
Revisión del 22:52 4 mar 2015
Contenido
- 1 Introducción e hipótesis iniciales
- 2 Resolución del PVI mediante los métodos de Euler y Trapecio
- 3 Resolución del modelo completo mediante Euler
- 4 Resolución del modelo completo mediante Runge-Kutta de orden 4
- 5 Resolución del modelo completo mediante Heun con a(t) una función dependiente del tiempo
- 6 Calibración del parámetro 'a' en base a una experiencia previa
1 Introducción e hipótesis iniciales
En el desarrollo de una epidemia, se distinguen dos tipos de individuos. Los que ya han contraído la enfermedad(infectados), que llamaremos I(t); y los que son susceptibles de contraerla, a los que llamaremos S(t). Donde t es la variable temporal. Se dan dos hipótesis para realizar este estudio:
1. La población de personas infectadas se altera por el fallecimiento o la cura de las mismas. En ambos casos, la tasa de cambio depende del número de personas infectadas.
2. La tasa de individuos que pasan de ser susceptibles a contraer la enfermedad a estar infectados es proporcional a la interacción entre el número de individuos en ambas clases.
Mediante el siguiente sistema de ecuaciones diferenciales se muestran las variaciones de ambas poblaciones respecto al tiempo t:
\begin{matrix}
\frac{dS}{dt} = -aSI \\
\frac{dI}{dt} = aSI-bI-cI
\end{matrix}
donde a,b y c son parámetros. Los interpretamos como:
\frac{dI}{dt} es la variacion de la población infectada. Esta población únicamente se ve alterada por:
- aSI ( con signo positivo), donde SI es la interacción entre ambas poblaciones.
- -bI (con signo negativo) que son los que fallecen. Por lo tanto, 'b' es la tasa de fallecimiento.
- -cI (con signo negativo) que son los que se han curado. Es decir, 'c' es la tasa de curación.
Para interpretar 'a' observamos el sistema y la segunda hipótesis:
\frac{dI}{dt} = (aS-b-c)I
aS serán los individuos susceptibles convertidos en infectados. Por tanto, a es la tasa de contagio en la población susceptible.
2 Resolución del PVI mediante los métodos de Euler y Trapecio
Antes de pasar a la resolución completa del sistema, lo simplificaremos asignando el valor "cero" a la variable que representa a los susceptibles de contraer la enfermedad (S=0). Con ello, reducimos el sistema a una única ecuación:
¡Germán!Mete ecuación!
Mediante sentencias Matlab y la asignación a los coeficientes a,b y c de los valores constantes requeridos, resolvemos la ecuación diferencial simplificada aplicando los métodos de Euler y Trapecio (para ello se ha utilizado, en ambos métodos, la espaciación requerida de h=0.1):
%Método de Euler
%I(t): Población de individuos infectados
%S(t): Población susceptible de ser infectada
%los parámetros a=Tasa de contagio entre S e I,b=Tasa de
%fallecimiento de personas infectadas,c=Tasa de cura de personas infectadas
clear all clc
h=0.1;
format long
%Inicializo los vectores t,i e y asignando a la primera posición los valores
%y0 ,t0
t0=0;
y0=2000;
t(1)=t0;
y(1)=y0;
i=1;
%Defino las ctes b y c
b=0.3;
c=0.01;
%Bucle while: calculara valores de y(i) mientras y(i)>500
while y(i)>500
y(i+1)=y(i)+h*(-y(i)*(b+c));%Euler
t(i+1)=t(i)+h;
i=i+1;
end
[t',y']
disp('Tiempo final:')
disp(t(end));
plot(t,y);
legend('Población infectada(Euler)','Location','best');%Optimizar
Resolvemos de nuevo la ecuación simplificada, en este caso, mediante el método del trapecio:
%Método del Trapecio
%I(t): Población de individuos infectados
%S(t: Población susceptible de ser infectada
%los parámetros a=Tasa de contagio,b=Tasa de
%fallecimiento de personas infectadas,c=Tasa de cura de personas infectadas
clear all clc
format long
h=0.1;
%Inicializo los vectores t,i,z asignando a la primera posición los valores
%z0 ,t0
t0=0;
z0=2000;
t(1)=t0;
z(1)=z0;
i=1;
%Defino las ctes.b,c
b=0.3;
c=0.01;
%Bucle while: calculara valores de z(i) mientras z(i)>500
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 final:')
disp(t(end));
plot(t,z,'r');
legend('Población infectada(Trapecio)','Location','best');%Optimizar
Como adelantábamos al principio de este apartado, se va a proceder a la resolución del sistema en su conjunto, con el matiz de que ninguna ecuación se anulará, ya que tomamos como valor "S" de los susceptibles de ser infectados un valor costante de "cien" (S=100). La resolución, al igual que en caso anterior, se ejecutará mediante los métodos de Euler y del Trapecio:
%Método de Euler
%I(t): Población de individuos infectados
%S(t): Población susceptible de ser infectada
%los parámetros a=Tasa de contagio,b=Tasa de
%fallecimiento de personas infectadas,c=Tasa de cura de personas infectadas
clear all clc
format long
h=0.1;
%Inicializo los vectores t,i e y asignando a la primera posición los valores
%y0 ,t0
i=1;
t0=0;
y0=2000;
t(1)=t0;
y(1)=y0;
%Defino las ctes. S,a,b,c
S=input('Introduce población susceptible de contagio:');
a=0.003;
b=0.3;
c=0.01;
%Bucle while: calculara valores de y(i) mientras y(i)>500
while y(i)>500
y(i+1)=y(i)+h*(-y(i)*(b+c-S*a));%Euler
t(i+1)=t(i)+h;
i=i+1;
end
[t',y']
disp('Tiempo final:')
disp(t(end));
plot(t,y);
De forma análoga aproximamos ese mismo sistema, sin variar el valor de las constantes, de la espaciación, ni de "S" (S=100), pero en este caso utilizando el método del trapecio:
%Método del Trapecio
%I(t): Población de individuos infectados
%S(t): Población susceptible de ser infectada
%los parámetros a=Tasa de contagio,b=Tasa de
%fallecimiento de personas infectadas,c=Tasa de cura de personas infectadas
clear all clc
format long
h=0.1;
%Inicializo los vectores t,i y z asignando a la primera posición los valores
%z0 ,t0
t0=0;
z0=2000;
t(1)=t0;
z(1)=z0;
i=1;
%Defino las ctes. S,a,b,c
S=input('Introduce población susceptible de contagio:');
a=0.003;
b=0.3;
c=0.01;
%Bucle while: calculara valores de z(i) mientras z(i)>500
while z(i)>500
z(i+1)=(z(i)*(1-(h/2)*(b+c-S*a)))/(1+h/2*(b+c-S*a));%trapecio
t(i+1)=t(i)+h;
i=i+1;
end
[t',z']
disp('Tiempo final:')
disp(t(end));
plot(t,z,'r');
legend('Población infectada(Trapecio)','Location','best');%Optimizar
3 Resolución del modelo completo mediante Euler
En este apartado resolveremos el sistema completo y analizaremos las gráficas para diferentes valores de las variables. Los valores iniciales que tomamos son (S0, I0)=(800, 20) y (S0, I0)=(10000, 40); el tiempo estará en el intervalo [0,40] y como pasos de discretización temporal (h) 10-1, 10-2, 10-3, 10-4. Por lo tanto, el código de matlab utilizado en este apartado es:
%Trabajo Ecuaciones diferenciales
% Considerando:
%I(t): Población de individuos infectados
%S(t): Población susceptible de ser infectada
%los parámetros a=Tasa de contagio,b=Tasa de
%fallecimiento de personas infectadas,c=Tasa de cura de personas infectadas
clear all clc
format long
h=input('Introduce tamaño de paso:');
S0=input('introduce un valor inicial S0:');
I0=input('introduce un valor inicial I0:');
a=0.003;
b=0.3;
c=0.01;
%Intervalo de tiempos
t0=0;
tN=40;
%Calculamos número de subintervalos
N=(tN-t0)/h;
%Creo un vector t
t=t0:h:tN;
%preparo vectores S,I para guardar las soluciones
S=zeros(1,N+1);%vector incognita 1
I=zeros(1,N+1);%vector incognita 2
S(1)=S0;
I(1)=I0;
for i=1:N
S(i+1)=S(i)+h*(-a*S(i)*I(i));%Euler para S
I(i+1)=I(i)+h*(a*S(i)*I(i)-b*I(i)-c*I(i));%Euler para I
end
[t',S',I']
hold on
plot (t,S)
plot(t,I,'r')
legend('Susceptibles','Infectados','Location','best');%Optimizar
hold off
Como hemos dicho aquí están las ocho gráficas referentes a los distintos valores. Como se puede observar, en las gráficas que tienen como valores iniciales (S0, I0)=(1000, 40) no vemos como descienden y ascienden los susceptibles y los infectados. Todo los contrario pasa con las gráficas que tienen como valores iniciales (S0, I0)=(800, 20), en ellas si que podemos observar que cuando aumenta el número de infectados, disminuye el número de susceptibles, hasta que llegan más o menos a los 20 días, donde el número de susceptibles y de infectados empieza a ser prácticamente el mismo y tienden los dos a cero.
4 Resolución del modelo completo mediante Runge-Kutta de orden 4
En esta sección se pretende resolver de nuevo el mismo problema que en el apartado anterior, sin embargo, esta vez lo haremos con el método de Runge-Kutta de cuarto orden con el objetivo de compararlo los resultados obtenidos con el método de Euler.El código matlab correspondiente queda expuesto a continuación:
%Trabajo Ecuaciones diferenciales
% Considerando:
%I(t): Población de individuos infectados
%S(t: Población susceptible de ser infectada
%los parámetros a=Tasa de contagio entre S e I,b=Tasa de
%fallecimiento de personas infectadas,c=Tasa de cura de personas infectadas
clear all clc
format long
S0=input('introduce un valor inicial S0:');
I0=input('introduce un valor inicial I0:');
h=input('introduce tamano de paso:');
a=0.003;
b=0.3;
c=0.01;
%Intervalo de tiempos
t0=0;
tN=40;
%Calculamos número de subintervalos
N=(tN-t0)/h;
%Creo un vector t
t=t0:h:tN;
%preparo vectores S,I para guardar las soluciones uso el comando zeros(filas,columnas)
S=zeros(1,N+1);%vector incognita 1
I=zeros(1,N+1);%vector incognita 2
S(1)=S0;
I(1)=I0;
for i=1:N
K1=-a*S(i)*I(i);
K2=-a*(S(i)+(1/2)*K1*h)*(I(i)+(1/2)*K1*h);
K3=-a*(S(i)+(1/2)*K2*h)*(I(i)+(1/2)*K2*h);
K4=-a*(S(i)+K3*h)*(I(i)+K3*h);
S(i+1)=S(i)+(h/6)*(K1+2*K2+2*K3+K4);%Runge-Kutta orden 4 para S
K5=a*S(i)*I(i)-b*I(i)-c*I(i);
K6=a*(S(i)+(1/2)*K5*h)*(I(i)+(1/2)*K5*h)-b*(I(i)+(1/2)*K5*h)-c*(I(i)+(1/2)*K5*h);
K7=a*(S(i)+(1/2)*K6*h)*(I(i)+(1/2)*K6*h)-b*(I(i)+(1/2)*K6*h)-c*(I(i)+(1/2)*K6*h);
K8=a*(S(i)+K7*h)*(I(i)+K7*h)-b*(I(i)+K7*h)-c*(I(i)+K7*h);
I(i+1)=I(i)+(h/6)*(K5+2*K6+2*K7+K8);%Runge-Kutta orden 4 para I
end
[t',S',I']
hold on
plot (t,S)
plot(t,I,'r')
legend('Susceptibles','Infectados','Location','best');%Optimizar
hold offA continuación, comparamos las gráficas obtenidas con (S0, I0)=(800, 20) entre Runge-Kutta y Euler:
Al trabajar con los valores (S0, I0)=(10000, 40) las aproximaciones mediante Euler y RK4 muestran gráficas incongruentes con la evolución de personas susceptibles e infectadas observada hasta ahora, en las aproximaciones anteriores hemos obtenido gráficas de tipo exponencial con un máximo de infectados y un descenso gradual de ambas poblaciones cuando el tiempo toma valores altos, sin embargo, en este caso mediante Euler hemos obtenido gráficas con una sucesión de máximos y mínimos en intervalos de tiempo muy pequeños( en cuestión de 1 día las cifras de infectados oscilan entre pocos cientos y más de 10000 alternativamente) y además las personas susceptibles adquieren valores negativos(Físicamente imposible) lo cual nos indica que el método de Euler es inadecuado para la aproximación requerida, seguramente debido a la inestabilidad que presenta en la resolución de sistemas.Con el método de Runge-Kutta se obtiene de la misma manera un resultado gráfico desconcertante con respecto a las gráficas expuestas anteriormente ofreciendo resultados que nos son equiparables a una evolución epidemiológica real. Otro posible motivo del error en las gráficas es la incompatibilidad de la pareja de valores iniciales con la aproximación de una epidemia expuesta en este trabajo, es decir, para valores muy altos de S(t) la modelización de una epidemia no puede ser aproximada por un sistema de ecuaciones diferenciales como el proporcionado para este trabajo(Obtenemos valores que no concuerdan con la realidad). A modo de aclaración se adjuntan dichas gráficas a continuación:
5 Resolución del modelo completo mediante Heun con a(t) una función dependiente del tiempo
Hasta ahora, el coeficiente a ( tasa de contagio en población susceptible) ha sido dada como un coeficiente constante. Ahora 'a' es una función dependiente del tiempo t (coeficiente variable) dada por:
a(t)=\frac{0.003}{1+t}
%Trabajo Ecuaciones diferenciales
% Considerando:
%I(t): Población de individuos infectados
%S(t: Población susceptible de ser infectada
%los parámetros a=Tasa de contagio,b=Tasa de
%fallecimiento de personas infectadas,c=Tasa de cura de personas infectadas
clear all clc
format long
S0=1600;
I0=40;
h=input('introduce tamano de paso:');
b=0.3;
c=0.01;
%Intervalo de tiempos
t0=0;
tN=40;
%Calculamos número de subintervalos
N=(tN-t0)/h;
%Creo un vector t
t=t0:h:tN;
%preparo vector y para guardar las soluciones uso el comando zeros(filas,columnas)
S=zeros(1,N+1);%vector incognita 1
I=zeros(1,N+1);%vector incognita 2
S(1)=S0;
I(1)=I0;
for i=1:N
K1=-(0.003./(1+t(i)))*(S(i)*I(i));
K2=-(0.003./(1+t(i))).*(S(i)+K1*h)*(I(i)+K1*h);
K3=(0.003./(1+t(i))).*(S(i)*I(i)-b*I(i)-c*I(i));
K4=(0.003./(1+t(i))).*(S(i)+K3*h)*(I(i)+K3*h)-b*(I(i)+K3*h)-c*(I(i)+K3*h);
S(i+1)=S(i)+(h/2)*(K1+K2);
I(i+1)=I(i)+(h/2)*(K3+K4);
end
[t',S',I']
hold on
plot(t,S);
plot(t,I,'r');
legend('Susceptibles','Infectados','Location','best');%Optimizar
hold off
El gráfico de a(t) es el siguiente:
a(t) se trata de una función no definida en t=-1. Esto no nos importa, ya que t es el tiempo, que siempre es positivo y concretamente toma valores entre 0 y 40. Como la t>0 siempre:
\lim_{t\to \infty } a(t))= 0
Esto quiere decir que, debido al coeficiente variable a(t), la variación de las poblaciones no se verá alterada con respecto a casos anteriores.
6 Calibración del parámetro 'a' en base a una experiencia previa
clear clc
format long
%Defino valores iniciales
S0=1600;
I0=40;
h1=0.0001;%tamaño de paso para intervalo de a
a=0.0005:h1:0.002;
b=0.3;
c=0.01;
h=0.1;%tamaño de paso para t
%Intervalo de tiempos
t0=0;
tN=40;
%Calculamos número de subintervalos
N=(tN-t0)/h;
%Creo un vector t
t=t0:h:tN;
%preparo vector S e I para guardar las soluciones
S=zeros(1,N+1);%vector incognita 1
I=zeros(1,N+1);%vector incognita 2
S(1)=S0;
I(1)=I0;
%preparo vectores m1(máximos de I) m2(tiempos de los máximos)
m1=zeros(1,length(a));
m2=zeros(1,length(a));
%Bucle for: calcula valores maximos de I para cada valor 'a' del intervalo y los almacena
%junto con su tiempo correspondiente en los vectores m1 y m2
%respectivamente.Metodo numérico de Heun
for j=1:length(a)
for i=1:N
K1=-(a(j))*(S(i)*I(i));
K2=-(a(j)).*(S(i)+K1*h)*(I(i)+K1*h);
K3=(a(j)).*(S(i)*I(i)-b*I(i)-c*I(i));
K4=(a(j)).*(S(i)+K3*h)*(I(i)+K3*h)-b*(I(i)+K3*h)-c*(I(i)+K3*h);
S(i+1)=S(i)+(h/2)*(K1+K2);
I(i+1)=I(i)+(h/2)*(K3+K4);
m1(j)=max(I);
m2(j)=t(max(I)==I);
end
%Bucle if: si la diferencia entre el valor de tiempo para el máximo y 5 es
%menor o igual que 0.1 muestra el valor de a, maximo de I y tiempo
%correspondiente para los que esto sucede
if abs(m2(j)-5)<=0.1
disp('Valor de a para el cual el máximo está mas próximo a 5:');
disp(a(j));
disp('Valor máximo de I correspondiente:');
disp((m1(j)));
disp('Valor de tiempo correspondiente:');
disp((m2(j)));
end
end
%Grafica con los máximos de I y los tiempos correspondientes
plot(m2,m1,'*');
legend('Máximos de infectados','Location','best');%Optimizar
