Ver el código de «Modelos epidemiológicos Grupo 6C»

{{ TrabajoED | Modelos epidemiológicos. Grupo 6-C | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | 
Manuel Morales López 1175
David Toledo Menéndez 1228
Sergio Rodríguez Torcal 994
Jose María Rodríguez Vicente 1213
Lourdes Sánchez-Ocaña Merino 1248
Jorge Villa Lobo 1237 }}
==Enunciado==

En el desarrollo de una epidemia se distinguen dos tipos
de individuos: los que ya han contraido la enfermedad o infectados I, y los que son susceptibles de
contraerla por encontrarse en zona de riesgo S. Supongamos que se dan las siguientes hip´otesis:
1. La poblaci´on de personas infectadas se altera por el fallecimiento o la cura de las mismas. En
ambos casos, la tasa de cambio depende del n´umero 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´on entre el n´umero de individuos en ambas clases.
Consideramos las variables: t tiempo, S(t) poblaci´on de individuos susceptibles a contraer la enfermedad,
I(t) poblaci´on de individuos infectados; y el sistema:
dS
dt = −aSI
dI
dt = aSI − bI − cI
donde a, b, c son parametros.
 
==Interpretación de parámetros==
En el problema: el coeficiente "a" es la tasa de infectados por contagio, "b" la de muertos y "c" la de curados.


