Diferencia entre revisiones de «Modelos epidemiológicos. Grupo C14»

Revisión del 10:49 6 mar 2015

En esta página se discutirá y analizará un modelo epidemiológico mediante métodos de cálculo numérico.

El modelo se basa en las siguientes hipótesis:

- El número de personas infectadas sólo se altera por el fallecimiento o cura de éstas, y se ve afectada por el número de contagios entre la población de riesgo.

- La tasa de individuos que pasan de la población de riesgo a estar infectados es proporcional a la interacción entre el número de individuos entre ambas clases.

1 . INTERPRETACIÓN DEL MODELO

Matemáticamente el modelo viene definido por las siguientes ecuaciones:

[math] \left\{\begin{matrix} \frac{dS}{dt}=-aSI\\ \frac{dI}{dt}=aSI-bI-cI \end{matrix}\right. [/math]

Donde:

-La función [math] S [/math] indica la cantidad de personas en riesgo de ser infectadas e [math] I [/math] el número de individuos que ya han contraído la enfermedad.

- El valor de [math] a [/math] es un coeficiente relacionado con la probabilidad de contagio cuando se juntan individuos de la población de riesgo con la de infectados. El coeficiente [math] b [/math] es el porcentaje de infectados que sanan en una unidad de tiempo, en este caso días, y el [math] c [/math] el porcentaje de infectados que fallece cada día.

2 . CASO DE POBLACIÓN DE RIESGO CONSTANTE

Conocido el dato de población de riesgo y siendo éste constante,el problema queda simplificado a una única ecuación diferencial lineal y de coeficientes constantes.

Para resolver dicha ecuación, se han utilizado los métodos de Euler y trapecio, dados los siguientes datos iniciales:

[math] \left\{\begin{matrix} I_{0}=2000\\ a=0.003\\ b=0.3\\ c=0.01 \end{matrix}\right.[/math]

Éste es el programa que resuelve la ecuación mediante el método de Euler, método del trapecio y a su vez calcula cuándo llegará el número de infectados a la cuarta parte por este último método, puesto que es el más preciso:

2.1 . [math] S=0 [/math]

Debido a que la ecuación lineal es de coeficientes constantes, es sencillo obtener la solución de una manera analítica, resultando:

[math] I=2000e^{-0.31t}[/math]

Para [math] S=0 [/math] el tiempo que tarda en reducirse el número de infectados a la cuarta parte es: [math]t=4,6 días[/math]

Las gráficas devueltas, son:

Gráfica para S constante de valor 0

2.2 . [math] S=100 [/math]

Resuelto analíticamente:
[math] I=2000e^{-0.01t}[/math]

Para [math] S=100 [/math] el tiempo en reducirse es: [math]t=40,1 días[/math]

Las gráficas devueltas, son:

Gráfica para S constante de valor 100

2.3 . [math] S=200 [/math]

Al resolver la ecuación de manera analítica, se obtiene:
[math] I=2000e^{0.29t}[/math]

La gráfica resultante es:
imagen

En la gráfica se puede observar que el número de personas infectadas crecerá exponencialmente a lo largo del tiempo.

Analíticamente se demuestra que ésto ocurrirá a partir de [math] S\geq 104 [/math].
[math]\frac{\mathrm{d} I}{\mathrm{d} t}= aSI-bI-cI[/math]
Sustituyendo en dicha ecuación los valores iniciales,
[math]\frac{\mathrm{d} I}{\mathrm{d} t}= 0.003SI-0.3I-0.01I[/math]
La ecuación anterior alcanza su punto crítico en [math] S=104 [/math] .

Por tanto, para [math] S=200 [/math] no se puede obtener el tiempo que tardará el número de infectados en reducirse a la cuarta parte, puesto que la función es creciente.

3 . ANÁLISIS DEL SISTEMA POR EL MÉTODO DE EULER

Para este análisis, se considera el modelo completo, teniendo en cuenta ambas ecuaciones del sistema y utilizando los siguientes datos:

[math]\left\{\begin{matrix} a=0.003\\ b=0.3\\ c=0.01\\ t\epsilon [0,40] \end{matrix}\right.[/math]

Además, se han analizado distintas soluciones a partir de la variación del paso:
[math]h=10^{-1}, h=10^{-2}, h=10^{-3}, h=10^{-4}[/math]

Éste ha sido el programa utilizado:
Insertar programa

3.1 . [math] (S_{0},I_{0})=(800,20) [/math]

Insertar gráficas para las distintas h

El programa además, devuelve el valor máximo de enfermos esperados [math]I_{max}=?[/math] y el momento en el que ésto ocurrirá [math]t=?[/math]. Se tienen en cuenta los datos obtenidos por el programa con el menor paso, ya que serán los más exactos.

3.2 . [math] (S_{0},I_{0})=(10000,40) [/math]

Insertar gráficas para las distintas h

Se puede observar, que a partir del paso [math]h=?[/math], el resultado obtenido por la gráfica es incongruente, debido a que el valor del número de personas en riesgo desciende rápidamente.

Ésto ocurre porque el número inicial de personas infectadas [math]I_{0}[/math] es muy elevado y al aumentar el valor del paso, el programa se hace cada vez mas inexacto, provocando que la derivada en el tramo sea muy grande, influyendo por tanto en la pendiente de la gráfica.

Utilizando el menor paso, para que el programa sea lo más exacto posible, se obtiene el número máximo de enfermos esperado, que es [math]I_{max}=?[/math] y cuándo se producirá [math]t=?[/math]. Al haber utilizado [math]h=10^{-4}[/math], el error mencionado anteriormente no influye.

4 . ANÁLISIS DEL SISTEMA POR RUNGE-KUTTA

A continuación se resuelve el sistema por el método de Runge-Kutta, de cuarto orden, partiendo de los mismos datos que los utilizados para el método de Euler:

[math]\left\{\begin{matrix} a=0.003\\ b=0.3\\ c=0.01\\ t\epsilon [0,40] \end{matrix}\right.[/math]

El programa por el que se obtiene la solución, a partir de los datos iniciales [math](S_{0},I_{0})[/math], es el siguiente:

clear all
t0=0;tN=20;
h=0.1;
N=(tN-t0)/h;
t=t0:h:tN; 
%VALORES
a=0.003; b=0.3; c=0.01;
%x0=[800,20]; 
S0=800;
I0=20;
S=zeros(1,N+1);
I=S;
S(1)=S0;
I(1)=I0;

x=zeros(2,N+1); %Es una matriz de dos filas y N+1 columnas
%Inicializamos el bucle
%x(:,1)=x0;
%xx=x0; 
for i=1:N
    k11=-a*S(i)*I(i);
    k12=a*S(i)*I(i)-b*I(i)-c*I(i);
    k21=-a*(S(i)+(h/2)*k11)*(I(i)+(h/2)*k12);
    k22=a*(S(i)+(h/2)*k11)*(I(i)+(h/2)*k12)-b*(I(i)+(h/2)*k12)-c*(I(i)+(h/2)*k12);
    k31=-a*(S(i)+(h/2)*k21)*(I(i)+(h/2)*k22);
    k32=a*(S(i)+(h/2)*k21)*(I(i)+(h/2)*k22)-b*(I(i)+(h/2)*k22)-c*(I(i)+(h/2)*k22);
    k41=-a*(S(i)+h*k31)*(I(i)+h*k32);
    k42=a*(S(i)+h*k31)*(I(i)+h*k32)-b*(I(i)+h*k32)-c*(I(i)+h*k32);
 
 S(i+1)=S(i)+(h/6)*(k11+2*k21+2*k31+k41);
 
 I(i+1)=I(i)+(h/6)*(k12+2*k22+2*k32+k42);
  
end
 hold on
 plot(t,S)
plot(t,I,'r')
legend('S','I')

Obteniendo las siguientes gráficas: imagen para (S0,I0)=(800,20) y para (S0,I0)=(10000,40)

Este método comparado con Euler, da lugar a las gráficas:
imagen para los distintos (S0,I0)

Usar el método trapezoidal, o cualquier otro método implícito, conlleva un excesivo trabajo operativo producido por su alta complejidad y duración. Ésto deriva en que la utilización de métodos como Euler o Runge-Kutta sea predominante a la hora de resolver este tipo de sistemas.

5 . ANÁLISIS DEL SISTEMA CON PARÁMETRO [math] a [/math] VARIABLE

En este caso, se tiene que el parámetro [math] a [/math], es una función que depende del tiempo:
[math] a(t)=\frac{0.003}{1+t} [/math]

Además se tiene unos datos iniciales de infectados [math] (I_{0}) [/math] y de población total [math] (S+I) [/math] de:
[math]\left\{\begin{matrix} I_{0}=40\\ S+I=1640 \end{matrix}\right.[/math]

Tomando dichos datos, se ha utilizado el método de Heun para aproximar la solución y así poder dibujar sendas gráficas, susceptibles e infectados: imagenes

clear all; clf

t0=0; tN=20; S0=1600; I0=40; h=0.1; %a=0.003; 
b=0.3; c=0.01;

N=(tN-t0)/h;
t=t0:h:tN; %vector n+1 componentes
S=zeros(1,N+1);
a=S;
S(1)=S0;
I=S;
I(1)=I0;


for i=1:N
    a(i)=0.003/(1+t(i));
    Ss=S(i)-h*a(i)*I(i)*S(i);
    Ii=I(i)+h*(a(i)*I(i)*S(i)-b*I(i)-c*I(i)); %aproximaciones intermedias de Heun
    
    S(i+1)=S(i)+h/2*(-a(i)*S(i)*I(i)-a(i)*Ii*Ss);
    
    I(i+1)=I(i)+h/2*(a(i)*S(i)*I(i)-b*I(i)-c*I(i)+a(i)*Ss*Ii-b*Ii-c*Ii);   
    
end
[fila,col]=find(I==max(I)); %hallar la posicion del maximo dentro del vector.
plot(t,S,'r')
hold on
plot(t,I)
legend('Susceptibles','Infectados')

El coeficiente [math] a(t)=\frac{0.003}{1+t} [/math] , varía con respecto al tiempo, tal y como indica la fórmula.
Se puede observar que a medida que pasa el tiempo el coeficiente disminuye, esto puede ser debido a la inmunidad, que se basa en aumentar la resistencia frente a una enfermedad infecciosa y es directamente proporcional al tiempo. Durante la primera fase de una infección el crecimiento es muy rápido, porque el sistema inmune ni conoce ni sabe combatir esa infección. Una vez alcanzado el máximo, el cuerpo humano ya ha identificado el anticuerpo y sabe como combatirlo debido a que el sistema inmune humano tiene memoria. Esto provoca que esa infección empezará a descender rápidamente, llegando incluso hasta la desaparición de la infección.

6 . CALIBRACIÓN DEL COEFICIENTE [math] a [/math]

Para calibrar el coeficiente [math] a [/math], se tomará una situación similar a la del apartado anterior, tomando [math] a [/math] constante en el intervalo [math]a\in [0.0005, 0.002][/math]. A partir de una experiencia previa se conoce que el máximo de personas infectadas se alcanza a los 5 días.
[math]\left\{\begin{matrix} I_{0}=40\\ S+I=1640\\ a(t)=\frac{0.003}{1+t} \end{matrix}\right.[/math]

El programa utilizado para resolver este apartado requiere que se divida el intervalo al que pertenece [math] a [/math] en un conjunto de valores equidistribuidos, para ello se ha tomado un paso de [math] h=10^{-4} [/math]
Se usará un bucle para calcular el máximo para cada valor de [math] a [/math], y se seleccionará aquel que esté más cerca de 5 días, que es cuando el número de personas infectadas alcanza su punto más alto.

ha=0.0001;
a0=0.0005;
aN=0.002;
a=a0:ha:aN;
Na=round((aN-a0)/ha);
ht=0.01;
t0=0;
tN=10;
t=t0:ht:tN;
Nt=round((tN-t0)/ht);
S0=1600;
I0=40;
b=0.3;
c=0.01;
S=zeros(1,Nt+1);
I=zeros(1,Nt+1);
maximos=zeros(1,Na);
I(1)=I0;
S(1)=S0;
for i=1:Na
aaux=a(i);
  for j=1:Nt
  I(j+1)=I(j)+ht*((aaux*S(j)*I(j))-b*I(j)-c*I(j));
  S(j+1)=S(j)+ht*(-aaux*S(j)*I(j));
  end
[maximo,posicion]=max(I(1,:));
tiempohastamax=posicion*ht;
maximos(1,i)=tiempohastamax;
end
maximos;
difemaximos=zeros(size(maximos));
for j=1:15;
    difemaximos(1,j)=abs(maximos(j)-5);
end
[minimo posicionminimo]=min(difemaximos(1,:));
posicionminimo;
aparamaximoen5dias=a(posicionminimo)

imagen