T4. Modelos Epidemiológicos
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Modelos Epidemiológicos, desarrollo de una epidemia. Grupo 28 |
| Asignatura | Ecuaciones Diferenciales |
| Curso | Curso 2014-15 |
| Autores | Luis Doreste Henríquez 127 Luis López Díaz 239 Enrique Martínez Mur 271 |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
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)
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
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
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
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
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.....
