Ver el código de «Modelos epidemiológicos. Grupo C14»

{{ TrabajoED | Modelos epidemiológicos. Grupo C14 | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | Palacios Pintor, Pedro (1188)
 
Pontiveros Bermejo, Diego (1197)

Reinoso Muñoz, Cristina (988)

Rojas Arranz, Almudena (942)

Torre Prado, Yago de la (1006)

Vidal Sánchez, Nieves (1236)}}

[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED14/15]]

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.

==. Interpretación del modelo==

Matemáticamente el modelo viene definido por las siguientes ecuaciones:<br />

<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 un término relacionado con el porcentaje de infectados que sanan en una unidad de tiempo, en este caso días, y análogamente, <math> c </math> guardará una relación con  el porcentaje de infectados que fallece cada día.

==. 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.<br/>

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

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


É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:<br/>
{{matlab| codigo=
%calculos numericos
clf
t0=0;
tN=40;
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
saux=ones(1,N+1);
s=input('Introduzca el valor de S: ');
figure(1)
hold on
%solucion por euler
yeuler=zeros(1,N+1);
%condición inicial
yeuler(1)=2000;
%bucle
for i=1:N
yeuler(i+1)=yeuler(i)+h*(0.003*s*yeuler(i)-0.31*yeuler(i));
end
clear ('i');
plot(t,yeuler,'r')
%trapecio
ytrape=zeros(1,N+1);
%condicion inicial
ytrape(1)=2000;
%bucle
for i=1:N
ytrape(i+1)=(ytrape(i)+(0.0015*s*ytrape(i)-0.155*ytrape(i))*h)/(1-0.0015*s*h+0.155*h);
end
clear ('i');
plot(t,ytrape,'k');
legend ('solucion euler','solucion por trapecio','location','best');
%el numero de infectados se habra reducido a la cuarta parte cuando sea igual a 500
%para este apartado usamos la aproximación del trapecio ya que es la mas precisa de las dos numericas calculadas.
hold off
figure(2);
plot(t,ytrape);
hold on
cuartaparteaux=ones(1,N+1);
cuartaparte=500*cuartaparteaux;
plot(t,cuartaparte,'r');
legend('solucion por trapecio','cuarta parte de infectados','location','best');
xlabel('t(días)')
ylabel('Población')
%calculamos numericamente la distancia a la recta
hold off
distancia=zeros(1,N+1);
for i=1:(N+1)
distancia(i)=abs((ytrape(i)-500));
end
%buscamos la posicion del minimo 
[minimo posicion]=min(distancia(1,:));
%posicion nos dice el valor de t para el que se da el minimo, lo transformamos a tiempo
tiempoacuarta=h*(posicion-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:<br/>

<math> I=2000e^{-0.31t}</math><br/>


El tiempo que tarda en reducirse el número
de infectados a la cuarta parte es:<br/> <math>t=4,6</math> días.<br/>


Las gráficas devueltas, son:<br/>

{|
|-
| [[Archivo:EulerS000.jpg|thumb|300px|left|Gráfica para S constante de valor 0]] || [[Archivo:EulerS0zoom.jpg|thumb|300px|left|Gráfica para S constante de valor 0 ampliada]] || [[Archivo:EulerS0cuarta.jpg|thumb|300px|left|Gráfica para la reducción de infectados a la cuarta parte]]
|}

En estas gráficas, se puede observar que el número de infectados descenderá de manera rápida, puesto que en ningún momento habrá personas susceptibles a contraer la enfermedad; lo que significa que los infectados sanarán o fallecerán, sin enfermar otras personas.

===. <math> S=100 </math>===
Resuelto analíticamente:<br/>
<math> I=2000e^{-0.01t}</math><br/>


El tiempo en reducirse es:<br/> <math>t=138.6</math> días.<br/>


Las gráficas devueltas, son:<br/>

{|
|-
| [[Archivo:EulerS0100.jpg|thumb|300px|left|Gráfica para S constante de valor 100]] || [[Archivo:EulerS100zoom.jpg|thumb|300px|left|Gráfica para S constante de valor 100 ampliada]] || [[Archivo:EulerS100cuarta.jpg|thumb|300px|left|Gráfica para la reducción de infectados a la cuarta parte]]
|}

Como se puede observar en las gráficas, el número de infectados descenderá más lentamente que en el caso anterior (<math>S=0</math>), puesto que aquí sí que habrá personas susceptibles a contraer la enfermedad, que se convertirán en infectadas. Ésto también queda reflejado en el tiempo que tarda en reducirse el número de personas infectadas a la cuarta parte, que en este caso es mucho mayor.

===. <math> S=200 </math>===
Al resolver la ecuación de manera analítica, se obtiene:<br/>
<math> I=2000e^{0.29t}</math> <br/>


Las gráficas resultantes son:<br/>
{|
|-
| [[Archivo:EulerS200.jpg|thumb|300px|left|Gráfica para S constante de valor 200]] || [[Archivo:EulerS200zoom.jpg|thumb|300px|left|Gráfica para S constante de valor 200 ampliada]]
|}


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


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


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.

==. 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:<br/>

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


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


Éste ha sido el programa utilizado:<br/>
[[Archivo:diego1.jpg|400px|thumb||right|Comparación entre las gráficas obtenidas para la población de enfermos para cada paso con (S0,I0)=(800,20)]]
[[Archivo:diego2.jpg|400px|thumb||right|Comparación entre las gráficas obtenidas para la población de susceptible para cada paso con (S0,I0)=(800,20)]]

{{matlab|codigo=
h1=0.1;
h2=0.01;
h3=0.001;
h4=0.0001;
t0=0;
tN=40;
s0=800;
I0=20;
N1=round((tN-t0)/h1);
N2=round((tN-t0)/h2);
N3=round((tN-t0)/h3);
N4=round((tN-t0)/h4);
t1=t0:h1:tN;
t2=t0:h2:tN;
t3=t0:h3:tN;
t4=t0:h4:tN;
seuler1=zeros(1,N1+1);
ieuler1=zeros(1,N1+1);
seuler2=zeros(1,N2+1);
ieuler2=zeros(1,N2+1);
seuler3=zeros(1,N3+1);
ieuler3=zeros(1,N3+1);
seuler4=zeros(1,N4+1);
ieuler4=zeros(1,N4+1);
%imponemos la condicion inicial
seuler1(1,1)=800;
seuler2(1,1)=800;
seuler3(1,1)=800;
seuler4(1,1)=800;
ieuler1(1,1)=20;
ieuler2(1,1)=20;
ieuler3(1,1)=20;
ieuler4(1,1)=20;
%construimos los bucles que nos den los resultados
for i=1:N1
seuler1(i+1)=seuler1(i)+h1*(-0.003*seuler1(i)*ieuler1(i));
ieuler1(i+1)=ieuler1(i)+h1*((0.003*seuler1(i)*ieuler1(i))-0.31*ieuler1(i));
end
figure(1);
hold on
plot(t1,seuler1,'k')
plot(t1,ieuler1,'g')
legend('poblacion de riesgo','enfermos','location','best');
hold off
clear('i');
for i=1:N2
seuler2(i+1)=seuler2(i)+h2*(-0.003*seuler2(i)*ieuler2(i));
ieuler2(i+1)=ieuler2(i)+h2*((0.003*seuler2(i)*ieuler2(i))-0.31*ieuler2(i));
end
figure(2);
hold on
plot(t2,seuler2)
plot(t2,ieuler2,'r')
legend('poblacion de riesgo','enfermos','location','best');
hold off
clear('i');
for i=1:N3
seuler3(i+1)=seuler3(i)+h3*(-0.003*seuler3(i)*ieuler3(i));
ieuler3(i+1)=ieuler3(i)+h3*((0.003*seuler3(i)*ieuler3(i))-0.31*ieuler3(i));
end
figure(3);
hold on
plot(t3,seuler3)
plot(t3,ieuler3,'r')
legend('poblacion de riesgo','enfermos','location','best');
hold off
clear('i');
for i=1:N4
seuler4(i+1)=seuler4(i)+h4*(-0.003*seuler4(i)*ieuler4(i));
ieuler4(i+1)=ieuler4(i)+h4*((0.003*seuler4(i)*ieuler4(i))-0.31*ieuler4(i));
end
figure(4);
hold on
plot(t4,seuler4)
plot(t4,ieuler4,'r')
legend('poblacion de riesgo','enfermos','location','best');
hold off
figure(5)
hold on
plot(t1,seuler1,'k');
plot(t2,seuler2,'g');
plot(t3,seuler3,'r');
plot(t4,seuler4);
legend('0.1','0.01','0.001','0.0001','location','best');
hold off 
figure(6)
hold on
plot(t1,ieuler1,'k');
plot(t2,ieuler2,'g');
plot(t3,ieuler3,'r');
plot(t4,ieuler4);
legend('0.1','0.01','0.001','0.0001','location','best');
hold off 
%calculamos el maximo y su posicion, usamos el paso mas pequeño dado que nos dara mayor precision 
[maximo posicion]=max(ieuler4(1,:));
tiempohastaelmaximo=(posicion-1)*h4
maximo
}}<br/>

===. <math> (S_{0},I_{0})=(800,20) </math>===
Ésta son dos de las gráficas resultado:
<br/>
{|
|-
| [[Archivo:8002001.jpg|thumb|300px|left|Gráfica para h=0.1]] || [[Archivo:80020001.jpg|thumb|300px|left|Gráfica para h=0.01]]
|}<br/>

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

<br/>

===. <math> (S_{0},I_{0})=(10000,40) </math>===
Se obtendrán las siguientes gráficas:<br/>
{|
|-
| [[Archivo:1401.jpg|thumb|300px|left|Gráfica para h=0.1]] || [[Archivo:14001.jpg|thumb|300px|left|Gráfica para h=0.01]]
|}<br/>

Se puede observar, que para el paso <math>h=0.1</math>, el resultado obtenido por la gráfica es incongruente, debido a que el método de Euler tiene deficiencias cuando se trabaja con derivadas elevadas como es el caso. Sin embargo, al disminuirlo a <math>h=0.01</math> o valores inferiores, la gráfica guarda más relación con la solución real.

Ésto ocurre porque el número inicial de personas infectadas <math>I_{0}</math> es muy elevado y al disminuir el valor del paso, el programa aproxima mejor la solución.<br/>


Ésta es la comparación de las distintas poblaciones para los diferentes valores de h:

[[Archivo:diego3.jpg|450px|thumb||left|Comparación entre las gráficas obtenidas para la población de enfermos para cada paso con (S0,I0)=(10000,40)]]
[[Archivo:diego4.jpg|450px|thumb||right|Comparación entre las gráficas obtenidas para la población de susceptible para cada paso con (S0,I0)=(10000,40)]]



































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}=9465</math> y cuándo se producirá <math>t=0.34</math> días. Al haber utilizado <math>h=10^{-4}</math>, el error mencionado anteriormente no influye.

==. 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:<br/>

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


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: 
[[Archivo:RK80020.jpg|350px|thumb||right|Runge-Kutta para (S0,I0)=(800,20)]]
[[Archivo:RK1000040.jpg|350px|thumb||right|Runge-Kutta para (S0,I0)=(1000,40)]]
{{matlab| codigo =
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=input('Introduzca el valor S0: ');
I0=input('Introduzca el valor I0: ');
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')
}}

Este método comparado con Euler, da lugar a las gráficas:<br/>
[[Archivo:RKcomp80020.jpg|400px|thumb||left|Runge-Kutta para (S0,I0)=(800,20)]]
[[Archivo:RKcomp1000040.jpg|400px|thumb||centre|Runge-Kutta para (S0,I0)=(10000,40)]]<br/>

En la primera gráfica, se observa que los métodos de Euler y Runge-Kutta dan resultados similares.<br/> 
Sin embargo, en la segunda gráfica, como ya se ha explicado anteriormente, el método de Euler para los valores iniciales <math> (S_{0},I_{0})=(10000,40)</math> conforme aumenta el paso, más impreciso se vuelve. En cambio, el método Runge-Kutta consigue aproximarlo de manera más exacta.<br/>


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.

==. 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:<br/> <math> a(t)=\frac{0.003}{1+t} </math> <br/>


Además se tiene unos datos iniciales: <br/>
<math>\left\{\begin{matrix}
I_{0}=40\\ 
S+I=1640
\end{matrix}\right.</math><br/>


Tomando dichos datos, se ha utilizado el método de Heun para aproximar la solución y así poder dibujar la gráfica:<br/>
[[Archivo:avariable.jpg|400px|thumb||right|]]
{{matlab| codigo=
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.<br/>


Se puede observar que a medida que pasa el tiempo el coeficiente variable <math>a</math> 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 el número de infectados empiece a descender rápidamente, llegando incluso hasta su desaparición, ya sea por fallecimiento o curación de los mismos. En consecuencia, el número de personas susceptibles a contraer la enfermedad también descenderá, llegando incluso a permanecer constante cuando la población de infectados desaparezca.

==. 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. <br/>


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>. <br/>
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, momento en el que el número de personas infectadas alcanza su punto más alto.<br/>


Se han utilizado los siguientes datos iniciales:<br/>
<math>\left\{\begin{matrix}
I_{0}=40\\ 
S+I=1640\\
a(t)=\frac{0.003}{1+t}
\end{matrix}\right.</math><br/>


{{matlab| codigo=
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)
}}

El programa devuelve el valor para el cual <math>a</math> alcanza su máximo para un tiempo <math>t=5</math> días. Éste valor será <math>a=0.0008</math>.