Diferencia entre revisiones de «T4. Modelos Epidemiológicos»

Revisión del 19:05 6 mar 2015

1 Interpretación de los diferentes parámetros

En el desarrollo de la epidemia encontramos dos tipos distintos de población que influyen a la evolución de la misma y que se rigen por las siguientes ecuaciones

[math]\frac{\mathrm{d} S}{\mathrm{d} T}=-aSI [/math]

[math]\frac{\mathrm{d} I}{\mathrm{d} T}=aSI-bI-cI [/math]

Identificando primero las variables y luego los parámetros diferenciamos:

• T: tiempo
• S: población que puede ser infectada por la enfermedad. Parámetro en función del tiempo.
• I: población infectada por la enfermedad. Parámetro en función del tiempo.
• a: interacción entre el número de individuos infectados y no infectados.
• b: fallecimientos a causa de la enfermedad.
• c: población que consigue recuperarse de la enfermedad.

Además, en el desarrollo de la epidemia consideraremos que la tasa de cambio depende tanto de los fallecimientos como de la población curada y que la tasa de población susceptible a ser infectada es proporcional a la interacción entre ambas poblaciones.

Si queremos entender las ecuaciones por las que se rige nuestra epidemia:

[math]\frac{\mathrm{d} S}{\mathrm{d} T}=-aSI [/math]

Representa la exposición a la que se encuentra la población frente a la epidemia, es decir, la variación de la población susceptible a ser infectada, [math]S[/math], a causa de la interacción directa de ésta con la población ya infectada, [math]I[/math], en función del paramtro [math]a[/math] antes descrito. Así observamos que es una función no lineal (producto de las variables [math]S[/math] y [math]I[/math]) y decreciente en el tiempo.

[math]\frac{\mathrm{d} I}{\mathrm{d} T}=aSI-bI-cI [/math]

En este caso se representa la variación de la población infectada con un primer término similar a la primera ecuación pero de signo contrario lo que significa que este término hace crecer el número de infectados por la epidemia y otro posterior ya únicamente dependiente de la población infectada, [math]I[/math], que esta acompañado por [math]-(b+c)[/math] que representas las muertes y curas a causa de la enfermedad.

2 Método de Euler y Trapecio

Consideramos que los parámetros del problema en el desarrollo de los metodos de Euler y del Trapecio son [math]I0 = 2000[/math], [math]b = 0.3[/math] y [math]c = 0.01[/math] y lo resolveremos en ambos casos para una discretización [math]h = 0.1[/math]

2.1 S=0

%Apartado2
clear all
t0=0; tN=40; 
I0=2000; 
h=0.01; 
N=(tN-t0)/h;
t=t0:h:tN;
I=zeros(1,N+1);
I(1)=I0;
for i=1:N
    I(i+1)=I(i)-h*(0.3+0.001)*I(i);
end
%Gráficas
[t',I'];
[minimo,indice]=min(abs(I-500))
f=I(indice)
p=t(find(I==(500-minimo)))
hold on
plot(p,I,'r')
plot(t,f,'b')
plot(t,I)
hold off
sprintf('Los 500 enfermos se alcanza al cabo de %d días', p)

Evolución de la población de los Infectados con ausencia de población susceptible a ser infectada. Se observa claramente la desaparición de la población al completo con el paso del tiempo.

La población de infectados se reduce hasta un cuarto del total inicial, siendo ese valor igual a [math]500[/math], en un tiempo en torno a [math]5[/math] segundos, información sacada directamente de la interpretación de la gráfica aportada. [math]5[/math]

2.2 S=100

En este segundo caso necesitamos también el valor del parametro [math]a=0.003[/math].

%Apartado 3
clear all
t0=0; tN=40; %Intervalo t
S0=100; 
I0=2000;
%Valores iniciales I y S 
h=0.1; 
t=t0:h:tN; 
a=0.003; %Dato problema
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0; 
I(1)=I0;
%Resolucion 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;
    I(n+1)=A(2); 
end 
%Gráfica 
hold on
figure(1) 
plot(t,S,'r')
plot(t,I,'b')
hold off

Evolución de la población de los Infectados en presencia de población susceptible a ser infectada (S=100). Se observa que la población susceptible se mantiene constante en el tiempo mientras que la población de infectados desciende linealmente en el tiempo.

Evolución de la población de los Infectados en presencia de población susceptible a ser infectada (S=200). En este caso se observa que la población susceptible acaba por desaparecer mientras que la poblacion infectada crece de forma exponencial.

3 Euler Completo

En el modelo de Euler completo consideramos ahora los datos iniciales [math]S0=800[/math], [math]I0=20[/math] y [math]S0=10000[/math] y [math]I0=40[/math] para tiempos menores de [math]40[/math] días y discretización temporal [math]h=0.1, h=0.01,h=0.001,h=0.0001[/math].

%Apartado 4
clear all
t0=0; tN=40; %Intervalo t
S0=input('N población susceptible: '); 
I0=input('N infectados: ');
%Valores iniciales I y S 
h=input('Valor del paso: '); 
t=t0:h:tN; 
a=0.003; 
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0; 
I(1)=I0;
%Resolucion 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);
    I(n+1)=A(2); 