==Método de Euler y Trapecio con S=0==
{{matlab|codigo=
%Apartado 2 trabajo1 2015
clear all
t0=0;
tN=20;
y0=2000;
h=0.1;
t=t0:h:tN;
N=(tN-t0)/h;
y=zeros(1,N+1); %euler 
y(1)=y0;
z=zeros(1,N+1); %trapecio 
z(1)=y0;
w=zeros(1,N+1);
w(1)=y0;
x=zeros(1,N+1);
x(1)=y0;
for i=1:N
y(i+1)=y(i)-h*(0.3+0.001)*y(i);
z(i+1)=(1/(1-h/2*(-0.31)))*(z(i))+h/2*(-0.31*z(i));
end
[t',y',z'];
i=1;
while abs(y0-w(i))<3/4*y0
  w(i+1)=(1/(1-h/2*(-0.31)))*(w(i))+h/2*(-0.31*w(i)); %trapecio, evitamos la variable auxiliar yy
  t(i+1)=t(i)+h;
  i=i+1;
end
disp('Tiempo final con trapecio:')
disp(t(i))
%gráfico
i=1;
while abs(y0-x(i))<3/4*y0
  x(i+1)=x(i)-h*(0.3+0.001)*x(i); %euler, evitamos la variable auxiliar yy
  t(i+1)=t(i)+h;
  i=i+1;
end
disp('Tiempo final con Euler:')
disp(t(i))
hold on
plot(t,y,'b')
plot(t,z,'r')
plot(t,w,'g')
plot(t,x,'k')
legend('Euler','Trapecio','Tiempo que tarda por Trapecio','Tiempo que tarda por Euler','Location','best'); 
hold off}}
Esta gráfica nos muestra que en t=4.60 se alcanza la condición final, que el número de infectados se reduzca a 500 mientras que con el trapecio el tiempo final es t=4.50
[[Archivo:euler1111.jpg|520px|thumb|left|Tiempo en llegar a 500 infectados]]




































==Método de Euler y Trapecio con S=100==
Se puede interpretar como que a partir de un valor limite de "S" entre 100 y 200, el numero de infectados se mantiene constante en el tiempo. De igual forma, si "S" es  menor que este valor limite, el numero de infectados desciende en el tiempo, y si es mayor que el valor limite "S", asciende. 
{{matlab|codigo=
%Apartado 3 trabajo 6-C 2015
clear all
t0=0; tN=40; %Intervalo de tiempo tomado
S0=input('introduce valor de población susceptible: '); 
I0=2000; %Valores iniciales de infectados(I) y susceptibles(S) 
h=0.1; %Determinación del paso 
t=t0:h:tN; %Desarrollo del tiempo desde t0 hasta tN tomando intervalos de paso h
a=0.003; %valores de los parametros
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0; %Asignacion del valor incial para la primera componente
I(1)=I0;
 
%Resolucion del sistema de forma matricial
for n=1:N
   A=[S0;I(n)]+h*[-a*S0*I(n);a*S0*I(n)-(b+c)*I(n)];
   S(n+1)=S0;%Asignacion de los distintos valores de S como la primera compnente de la matriz A
   I(n+1)=A(2); %Asignacion de los distintos valores de I como la primera componente de la matriz A
end

%Interpretación gráfica 
hold on
figure(1) 
plot(t,S,'r')
plot(t,I,'b')
hold off}}

[[Archivo:s100.jpg|520px|thumb|left|Gráfica para 100 supceptibles]]   [[Archivo:s200.jpg|520px|thumb|right|Gráfica para 200 supceptibles]]

Podemos llegar a la conclusión de que el número de personas infectadas varía notablemente cuando el número de susceptibles se mantiene constante a lo largo del tiempo. Esto se debe a que cuanto mayor sea el número de susceptibles, mayor será el número de infectados, es decir, con 200 susceptibles, mantiene un ritmo constante hasta que llega un momento en el cual el número de infectados se dispara, debido a que la tasa de personas susceptibles es mayor que las personas curadas o fallecidas. En conclusión, debe de haber un valor de personas susceptibles para que la función de infectados sea lineal, es decir, sea constante en el tiempo.




==Modelo completo==

{{matlab|codigo=
%Apartado 4 trabajo 6-C 2015
clear all
t0=0; tN=40; %Intervalo de tiempo tomado
S0=input('introduce valor de población susceptible: '); 
I0=input('introduce valor de población infectada inicial: '); %Valores iniciales de infectados(I) y susceptibles(S) 
h=input('introduce valor del paso de tiempo: '); %Determinación del paso 
t=t0:h:tN; %Desarrollo del tiempo desde t0 hasta tN tomando intervalos de paso h
a=0.003; %valores de los parametros
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0; %Asignacion del valor incial para la primera componente
I(1)=I0;
 
%Resolucion del sistema de forma matricial
for n=1:N
   A=[S(n);I(n)]+h*[-a*S(n)*I(n);a*S(n)*I(n)-(b+c)*I(n)];
   S(n+1)=A(1);%Asignacion de los distintos valores de S como la primera compnente de la matriz A
   I(n+1)=A(2); %Asignacion de los distintos valores de I como la primera componente de la matriz A
end

%Interpretación gráfica 
hold on
figure(1) 
plot(t,S,'r')
plot(t,I,'b')
hold off
maximo=max(I)
posicion=find(I==maximo)
tiempo=t(posicion)}}

Como podemos apreciar en las gráficas el número máximo de enfermos esperados con 800 susceptibles y 20 infectados inicialmente, es aproximadamente 500, mientras que para 10000 susceptibles y 40 infectados inicialmente, observamos un aumento de personas enfermas, aproximadamente de 9500.



[[Archivo:411.jpg|360px|thumb|left]] [[Archivo:412.jpg|360px|thumb|right]]
[[Archivo:413.jpg|360px|thumb|left]] [[Archivo:414.jpg|360px|thumb|right]]
[[Archivo:421.jpg|360px|thumb|left]] [[Archivo:422.jpg|360px|thumb|right]]
[[Archivo:423.jpg|360px|thumb|left]] [[Archivo:424.jpg|360px|thumb|right]]






































































































==Comparación Runge-Kutta con Euler==
{{matlab|codigo=
%Apartado 5 trabajo 6-C 2015
clear all
t0=0; tN=30; %Intervalo de tiempo tomado
y1=input('introduce valor de población susceptible inicial: ');
y2=input('introduce valor de población infectada inicial: ');
y0=[y1;y2]; %Valores de población iniciales (Susceptibles e Infectados)
h=input('introduce valor del paso del tiempo: ');
N=(tN-t0)/h;
y=y0;
S(1)=y(1); %Asignación del valor incial para la primera componente de S
I(1)=y(2); %Asignación del valor incial para la primera componente de I
a=0.003; % Valores de los parámetros
b=0.3; c=0.01; %Resolución empleando el método Runge Kutta
 
for n=1:N  
    k1=[-a*y(1)*y(2);a*y(1)*y(2)-(b+c)*y(2)];    
    k2=[-a*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2));a*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2)-(b+c)*(y(2)+1/2*h*k1(2)))];         
    k3=[-a*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2));a*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2)-(b+c)*(y(2)+1/2*h*k2(2)))];
    k4=[-a*(y(1)+h*k3(1))*(y(2)+h*k3(2));a*(y(1)+h*k3(1))*(y(2)+h*k3(2)-(b+c)*(y(2)+h*k3(2)))];    
    y=y+h/6*(k1+2*k2+2*k3+k4);    
    S(n+1)=y(1);    
    I(n+1)=y(2);
