Diferencia entre revisiones de «Discusión:Modelos Epidemiológicos Grupo 3A»

De MateWiki
Saltar a: navegación, buscar
(→‎Resolución del modelo completo mediante Runge-Kutta de orden 4)
(→‎Resolución del PVI mediante los métodos de Euler y Trapecio)
Línea 32: Línea 32:
 
[[Archivo:2Euler|100px|miniaturadeimagen|derecha|Aproximación de la función simplificada por el método de Euler]]
 
[[Archivo:2Euler|100px|miniaturadeimagen|derecha|Aproximación de la función simplificada por el método de Euler]]
 
  {{matlab|codigo=
 
  {{matlab|codigo=
 +
[[Archivo:2Euler|miniaturadeimagen|derecha]]
 
%Método de Euler
 
%Método de Euler
 
%I(t): Población de individuos infectados
 
%I(t): Población de individuos infectados

Revisión del 14:37 2 mar 2015

1 Introducción, hipótesis inciciales

En el desarrollo de una epidemia, se distinguen dos tipos de individuos. Los que ya han contraído la enfermedad, que llamaremos I(t); y los que son susceptibles de contraerla, a los que llamaremos S(t). Donde t es el tiempo. 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.

Realizaremos el estudio gracias al siguiente sistema de ecuaciones en el que se muestran las variaciones de ambos individuos 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:

Archivo:2Euler
Aproximación de la función simplificada por el método de Euler
[[Archivo:2Euler|miniaturadeimagen|derecha]]
%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,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


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. 

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


4 Resolución del modelo completo mediante Runge-Kutta de orden 4

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

Rungekutta3A15.jpg

5 Resolución del modelo completo mediante Heun con a(t) una función dependiente del tiempo

%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


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
Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Discusión:Modelos_Epidemiológicos_Grupo_3A&oldid=24388»