end 
%Gráfica 
hold on
figure(1) 
plot(t,S,'r')
plot(t,I,'b')
hold off
maxi=max(I);
pos=find(I==maxi);
tiempo=t(pos);
t1=sprintf('El número máximo de infectados será %d',maxi);
t2=sprintf(' y se dará transcurridos %d días',tiempo);
t3=[t1 t2];
disp(t3)

3.1 S0=800, I0=20

h=0.1. El máximo de enfermos esperado en este caso es de [math]517.44[/math] y se va a producir transcurridos [math]2[/math] días y [math]19[/math] horas.

h=0.01. El máximo de enfermos esperado en este caso es de [math]506.3[/math] y se va a producir transcurridos [math]2[/math] días y [math]17[/math] horas.

h=0.001.El máximo de enfermos esperado en este caso es de [math]505.3[/math] y se va a producir transcurridos [math]2[/math] días y [math]17[/math] horas.

h=0.0001. El máximo de enfermos esperado en este caso es de [math]505.2[/math] y se va a producir transcurridos [math]2[/math] días y [math]17[/math] horas.

Por tanto en esta primera simulación de poblaciones observamos que cuando reducimos el paso [math]h[/math] obtenemos mayor precisión en la población así pues podemos indicar que el máximo alcanzado será entonces de [math]505[/math], obtenido a los [math]2[/math] días y [math]17[/math] horas aproximadamente.

3.2 S0=10000 y I0=40

h=0.1. El máximo de enfermos esperado en este caso es de [math]12680.8[/math] y se va a producir transcurridas [math]12[/math] horas.

h=0.01. El máximo de enfermos esperado en este caso es de [math]9521.1[/math] y se va a producir transcurridas [math]8[/math] horas y [math]24[/math] minutos.

h=0.001. El máximo de enfermos esperado en este caso es de [math]9469.6[/math] y se va a producir transcurridas [math]8[/math] horas y [math]12[/math] minutos.

h=0.0001.El máximo de enfermos esperado en este caso es de [math]9464.7[/math] y se va a producir transcurridas [math]8[/math] horas y [math]11[/math] minutos.

El máximo de enfermos esperado será en esta ocasión de [math]9465[/math] y se va a producir al cabo de [math]8[/math] horas y [math]11[/math] minutos.

4 Runge-Kutta 4º Orden

4.1 Comparativa con Euler

Si realizamos una comparativa entre ambos métodos, encontramos la similitud de que ambos métodos están basados en iteraciones. El método de Euler busca la pendiente de la recta tangente empleando un paso muy pequeño. Por su parte, el método de Runge Kutta emplea cuatro parámetros, que son las pendientes de las cuatro tangentes a la curva integral, y empleando un paso mayor. Podemos concluir, por tanto, que el método de Runge Kutta supone una mejora con respecto al de Euler.

4.2 Dificultad Método del Trapecio

Cuando se emplea un método implícito, como podría ser el método trapezoidal, para resolver un problema, la mayor dificultad que podemos encontrar es el hecho de que deberemos despejar de una forma manual aquellos parámetros necesarios para poder posteriormente programar la solución. Puede que, en algún caso, el proceso sea bastante complicado, llegando incluso a ser imposible despejar.

%Apartado5
clear all
t0=0; tN=30; %Intervalo t
y1=input('N población susceptible: ');
y2=input('N infectados de inicio: ');
y0=[y1;y2]; %Valores iniciales S e I
h=input('valor del paso: ');
N=(tN-t0)/h;
t=t0:h:tN;
y=y0;
S(1)=y(1); 
I(1)=y(2); 
%Parámetros
a=0.003; 
b=0.3; 
c=0.01;
%Euler
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);
    I(n+1)=A(2); 
end
%Gráfica
figure(1)
hold on
plot(t,S,'g')
plot(t,I,'r')
legend('Susceptible','Infectados','Location','best')
hold off
%R-k
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);
    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);
end
hold on
figure(2)
plot(t,S,'g')
plot(t,I,'r')
legend('Susceptible','Infectados','Location','best')
hold off

h=0.1. Gráfica de euler con S0=800 y I0=20 con el paso mayor.

h=0.1. Gráfica de Runge-Kutta 4º con S0=800 y I0=20 con el paso mayor.

h=0.0001. Gráfica de euler con S0=800 y I0=20 con el paso menor.

h=0.0001.Gráfica de Runge-Kutta 4º con S0=800 y I0=20 con el paso menor.

GRAFICAS FALTA SACAR

5 Método Heun, Variable a dependiente de t, a(t)

%Apartado6
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; 
a=zeros(1,N+1);
a(1)=0.003/(1+t0); %a(t)
b=0.3; c=0.01;
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);
    l1=a(i)*S(i)*I(i)-I(i)*(b+c);
    l2 =(a(i)*S(i)*I(i))-I(i)*(b+c)+l1*h;
    I(i+1)=I(i)+h/2*(l1+l2);
end
%Gráficas
hold on
plot (t,S,'b')
plot (t,I,'r')
hold off

El coeficiente de variable a(t) lo interpretamos.....