end
t=t0:h:tN;
figure(1)
hold on
plot(t,S,'b')
plot(t,I,'r')
hold off

for n=1:N
    A=[S(n);I(n)]+h*[-a*S(n)*I(n);a*S(n)*I(n)-(b+c)*I(n)];
    S(n+1)=A(1);%Asignacion de los distintos valores de S como la primera compnente de la matriz A
    I(n+1)=A(2); %Asignacion de los distintos valores de I como la primera componente de la matriz A
end
%Interpretación gráfica
figure(2)
hold on
plot(t,S,'b')
plot(t,I,'r')
hold off}}
[[Archivo:runge.jpg|left|360px|thumb|Texto de la leyenda]]
[[Archivo:runge22.jpg|rigth|360px|thumb|Texto de la leyenda]]



























==Tasa de infectados por contagio dependiente del tiempo. Método de Heun==

Aquí se suele poner algo...

{{matlab|codigo=

clear all
%Datos
t0 = 0;
tN = 40;
h = 0.1;
N = round((tN-t0)/h); 
t = t0:h:tN;

S = zeros(1,N+1);
I = zeros(1,N+1);
S(1) = 1600;
I(1) = 40;
%Valores de los parámetros

a = zeros(1,N+1);
a(1) = 0.003/(1+t0);
b = 0.3;
c = 0.01;
 
%Calcular los vectores infectados y susceptibles.

for i = 1:N
    a(i+1) = 0.003/(1+t(i));
 
    K1 = -a(i)*S(i)*I(i);  
    K2 = -a(i)*S(i)*I(i) + K1*h;
    S(i+1) = S(i) + h/2*(K1+K2); 
    Z1 = (a(i)*S(i)*I(i)) -I(i)*(b+c);
    Z2 = (a(i)*S(i)*I(i)) -I(i)*(b+c) + Z1*h;
    I(i+1) = I(i) + h/2*(Z1+Z2);
 
end
 
%Gráficas
hold on
plot (t,S,'-')
plot (t,I,'-r')
hold off
}}


==Calibración del parámetro a==


{{matlab|codigo=
clear all
%Datos

t0 = 0;
tN = 8;
h = 0.0001;

N = (tN-t0)/h;
t = t0:h:tN;

S = zeros(1,N+1);
I = zeros(1,N+1);

S(1) = 1600;
I(1) = 40;


%Valores de los parámetros

a = 0.0005:h:0.0020;
b = 0.3;
c = 0.01;

%Definimos los vectores M, T y D, donde guardaremos respectivamente
%los máximos de infectados, el tiempo en el que se producen y la 
%diferencia entre ese tiempo y los 5 días  

M=zeros(1,length(a));
T=zeros(1,length(a));
D=zeros(1,length(a));

%Calcular los vectores infectados y susceptibles

for j=1:length(a)
A=a(j);
  for i = 1:N
    K1 = -A*S(i)*I(i);  
    K2 = -A*S(i)*I(i) + K1*h;
    S(i+1) = S(i) + h/2*(K1+K2); 
    Z1 = (A*S(i)*I(i)) -I(i)*(b+c);
    Z2 = (A*S(i)*I(i)) -I(i)*(b+c) + Z1*h; 
    I(i+1) = I(i) + h/2*(Z1+Z2);
  end
M(j)=max(I)
posicion=find(I==M(j));
tiempo=t(posicion);
T(j)=tiempo
D(j)=abs(T(j)-5);
end
[M']
[D']
[T']
T1=min(D);
posicion=find(D==T1);

%soluciones obtenidas

A_buscado=a(posicion)
Maximo_en_ese_caso=M(posicion)


}}

Cuando calibramos el coeficiente de a con una experiencia anterior en el que el máximo de personas infectadas se alcanza a los 5 días obtenemos el resultado de que a=0.0008 es el valor que mas se ajusta a esa experiencia, tomando como paso h=10^-4.



[[Categoría:Ecuaciones Diferentiales]] 
[[Categoría:ED14/15]]
[[Categoría:Trabajos 2014-15]]