MateWiki - Contribuciones del usuario [es]

Impacto acústico de la carretera LO-20

2016-12-09T18:54:33Z

Trabajocampos: /* Conclusiones */

{{ TrabajoSIG | IMPACTO ACÚSTICO DE LA CARRETERA L0-20 | Pablo Revuelta Aragón, Pedro Torrecilla Sánchez, Gonzalo Bolea | [[:Categoría:SIGAIC_16/17|Curso 16/17]] }}

El objeto de trabajo es un estudio del impacto acústico provocado por el tráfico rodado de la carretera LO-20 que pasa por el medio de Logroño. Para ello hemos recopilado diversos datos acerca de la carretera, los parámetros del tráfico y tipo de vehículos. Basándonos en la tecnología que nos posibilitan los programas de SIG, en especial QGIS, que es un programa de código libre, y con la ayuda del complemento OpeNoise hemos podido realizar una serie de mapas comparativos indicándonos la afección que sufren los vecinos de Logroño en función de la distancia y variando distintos parámetros de la carretera como son el tipo de superficie, la velocidad y probando con limitaciones como la de prohibir el tráfico de vehículos pesados. También hemos calculado el impacto de dicho ruido sobre una futura guardería con y sin barreras acústicas.

Este trabajo tiene el claro objetivo de verificar que Logroño se trata de una ciudad segura en el ámbito acústico para todos los vecinos. Por ello hemos hecho un breve resumen de como afecta dicho impacto a la salud de las personas.

== Introducción ==
Para empezar a plantear el problema lo primero es situar la carretera, que es la siguiente:

[[Archivo:Foto1R.jpg|500px|miniaturadeimagen|thumb|centro|figure 1]]

Como se puede observar, la carretera, pasa muy cerca del centro de Logroño por lo que tiene mucho interés el estudio del impacto acústico ya que es un problema, sufrido por una gran parte de la población en la actualidad. Además observando la distribución de la ciudad nos hace pensar que el desarrollo urbanístico de Logroño se expandirá hacia el lado Sur de esta carretera por lo que será indispensable realizar un estudio del impacto acústico.

En España parte de la población vive afectada por este problema. A continuación podemos ver en la gráfica el número de personas que sufren a consecuencia de dicha molestia.

[[Archivo:Foto2.png|500px|miniaturadeimagen|thumb|centro|figure 2]]

Con esta imagen queda de manifiesto la existencia de un grave problema con el ruido producido por el transporte, ya que como se puede ver las carreteras son el principal causante de contaminación acústica en cuanto al transporte de nuestro país, ya que 2.000.000 de personas sufren las consecuencias de lo que posiblemente pueda considerarse una mala gestión y planificación de este problema. A partir de 55 dB las molestias ocasionadas por el ruido ya se pueden considerar molestias moderadas.
En nuestra carretera estudiaremos todas estas zonas en las que se supera dicha cantidad en distintas situaciones que serán:

- Según el tipo de superficie de la carretera.

- Con la restricción del tráfico de vehículos pesados.

- Aumentando la velocidad de 80 km/h a 100 km/h.

También diseñaremos una guardería en las inmediaciones de la carretera en la que estudiaremos dicho impacto acústico y realizaremos una barrera acústica para paliar estos problemas que pueden causar sobre el desarrollo intelectual y psicológico de los niños.

== Metodología ==

Para realizar el estudio principalmente hemos utilizado la herramienta OpeNoise de QGIS que se basa en el método nacional de cálculo francés NMPB-Routes-96, como se indica en la Directiva 2002/49/CE DEL PARLAMENTO EUROPEO, sobre evaluación y gestión del ruido ambiental. El modelo considera simplificaciones como el modelamiento es en 2D (sólo se considera que los puntos de recepción están a 4m del suelo), no se calculan reflexiones o difracciones y el terreno es plano.

Ha habido que buscar una serie de datos como la IMD de vehículos ligeros, la IMD de vehículos pesados, la velocidad de dicha carretera para vehículos ligeros y pesados, el tipo de tráfico y la superficie. Todos estos datos han sido extraídos de la página del Ministerio de Fomento.

- IMD de vehículos ligeros: 24.3990
- IMD de vehículos pesados: 3.328
- Velocidad de vehículos ligeros: 80
- Velocidad de vehículos pesados: 70
- Tipo de superficie de la carretera: lisa
- Tipo de tráfico: continuo

A continuación podemos ver la ventana en la que se incluyen los datos dentro del programa.

[[Archivo:Foto_3.png|500px|miniaturadeimagen|thumb|centro|figure 3]]

[[Archivo:Foto4.png|500px|miniaturadeimagen|thumb|centro|figure 4]]

En primer lugar tuvimos que buscar tanto el mapa ráster como el mapa vectorial de la zona para poder trabajar con ellos dentro del programa. Además conseguimos un mapa vectorial de todos los edificios de la zona, imprescindibles para el cálculo del impacto del sonido ya que harán de barreras del sonido.

Dentro del programa introducimos la carreta que deseamos estudiar como una línea. Para comenzar a ejecutar el complemento es necesario realizar una capa de puntos aleatorios en la zona que vamos a trabajar, entonces OpeNoise calculará el grado de afección del sonido a cada uno de estos puntos. El programa cada cierta distancia hace un punto en la línea, y desde este emite un haz de sonido que va en todas las direcciones, que disminuye con la distancia y se anula al chocar contra un objeto, como pueden ser los edificios. El OpeNoise va haciendo este proceso en cada diferencial de longitud (o intervalo escogido), y le transmite a una capa previamente creada de puntos el valor de los decibelios como atributo. Este complemento trabaja con los edificios de forma homogénea, es decir, no hace diferenciación entre un edificio y otro en función de los materiales, ni de la altura valor que no tendrá en cuenta para calcular el impacto sonoro. Una vez introducidos todos los datos necesarios para la creación del mapa obtendremos los mapas finales mediante una interpolación de los puntos anteriormente creados.
== Resultados ==

Una vez ejecutado el programa obtendremos una capa de puntos en la que se representa la afección del ruido a cada uno de estos puntos. El cálculo del impacto y la obtención de este resultado suele tardar un rato ya que la capa que introducíamos al complemento OpeNoise contenía 8000 mil puntos. Una vez obtenidos estos resultados realizamos una interpolación de todos los puntos lo cual nos permitía llegar a nuestro resultado final, que son unos mapas en los que quedan reflejados el impacto del sonido sobre el centro de la ciudad de Logroño. Los mapas de los puntos no vamos a reflejarlos en el trabajo ya que no reflejan ninguna información relevante que no nos proporcionen los mapas una vez interpolados. Los diferentes mapas obtenidos cambian en función de variaciones que hemos introducido en los distintos factores que afectan al sonido: velocidad, tipo de vehículos, superficie de rodadura, etc.

[[Archivo:Condiciones_actuales.png|500px|miniaturadeimagen|thumb|centro|Condiciones actuales]]

[[Archivo:Condiciones_actuales_mas_50.png|500px|miniaturadeimagen|thumb|centro|Condiciones actuales con solo la zona que supera los 50dB]]

[[Archivo:Sin_pesados.png|500px|miniaturadeimagen|thumb|centro|Sin vehículos pesados]]

[[Archivo:Sin_pesados_mas_50.png|500px|miniaturadeimagen|thumb|centro|Sin vehículos pesados con solo la zona que supera los 50dB]]

[[Archivo:Superficie_porosa.png|500px|miniaturadeimagen|thumb|centro|Con superficie porosa]]

[[Archivo:Superficie_porosa_mas_50.png|500px|miniaturadeimagen|thumb|centro|Con superficie porosa con solo la zona que supera los 50dB]]

== Conclusiones ==
El ruido generado por el tráfico es un grave problema que afecta a la vida de las personas. Este impacto puede llegar a generar problemas de salud de las personas afectadas. En caso de que el ruido supere los 50dB comenzará a poder provocar daños en la salud de las personas.

En el primer mapa hemos representado la zona como se encuentra en la actualidad, con una velocidad de 80km/h tanto para vehículos ligeros como para vehículos pesados. Como podemos observar el ruido supera los 50dB pero hemos de recordar que no hemos tenido en cuenta que los edificios tienen una altura y por lo tanto el ruido se absorbe en mayor medida. Para ver mejor las afecciones del impacto hemos realizado este mismo mapa en el que se representa solo las zonas en las que se supera dicho umbral. A continuación hemos realizado este mismo caso si se aumentara la velocidad máxima en 10km/h, y observamos supone un aumento de 2dB el máximo nivel sonoro.

Otro de los casos que hemos estudiado es la limitación de circulación únicamente a vehículos ligeros, es decir, sin vehículos pesados y el resultado obtenido revela que el ruido es generado por un tráfico continuo y por las condiciones de la carretera ya que este aumento de peso en los vehículos no genera un aumento sustancial de ruido.

Otra característica de las que hemos modificado es la superficie de rodadura, en el resto de casos hemos realizado el estudio con una superficie lisa, y en este caso hemos puesto una superficie rugosa que en teoría debería generar menos ruido que con una superficie lisa, ya que el ruido es absorbido por la propia carretera. Con los resultados obtenidos podemos observar que efectivamente el ruido disminuye en 2dB respecto a las condiciones iniciales que serán las mismas que en este caso pero con otro tipo de superficie.

En cuanto al colegio, dada su proximidad a la carretera, hay un ruido intolerable, por lo que hay que hacer medidas contra este. Esto debe a que es una carretera con una IMD bastante elevada, lo que imposibilitaría el aprendizaje a los estudiantes, además que atentaría gravemente contra su salud y capacidades cognitivas. Este impacto lo podemos observar en el mapa del ruido del colegio sin barrera, y vemos que en todo el área del centro hay niveles de ruido muy superiores a los recomendados estando incluso restringidos por la ley cerca de estos lugares.

Debido a lo comentado anteriormente, hemos diseñado una barrera contra el sonido. La primera abarcaría la extensión del centro, unos 400m, pero podemos observar que no es suficiente, ya que en las zonas que no están entre los edificios hay niveles muy altos de ruido, alcanzando los 67dB, lo que imposibilita el uso de este tipo de barrera, puesto que no cumple con sus objetivos.

La segunda barrera abarca una extensión mayor, de unos 700m. Esta barrera se extiende sobre todo hacia el Oeste, ya que es un descampado por lo que el ruido se propagaría mucho mas fácilmente. En el lado opuesto hay edificios, y carreteras, lo que nos imposibilita prolongarla más, aunque los edificios de este lado también favorecen nuestro propósito, puesto que actúa como barrera. Con esta disposición tendríamos que disponer el lugar de descanso de los estudiantes hacia el Oeste, con el objetivo de que no les importune el ruido causado por la LO-20.

[[Categoría:Sistemas de Información Geográfica Aplicados a la Ingeniería Civil]]
[[Categoría:SIGAIC_16/17]]

Impacto acústico de la carretera LO-20

2016-12-09T18:52:49Z

Trabajocampos: /* Resultados */

Impacto acústico de la carretera LO-20

2016-12-09T18:51:22Z

Trabajocampos: /* Resultados */

Archivo:Superficie porosa mas 50.png

2016-12-09T18:51:09Z

Trabajocampos:

Archivo:Superficie porosa.png

2016-12-09T18:49:26Z

Trabajocampos:

Archivo:Sin pesados mas 50.png

2016-12-09T18:48:08Z

Trabajocampos:

Impacto acústico de la carretera LO-20

2016-12-09T18:45:46Z

Trabajocampos:

Archivo:Sin pesados.png

2016-12-09T18:45:01Z

Trabajocampos:

Archivo:Condiciones actuales mas 50.png

2016-12-09T18:42:51Z

Trabajocampos:

Archivo:Condiciones actuales.png

2016-12-09T18:40:16Z

Trabajocampos:

Impacto acústico de la carretera LO-20

2016-12-09T15:09:16Z

Trabajocampos:

Impacto acústico de la carretera LO-20

2016-12-09T15:08:46Z

Trabajocampos:

Impacto acústico de la carretera LO-20

2016-12-09T15:04:39Z

Trabajocampos: /* Metodología */

Archivo:Foto4.png

2016-12-09T15:04:00Z

Trabajocampos:

Archivo:Foto 3.png

2016-12-09T15:02:56Z

Trabajocampos:

Impacto acústico de la carretera LO-20

2016-12-09T14:58:15Z

Trabajocampos: /* Introducción */

Impacto acústico de la carretera LO-20

2016-12-09T14:57:22Z

Trabajocampos: /* Introducción */

Archivo:Foto2.png

2016-12-09T14:56:27Z

Trabajocampos:

Impacto acústico de la carretera LO-20

2016-12-09T14:53:25Z

Trabajocampos: Página creada con «{{ TrabajoSIG | IMPACTO ACÚSTICO DE LA CARRETERA L0-20 | Pablo Revuelta Aragón, Pedro Torrecilla Sánchez, Gonzalo Bolea | Curso 16/17 }} E...»

Archivo:Foto1R.jpg

2016-12-09T14:50:53Z

Trabajocampos:

Modelos epidemiológicos (Grupo 3A)

2016-05-12T06:48:03Z

Trabajocampos:

{{ TrabajoED | Modelos Epidemológicos (Grupo 3A) | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED15/16|Curso 2015-16]] |
Ignacio Mollá Carcaño 
Pablo Revuelta Aragón 
David González Hernández 
Jose María García Rodríguez 
Alejandro Martínez Gamonal }}

En el desarrollo de una epidemia, se distinguen dos tipos de individuos. Los que ya han contraído la enfermedad (infectados), que llamaremos <math>I(t)</math> ; y los que son susceptibles de contraerla (sanos), a los que llamaremos <math>S(t)</math>. Donde <math>t</math> es la variable del tiempo. Se presentan dos hipótesis para realizar este estudio:

1. La población de personas infectadas se altera por el fallecimiento o la cura de las mismas. En ambos casos, la tasa de cambio depende del número 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ón entre el número de individuos en ambas clases.

A través del siguiente sistema de ecuaciones diferenciales en función del tiempo observamos las variaciones de ambas poblaciones:

<math>\begin{cases}\frac{dS}{dt}=-aSI\\\frac{dI}{dt}=aSI-bI-cI\end{cases}</math>

Donde <math>a</math>, <math>b</math> y <math>c</math> son parámetros.

* <math>\frac{dS}{dt}</math> : Es la variación de la población susceptible de contraer la enfermedad.

* <math>\frac{dI}{dt}</math> : Es la variación de la población infectada.

* <math>aSI</math> : Es la interacción entre ambas poblaciones que dependen de <math>a</math>.

* <math>bI</math> : Indica los individuos que fallecen. Por lo tanto <math>b</math> podría considerarse como la '''tasa de curación'''.

* <math>cI</math> : Indica los individuos que se han curado. Por lo tanto <math>c</math> podría considerarse como la '''tasa de mortalidad'''.

Para interpretar <math>a</math> observamos el sistema y la segunda hipótesis. Como los individuos que pasan de ser susceptibles a estar infectados son proporcionales a su interacción <math>(SI)</math>, <math>a</math> podría considerarse la '''tasa de contagio''' en la población susceptible.

==Resolución del problema con una sola ecuación diferencial==
Primero vamamos a estudiar el modelo epidemiológico de una población en la que incialmente las 2000 personas que conforman la población se encuentran infectadas. Para hacer el seguimiento de la enfermedad contamos tanto con la constante de proporcionalidad de personas que se curan (b=0.3) y la constante de proporcionalidad del número de personas que fallecen (c=0.01). De esta forma nos quedará el siguiente PVI:

<math> \begin{cases} \frac{dI}{dt} = -0.31 I \\I_{0} = 2000\end{cases} </math>

A continuación vamos a proceder a resolverlo por el método de Euler y trapecio asignandole un tamaño de paso de h=0.1 y veremos cuanto tiempo tardará en reducirse el número de infectados a la cuarta parte.

<math>\Rightarrow</math> Método de Euler
[[Archivo:Ejerciciodoseuler.png|500px|thumb|rigth|Aproximación mediante Euler]]
{{matlab|codigo=

clear all
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
y0=2000;
b=0.3;
c=0.01;

t(1)=t0;
ye(1)=y0; %Vector solución Euler

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de y(i) hasta que y(i)>500

%EULER

while ye(i)>500

ye(i+1)=ye(i)+h*(-b*ye(i)-c*ye(i));
t(i+1)=t(i)+h;
i=i+1;
end

[t',ye'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,ye)
legend('Población infectada(Método de Euler)','Location','best');

}}

<math>\Rightarrow</math> Método del trapecio
[[Archivo:Ejerciciodostrapecio.png|500px|thumb|rigth|Aproximación mediante el método del trapecio]]
{{matlab|codigo=
clear all
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
z0=2000;
b=0.3;
c=0.01;

t(1)=t0;
z(1)=z0; %Vector solución trapecio

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de z(i) hasta que z(i)>500

%TRAPECIO
while z(i)>500
z(i+1)=(z(i)*(1-(h/2)*(b+c)))/(1+h/2*(b+c));%trapecio
t(i+1)=t(i)+h;
i=i+1;
end

[t',z'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,z)
legend('Población infectada(Método del trapecio)','Location','best');

}}

Con todo esto podemos determinar que el número de infectados se reducirá a la cuarta parte en 4.5 días, esto ocurrirá o bien por el fallecimiento o por la curación de dicha enfermedad. Cabe señalar que el número de curados será mayor que el de defunciones, esto se puede ver a simple vista viendo que b>c.
De cara a nuestro estudio de esta enfermedad vamos a ver como afectaría al desarrollo de la misma si introducimos 100 sujetos sanos en la población, tomando una constante de interacción entre personas infectadas y sanas de a=0.003. Procederemos igual que antes a resolverlo por el método de Euler y trapecio con un tamaño de paso de 0.1.
Primero vamos a definir nuestro nuevo PVI:
<math> \begin{cases} \frac{dI}{dt} = -0.01 I \\I_{0} = 2000\end{cases} </math>
<math>\Rightarrow</math> Método de Euler
[[Archivo:Ejerciciotreseuler.png|500px|thumb|rigth|Aproximación mediante Euler para S=100]]
{{matlab|codigo=

%Euler
clear all
%Datos
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%a: Es la constante de proporcionalidad entre los infectados por cada
%interacción
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
y0=2000;
S=100;
a=0.003;
b=0.3;
c=0.01;

t(1)=t0;
y(1)=y0;

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de y(i) hasta que y(i)>500
while y(i)>500

y(i+1)=y(i)+h*(a*S*y(i)-b*y(i)-c*y(i));
t(i+1)=t(i)+h;
i=i+1;
end

[t',y'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,y)
legend('Población infectada para S=100(Método de Euler)','Location','best');

}}

<math>\Rightarrow</math> Método del trapecio
[[Archivo:Ejerciciotrestrapecio.png|500px|thumb|rigth|Aproximación mediante el método del trapecio para S=100]]
{{matlab|codigo=

%Trapecio
clear all
%Datos
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%a: Es la constante de proporcionalidad entre los infectados por cada
%interacción
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
z0=2000;
S=100;
a=0.003;
b=0.3;
c=0.01;

t(1)=t0;
z(1)=z0;

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de z(i) hasta que z(i)>500

%TRAPECIO
while z(i)>500
z(i+1)=(z(i)*(1-(h/2)*(b+c-S*a)))/(1+h/2*(b+c-S*a));
t(i+1)=t(i)+h;
i=i+1;
end

[t',z'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,z)
legend('Población infectada para S=100(Método del trapecio)','Location','best');

}}

Ahora con estos nuevos datos vemos que para que se reduzca a la cuarta parte el número de infectados tendrán que pasar 138.7 días, algo lógico al fijarnos en la ED, ambas son ecuaciones decrecientes pero para el caso de S=100 la pendiente será mucho menor.
Si continuamos haciendo pruebas vemos para valores superiores de s, por ejemplo, S=200 el matlab no nos da ninguna respuesta. Esto ocurre porque al aumentar S la ecuación pasa a ser creciente, por lo que lógicamente si intentas buscar un valor del número de infectados menor al inicial no lo podrás encontrar.

==Resolución del problema completo==
Vamos a resolver el problema completo para ver como evolucionarían las poblaciones en un periodo de 40 días.
Primero supongamos una población infectada inicial de 20 individuos y una población en riesgo de contagio de 800.
Después supondremos el caso de 40 individuos enfermos y 10000 en riesgo de contagio.
Para resolverlo usaremos el método de Euler y más tarde el método de Runge-Kutta de orden cuatro con distintas discretizaciones para ver la influencia de estas.

<math>\Rightarrow</math> Método de Euler
{{matlab|codigo=

%Euler
clear all
%Datos
%Introducimos los valores de las constantes
a=0.003;
b=0.3;
c=0.01;
h=input('Introducir valores de h:');
t0=0;
y0=input('Introducir vector [So,Io]:');
tN=40;
% calculamos los subintervalos
t=t0:h:tN;
N=(tN-t0)/h;
y=zeros(2,length(t));
y(:,1)=y0';
for i= 1:N;
y(:,i+1)=y(:,i)+h*[-a*y(1,i).*y(2,i);a*y(1,i).*y(2,i)-b*y(2,i)-c*y(2,i)];
end

% grafico
plot(t,y)
% vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h
legend('Población sana','Location','best','Población enferma','Location','best');

}}

<math>\Rightarrow</math> Método de Runge-Kutta

{{matlab|codigo=

%Runge-Kutta
clear all

%Datos
S0=input('Introducir valor inicial SO: ');
I0=input('Introducir valor inicial IO: ');
a=0.003;
b=0.3;
c=0.01;
%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector solucion
y0=[S0;I0];
y=zeros(2,N+1);
y(:,1)=y0;
for i=1:N
K1S=-a*y(1,i)*y(2,i);
K2S=-a*(y(1,i)+(1/2)*K1S*h)*(y(2,i)+(1/2)*K1S*h);
K3S=-a*(y(1,i)+(1/2)*K2S*h)*(y(2,i)+(1/2)*K2S*h);
K4S=-a*(y(1,i)+K3S*h)*(y(2,i)+K3S*h);
y(1,i+1)=y(1,i)+(h/6)*(K1S+2*K2S+2*K3S+K4S);
K1I=a*y(1,i)*y(2,i)-b*y(2,i)-c*y(2,i);
K2I=a*(y(1,i)+(1/2)*K1I*h)*(y(2,i)+(1/2)*K1I*h)-b*(y(2,i)+(1/2)*K1I*h)-c*(y(2,i)+(1/2)*K1I*h);
K3I=a*(y(1,i)+(1/2)*K2I*h)*(y(2,i)+(1/2)*K2I*h)-b*(y(2,i)+(1/2)*K2I*h)-c*(y(2,i)+(1/2)*K2I*h);
K4I=a*(y(1,i)+K3I*h)*(y(2,i)+K3I*h)-b*(y(2,i)+K3I*h)-c*(y(2,i)+K3I*h);
y(2,i+1)=y(2,i)+(h/6)*(K1I+2*K2I+2*K3I+K4I);
end

%Tabla de resutados
[t',y(1,:)',y(2,:)']
%Gráfico
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'k')
legend('Poblacion sana','Poblacion enferma','Location','best');
hold off

% vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h

}}

La dificultad que hay en el uso del método del trapecio o cualquier otro método implícito está en que las incógnitas de nuestro sistema <math>(S,I)</math> quedan implícitas en ambas ecuaciones y necesitan que se despejen. Al comenzar a realizar cálculos nos damos cuenta que el despeje es imposible, ya que las dos variables incógnita aparecen multiplicándose entre sí, quedando una variable dependiendo de la otra siempre.

===Situación inicial de <math>(S _{0} , I_{0} )=(800,20)</math>===
En estos programas, sustituyendo las <math>h</math> para 0.1, 0.01, 0.001 y 0.0001 obtendremos los siguientes resultados:

-Para <math>h=0.1</math>

[[Archivo:Tborraré.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.1rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

Como se puede ver la población sana se infecta rápidamente y a partir del día 15 prácticamente toda la población ha fallecido.
* '''Método de Euler''': nuestro programa nos indica que en este caso el momento en el que los enfermos son máximos se produce al '''tercer día (2.8)''' y el '''número de enfermos en ese momento es 517'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce al '''tercer día (2.6)''' y el '''número de enfermos en ese momento es 650'''.

-Para <math>h=0.01</math>
[[Archivo:Tborraré2.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.01rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

En el gráfico apenas se encuentra diferencia, lo cual es lógico porque es el mismo problema, pero si sobrepusiéramos las gráficas, lograríamos ver como se va ajustando a la realidad en función del aumento de h.
* '''Método de Euler''': en este caso el momento en el que los enfermos son máximos es también en el '''tercer día (2.7)''' y el número de '''enfermos es 506'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce al '''tercer día (2.69)''' y el número de '''enfermos en ese momento es 517'''.

-Para <math>h=0.001</math>
[[Archivo:Tborraré3.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.001rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': en este caso el momento en el que los enfermos son máximos se produce el '''tercer día (2.696)''' con un valor de '''505 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''tercer día (2.694)''' con un valor de '''506 enfermos'''.

-Para <math>h=0.0001</math>
[[Archivo:Tborraré4.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.0001rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': al igual que antes, se repite el momento en el que los enfermos son máximos en el '''tercer día (2.695)''' con un valor de '''505 enfermos'''.
* '''Método de Runge-Kutta''': se repite el momento en el que los enfermos son máximos en el '''tercer día (2.6948)''' con un valor de '''505 enfermos'''.

Observamos en las distintas gráficas obtenidas, que según vamos empleando un tamaño de paso <math>h</math> más pequeño, nos vamos acercando cada vez más a la población real que pretendemos reflejar.
Nos damos cuenta de que las gráficas tienen bastante lógica, pues a la vez que desciende el número de población sana va aumentando la población infectada.
También observamos que para el método de Runge-Kutta cuando utilizamos tamaño de paso <math>h=0.1</math> pierde precisión en relación a cuando utilizamos tamaños de paso más pequeños.

===Situación inicial de <math>(S _{0} , I_{0} )=(10000,40)</math>===
En estos programas, sustituyendo las <math>h</math> para 0.1, 0.01, 0.001 y 0.0001 obtendremos los siguientes resultados:

-Para <math>h=0.1</math>

[[Archivo:H0.11000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.5)''' con un valor de '''12600 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.4)''' con un valor de '''2.6e74 enfermos'''.

-Para <math>h=0.01</math>
[[Archivo:H0.0110000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.0110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.35)''' con un valor de '''9521 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.31)''' con un valor de '''13550 enfermos'''.

-Para <math>h=0.001</math>
[[Archivo:H0.00110000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.00110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.342)''' con un valor de '''9470 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.338)''' con un valor de '''9758 enfermos'''.

-Para <math>h=0.0001</math>
[[Archivo:H0.000110000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.000110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.3411)''' con un valor de '''9465 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.3407)''' con un valor de '''9493 enfermos'''.

Para este segundo caso, nos damos cuenta que al utilizar tamaño de paso <math>h=0.1</math> , llegamos a un resultado sin sentido, pues nos da valores negativos y es imposible que una población tenga tales valores.
En las demás gráficas, los resultados obtenidos no parecen cuadrar con la realidad. Esto podría ser debido a que la alta diferencia entre los valores iniciales <math>(S _{0} , I_{0} )=(10000,40)</math> no permiten un buen cálculo de su variación respecto del tiempo.
El número inicial de infectados es muy elevado y al disminuir el tamaño de paso, el programa aproxima mejor la solución.

==Resolución del problema completo suponiendo <math> a </math> una función variable con el tiempo==
Para el siguiente caso supondremos que el coeficiente <math> a </math> (tasa de contagio en población sana) es una función dada por :

<math> a(t)= \frac{0.003}{(1+t)} </math>
Suponiendo además las condiciones iniciales <math> (S_{0},I_{0})=(1640,40) </math> , vamos a dibujar las gráficas de infectados y sanos por el método de Heun, tomando un tamaño de paso <math> h=0.01 </math>
[[Archivo: H0.01heunn.png|1000px|thumb|rigth| Población infectada-Población sana para a=a(t)]]
{{matlab|codigo=
%Heun. Coeficiente variable a(t)

%Datos
S0=input('Introducir valor inicial SO: ');
I0=input('Introducir valor inicial IO: ');
b=0.3;
c=0.01;
%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector solucion
y0=[S0;I0];
y=zeros(2,N+1);
y(:,1)=y0;
for i=1:N
K1S=-(0.003/(1+t(i)))*y(1,i)*y(2,i);
K2S=-(0.003/(1+t(i)+h))*(y(1,i)+K1S*h)*(y(2,i)+K1S*h);
y(1,i+1)=y(1,i)+(h/2)*(K1S+K2S);
K1I=(0.003/(1+t(i)))*y(1,i)*y(2,i)-b*y(2,i)-c*y(2,i);
K2I=(0.003/(1+t(i)+h))*(y(1,i)+K1I*h)*(y(2,i)+K1I*h)-b*(y(2,i)+K1I*h)-c*(y(2,i)+K1I*h);
y(2,i+1)=y(2,i)+(h/2)*(K1I+K2I);
end

%Tabla de resutados
[t',y(1,:)',y(2,:)']
%Graficos
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'k')
title('Método de Heun');
legend('Poblacion sana','Población infectada','Location','best');
hold off

% Vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h

}}

Como se puede ver la población sana se infecta rápidamente y a partir del día 15 casi toda la población ha fallecido. Al segundo día <math>(t=2,15)</math> observamos que el número de infectados es de 994 individuos <math>(I=994,3)</math>, siendo este el máximo que se alcanzará.

A medida que va pasando el tiempo el coeficiente <math> a(t) </math> (tasa de contagio en población sana) va a ir decreciendo poco a poco hasta aproximarse a cero según podemos ver en el siguiente gráfico:
<math> \lim_{t\to \infty } a(t)=0</math>
[[Archivo: A(t).png|450px|thumb|right|a(t) para t=[0,40]]]
{{matlab|codigo=

%tasa de contagio en poblacion sana a(t)

%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector a
a=zeros(1,N+1);
for i=1:N+1
a(i)=0.003/(1+t(i));
end

%Grafica
plot(t,a)

}}

[[Archivo: Comparacionheun.png|border|1500px|center|frame|Poblaciones según la tasa de contagio a(t)]]

Ambas poblaciones parten del mismo inicio (para <math>t=0</math> las poblaciones son iguales). Podemos observar que la población sana (susceptible de contraer la enfermedad) en el caso en el que la tasa de contagio en población sana <math>(a)</math> es constante <math>(a=0.003)</math> se ve reducida antes que en el caso en el que depende del tiempo <math>(a=a(t))</math>.

<math>a(t)</math> disminuye a la vez que transcurre el tiempo, por lo que en general las pendientes tanto de la población sana como de la infectada serán menores para este caso que para el caso en el que <math>a=0.003</math>

==Calibración del parámetro <math>a</math> en base a una experiencia previa==
Ahora vamos a intentar calibrar el párametro <math>a</math> para que el número de infectados tenga el máximo en el quinto día, para ello escribimos el siguiente código:

{{matlab|codigo=
%Aplicamos el Método de Heunn

%Datos
t0=0;
tN=40;
h=0.1;
t=t0:h:tN;
a=0.0005:0.0001:0.002;
b=0.3;
c=0.01;
N1=round((tN-t0)/h);
N2=length(a);
Tmax=zeros(1,N2);

%Vectores iniciales de Población e Infectados
y=zeros(2,N1+1);

%Condiciones iniciales
y(1,1)=1600;
y(2,1)=40;

%Bucle que nos calcula todas las posibles soluciones
for j=1:N2
for i=1:N1
K1S=-a(j)*y(1,i)*y(2,i);
K2S=-a(j)*y(1,i)*y(2,i)+K1S*h;
y(1,i+1)=y(1,i)+h/2*(K1S+K2S);

K1I=y(2,i)*(a(j)*y(1,i)-(b+c));
K2I=y(2,i)*(a(j)*y(1,i)-(b+c))+K1I*h;
y(2,i+1)=y(2,i)+h/2*(K1I+K2I);

if i>1 && y(2,i+1)<y(2,i) && y(2,i)>y(2,i-1)
Tmax(j)=t(i);
end
end
end
distancia=(abs(Tmax-5));
[valormin,posmin]=min(distancia);
a=a(posmin)
}}

Hemos encontrado la solución utilizando un bucle de tipo <math>for</math> dos veces, esto nos permite obtener todas los valores de <math>I</math>, infectados, en el tiempo transcurrido entre [0,40] días para cualquier <math>a</math> en el intervalo de [0.0005,0.002]. Una vez tenemos todas las posibilidades aplicamos una condición del tipo <math>if</math> para localizar donde se sitúa el máximo de infectados en los distintos <math>a</math> analizados, después vemos cuál de los máximos esta más próximo al quinto día, para ello generamos un vector <math>distancia</math> donde hacemos el valor absoluto de la resta entre el valor obtenido y el quinto día tras esto conseguimos el valor mas cercano. Ahora analizamos el valor de <math>a</math> para dicha posición obteniendo la solución pedida. En nuestro caso el valor pedida daba '''a = 8.0000e-04''' que era el cuarto <math>a</math> analizado en el bucle.

[[Categoría:Ecuaciones Diferentiales]]
[[Categoría:ED15/16]]
[[Categoría:Trabajos 2015-16]]

Modelos epidemiológicos (Grupo 3A)

2016-05-02T19:26:47Z

Trabajocampos: /* Calibración del parámetro a en base a una experiencia previa */

{{ TrabajoED | Modelos Epidemológicos (Grupo 3A) | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED15/16|Curso 2015-16]] |
Ignacio Mollá Carcaño 
Pablo Revuelta Aragón 
David González Hernández 
Jose María García Rodríguez 
Alejandro Martínez Gamonal }}

En el desarrollo de una epidemia, se distinguen dos tipos de individuos. Los que ya han contraído la enfermedad (infectados), que llamaremos <math>I(t)</math> ; y los que son susceptibles de contraerla (sanos), a los que llamaremos <math>S(t)</math>. Donde <math>t</math> es la variable del tiempo. Se presentan dos hipótesis para realizar este estudio:

1. La población de personas infectadas se altera por el fallecimiento o la cura de las mismas. En ambos casos, la tasa de cambio depende del número 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ón entre el número de individuos en ambas clases.

A través del siguiente sistema de ecuaciones diferenciales en función del tiempo observamos las variaciones de ambas poblaciones:

<math>\begin{cases}\frac{dS}{dt}=-aSI\\\frac{dI}{dt}=aSI-bI-cI\end{cases}</math>

Donde <math>a</math>, <math>b</math> y <math>c</math> son parámetros.

* <math>\frac{dS}{dt}</math> : Es la variación de la población susceptible de contraer la enfermedad.

* <math>\frac{dI}{dt}</math> : Es la variación de la población infectada.

* <math>aSI</math> : Es la interacción entre ambas poblaciones que dependen de <math>a</math>.

* <math>bI</math> : Indica los individuos que fallecen. Por lo tanto <math>b</math> podría considerarse como la '''tasa de mortalidad'''.

* <math>cI</math> : Indica los individuos que se han curado. Por lo tanto <math>c</math> podría considerarse como la '''tasa de curación'''.

Para interpretar <math>a</math> observamos el sistema y la segunda hipótesis. Como los individuos que pasan de ser susceptibles a estar infectados son proporcionales a su interacción <math>(SI)</math>, <math>a</math> podría considerarse la '''tasa de contagio''' en la población susceptible.

==Resolución del problema con una sola ecuación diferencial==
Primero vamamos a estudiar el modelo epidemiológico de una población en la que incialmente las 2000 personas que conforman la población se encuentran infectadas. Para hacer el seguimiento de la enfermedad contamos tanto con la constante de proporcionalidad de personas que se curan (b=0.3) y la constante de proporcionalidad del número de personas que fallecen (c=0.01). De esta forma nos quedará el siguiente PVI:

<math> \begin{cases} \frac{dI}{dt} = -0.31 I \\I_{0} = 2000\end{cases} </math>

A continuación vamos a proceder a resolverlo por el método de Euler y trapecio asignandole un tamaño de paso de h=0.1 y veremos cuanto tiempo tardará en reducirse el número de infectados a la cuarta parte.

<math>\Rightarrow</math> Método de Euler
[[Archivo:Ejerciciodoseuler.png|500px|thumb|rigth|Aproximación mediante Euler]]
{{matlab|codigo=

clear all
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
y0=2000;
b=0.3;
c=0.01;

t(1)=t0;
ye(1)=y0; %Vector solución Euler

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de y(i) hasta que y(i)>500

%EULER

while ye(i)>500

ye(i+1)=ye(i)+h*(-b*ye(i)-c*ye(i));
t(i+1)=t(i)+h;
i=i+1;
end

[t',ye'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,ye)
legend('Población infectada(Método de Euler)','Location','best');

}}

<math>\Rightarrow</math> Método del trapecio
[[Archivo:Ejerciciodostrapecio.png|500px|thumb|rigth|Aproximación mediante el método del trapecio]]
{{matlab|codigo=
clear all
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
z0=2000;
b=0.3;
c=0.01;

t(1)=t0;
z(1)=z0; %Vector solución trapecio

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de z(i) hasta que z(i)>500

%TRAPECIO
while z(i)>500
z(i+1)=(z(i)*(1-(h/2)*(b+c)))/(1+h/2*(b+c));%trapecio
t(i+1)=t(i)+h;
i=i+1;
end

[t',z'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,z)
legend('Población infectada(Método del trapecio)','Location','best');

}}

Con todo esto podemos determinar que el número de infectados se reducirá a la cuarta parte en 4.5 días, esto ocurrirá o bien por el fallecimiento o por la curación de dicha enfermedad. Cabe señalar que el número de curados será mayor que el de defunciones, esto se puede ver a simple vista viendo que b>a.
De cara a nuestro estudio de esta enfermedad vamos a ver como afectaría al desarrollo de la misma si introducimos 100 sujetos sanos en la población, tomando una constante de interacción entre personas infectadas y sanas de a=0.003. Procederemos igual que antes a resolverlo por el método de Euler y trapecio con un tamaño de paso de 0.1.
Primero vamos a definir nuestro nuevo PVI:
<math> \begin{cases} \frac{dI}{dt} = -0.01 I \\I_{0} = 2000\end{cases} </math>
<math>\Rightarrow</math> Método de Euler
[[Archivo:Ejerciciotreseuler.png|500px|thumb|rigth|Aproximación mediante Euler para S=100]]
{{matlab|codigo=

%Euler
clear all
%Datos
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%a: Es la constante de proporcionalidad entre los infectados por cada
%interacción
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
y0=2000;
S=100;
a=0.003;
b=0.3;
c=0.01;

t(1)=t0;
y(1)=y0;

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de y(i) hasta que y(i)>500
while y(i)>500

y(i+1)=y(i)+h*(a*S*y(i)-b*y(i)-c*y(i));
t(i+1)=t(i)+h;
i=i+1;
end

[t',y'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,y)
legend('Población infectada para S=100(Método de Euler)','Location','best');

}}

<math>\Rightarrow</math> Método del trapecio
[[Archivo:Ejerciciotrestrapecio.png|500px|thumb|rigth|Aproximación mediante el método del trapecio para S=100]]
{{matlab|codigo=

%Trapecio
clear all
%Datos
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%a: Es la constante de proporcionalidad entre los infectados por cada
%interacción
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
z0=2000;
S=100;
a=0.003;
b=0.3;
c=0.01;

t(1)=t0;
z(1)=z0;

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de z(i) hasta que z(i)>500

%TRAPECIO
while z(i)>500
z(i+1)=(z(i)*(1-(h/2)*(b+c-S*a)))/(1+h/2*(b+c-S*a));
t(i+1)=t(i)+h;
i=i+1;
end

[t',z'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,z)
legend('Población infectada para S=100(Método del trapecio)','Location','best');

}}

Ahora con estos nuevos datos vemos que para que se reduzca a la cuarta parte el número de infectados tendrán que pasar 138.7 días, algo lógico al fijarnos en la ED, ambas son ecuaciones decrecientes pero para el caso de S=100 la pendiente será mucho menor.
Si continuamos haciendo pruebas vemos para valores superiores de s, por ejemplo, S=200 el matlab no nos da ninguna respuesta. Esto ocurre porque al aumentar S la ecuación pasa a ser creciente, por lo que lógicamente si intentas buscar un valor del número de infectados menor al inicial no lo podrás encontrar.

==Resolución del problema completo==
Vamos a resolver el problema completo para ver como evolucionarían las poblaciones en un periodo de 40 días.
Primero supongamos una población infectada inicial de 20 individuos y una población en riesgo de contagio de 800.
Después supondremos el caso de 40 individuos enfermos y 10000 en riesgo de contagio.
Para resolverlo usaremos el método de Euler y más tarde el método de Runge-Kutta de orden cuatro con distintas discretizaciones para ver la influencia de estas.

<math>\Rightarrow</math> Método de Euler
{{matlab|codigo=

%Euler
clear all
%Datos
%Introducimos los valores de las constantes
a=0.003;
b=0.3;
c=0.01;
h=input('Introducir valores de h:');
t0=0;
y0=input('Introducir vector [So,Io]:');
tN=40;
% calculamos los subintervalos
t=t0:h:tN;
N=(tN-t0)/h;
y=zeros(2,length(t));
y(:,1)=y0';
for i= 1:N;
y(:,i+1)=y(:,i)+h*[-a*y(1,i).*y(2,i);a*y(1,i).*y(2,i)-b*y(2,i)-c*y(2,i)];
end

% grafico
plot(t,y)
% vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h
legend('Población sana','Location','best','Población enferma','Location','best');

}}

<math>\Rightarrow</math> Método de Runge-Kutta

{{matlab|codigo=

%Runge-Kutta
clear all

%Datos
S0=input('Introducir valor inicial SO: ');
I0=input('Introducir valor inicial IO: ');
a=0.003;
b=0.3;
c=0.01;
%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector solucion
y0=[S0;I0];
y=zeros(2,N+1);
y(:,1)=y0;
for i=1:N
K1S=-a*y(1,i)*y(2,i);
K2S=-a*(y(1,i)+(1/2)*K1S*h)*(y(2,i)+(1/2)*K1S*h);
K3S=-a*(y(1,i)+(1/2)*K2S*h)*(y(2,i)+(1/2)*K2S*h);
K4S=-a*(y(1,i)+K3S*h)*(y(2,i)+K3S*h);
y(1,i+1)=y(1,i)+(h/6)*(K1S+2*K2S+2*K3S+K4S);
K1I=a*y(1,i)*y(2,i)-b*y(2,i)-c*y(2,i);
K2I=a*(y(1,i)+(1/2)*K1I*h)*(y(2,i)+(1/2)*K1I*h)-b*(y(2,i)+(1/2)*K1I*h)-c*(y(2,i)+(1/2)*K1I*h);
K3I=a*(y(1,i)+(1/2)*K2I*h)*(y(2,i)+(1/2)*K2I*h)-b*(y(2,i)+(1/2)*K2I*h)-c*(y(2,i)+(1/2)*K2I*h);
K4I=a*(y(1,i)+K3I*h)*(y(2,i)+K3I*h)-b*(y(2,i)+K3I*h)-c*(y(2,i)+K3I*h);
y(2,i+1)=y(2,i)+(h/6)*(K1I+2*K2I+2*K3I+K4I);
end

%Tabla de resutados
[t',y(1,:)',y(2,:)']
%Gráfico
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'k')
legend('Poblacion sana','Poblacion enferma','Location','best');
hold off

% vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h

}}

La dificultad que hay en el uso del método del trapecio o cualquier otro método implícito está en que las incógnitas de nuestro sistema <math>(S,I)</math> quedan implícitas en ambas ecuaciones y necesitan que se despejen. Al comenzar a realizar cálculos nos damos cuenta que el despeje es imposible, ya que las dos variables incógnita aparecen multiplicándose entre sí, quedando una variable dependiendo de la otra siempre.

===Situación inicial de <math>(S _{0} , I_{0} )=(800,20)</math>===
En estos programas, sustituyendo las <math>h</math> para 0.1, 0.01, 0.001 y 0.0001 obtendremos los siguientes resultados:

-Para <math>h=0.1</math>

[[Archivo:Tborraré.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.1rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

Como se puede ver la población sana se infecta rápidamente y a partir del día 15 prácticamente toda la población ha fallecido.
* '''Método de Euler''': nuestro programa nos indica que en este caso el momento en el que los enfermos son máximos se produce al '''tercer día (2.8)''' y el '''número de enfermos en ese momento es 517'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce al '''tercer día (2.6)''' y el '''número de enfermos en ese momento es 650'''.

-Para <math>h=0.01</math>
[[Archivo:Tborraré2.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.01rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

En el gráfico apenas se encuentra diferencia, lo cual es lógico porque es el mismo problema, pero si sobrepusiéramos las gráficas, lograríamos ver como se va ajustando a la realidad en función del aumento de h.
* '''Método de Euler''': en este caso el momento en el que los enfermos son máximos es también en el '''tercer día (2.7)''' y el número de '''enfermos es 506'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce al '''tercer día (2.69)''' y el número de '''enfermos en ese momento es 517'''.

-Para <math>h=0.001</math>
[[Archivo:Tborraré3.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.001rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': en este caso el momento en el que los enfermos son máximos se produce el '''tercer día (2.696)''' con un valor de '''505 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''tercer día (2.694)''' con un valor de '''506 enfermos'''.

-Para <math>h=0.0001</math>
[[Archivo:Tborraré4.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.0001rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': al igual que antes, se repite el momento en el que los enfermos son máximos en el '''tercer día (2.695)''' con un valor de '''505 enfermos'''.
* '''Método de Runge-Kutta''': se repite el momento en el que los enfermos son máximos en el '''tercer día (2.6948)''' con un valor de '''505 enfermos'''.

Observamos en las distintas gráficas obtenidas, que según vamos empleando un tamaño de paso <math>h</math> más pequeño, nos vamos acercando cada vez más a la población real que pretendemos reflejar.
Nos damos cuenta de que las gráficas tienen bastante lógica, pues a la vez que desciende el número de población sana va aumentando la población infectada.
También observamos que para el método de Runge-Kutta cuando utilizamos tamaño de paso <math>h=0.1</math> pierde precisión en relación a cuando utilizamos tamaños de paso más pequeños.

===Situación inicial de <math>(S _{0} , I_{0} )=(10000,40)</math>===
En estos programas, sustituyendo las <math>h</math> para 0.1, 0.01, 0.001 y 0.0001 obtendremos los siguientes resultados:

-Para <math>h=0.1</math>

[[Archivo:H0.11000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.5)''' con un valor de '''12600 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.4)''' con un valor de '''2.6e74 enfermos'''.

-Para <math>h=0.01</math>
[[Archivo:H0.0110000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.0110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.35)''' con un valor de '''9521 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.31)''' con un valor de '''13550 enfermos'''.

-Para <math>h=0.001</math>
[[Archivo:H0.00110000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.00110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.342)''' con un valor de '''9470 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.338)''' con un valor de '''9758 enfermos'''.

-Para <math>h=0.0001</math>
[[Archivo:H0.000110000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.000110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.3411)''' con un valor de '''9465 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.3407)''' con un valor de '''9493 enfermos'''.

Para este segundo caso, nos damos cuenta que al utilizar tamaño de paso <math>h=0.1</math> , llegamos a un resultado sin sentido, pues nos da valores negativos y es imposible que una población tenga tales valores.
En las demás gráficas, los resultados obtenidos no parecen cuadrar con la realidad. Esto podría ser debido a que la alta diferencia entre los valores iniciales <math>(S _{0} , I_{0} )=(10000,40)</math> no permiten un buen cálculo de su variación respecto del tiempo.
El número inicial de infectados es muy elevado y al disminuir el tamaño de paso, el programa aproxima mejor la solución.

==Resolución del problema completo suponiendo <math> a </math> una función variable con el tiempo==
Para el siguiente caso supondremos que el coeficiente <math> a </math> (tasa de contagio en población sana) es una función dada por :

<math> a(t)= \frac{0.003}{(1+t)} </math>
Suponiendo además las condiciones iniciales <math> (S_{0},I_{0})=(1640,40) </math> , vamos a dibujar las gráficas de infectados y sanos por el método de Heun, tomando un tamaño de paso <math> h=0.01 </math>
[[Archivo: H0.01heunn.png|1000px|thumb|rigth| Población infectada-Población sana para a=a(t)]]
{{matlab|codigo=
%Heun. Coeficiente variable a(t)

%Datos
S0=input('Introducir valor inicial SO: ');
I0=input('Introducir valor inicial IO: ');
b=0.3;
c=0.01;
%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector solucion
y0=[S0;I0];
y=zeros(2,N+1);
y(:,1)=y0;
for i=1:N
K1S=-(0.003/(1+t(i)))*y(1,i)*y(2,i);
K2S=-(0.003/(1+t(i)+h))*(y(1,i)+K1S*h)*(y(2,i)+K1S*h);
y(1,i+1)=y(1,i)+(h/2)*(K1S+K2S);
K1I=(0.003/(1+t(i)))*y(1,i)*y(2,i)-b*y(2,i)-c*y(2,i);
K2I=(0.003/(1+t(i)+h))*(y(1,i)+K1I*h)*(y(2,i)+K1I*h)-b*(y(2,i)+K1I*h)-c*(y(2,i)+K1I*h);
y(2,i+1)=y(2,i)+(h/2)*(K1I+K2I);
end

%Tabla de resutados
[t',y(1,:)',y(2,:)']
%Graficos
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'k')
title('Método de Heun');
legend('Poblacion sana','Población infectada','Location','best');
hold off

% Vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h

}}

Como se puede ver la población sana se infecta rápidamente y a partir del día 15 casi toda la población ha fallecido. Al segundo día <math>(t=2,15)</math> observamos que el número de infectados es de 994 individuos <math>(I=994,3)</math>, siendo este el máximo que se alcanzará.

A medida que va pasando el tiempo el coeficiente <math> a(t) </math> (tasa de contagio en población sana) va a ir decreciendo poco a poco hasta aproximarse a cero según podemos ver en el siguiente gráfico:
<math> \lim_{t\to \infty } a(t)=0</math>
[[Archivo: A(t).png|450px|thumb|right|a(t) para t=[0,40]]]
{{matlab|codigo=

%tasa de contagio en poblacion sana a(t)

%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector a
a=zeros(1,N+1);
for i=1:N+1
a(i)=0.003/(1+t(i));
end

%Grafica
plot(t,a)

}}

[[Archivo: Comparacionheun.png|border|1500px|center|frame|Poblaciones según la tasa de contagio a(t)]]

Ambas poblaciones parten del mismo inicio (para <math>t=0</math> las poblaciones son iguales). Podemos observar que la población sana (susceptible de contraer la enfermedad) en el caso en el que la tasa de contagio en población sana <math>(a)</math> es constante <math>(a=0.003)</math> se ve reducida antes que en el caso en el que depende del tiempo <math>(a=a(t))</math>.

<math>a(t)</math> disminuye a la vez que transcurre el tiempo, por lo que en general las pendientes tanto de la población sana como de la infectada serán menores para este caso que para el caso en el que <math>a=0.003</math>

==Calibración del parámetro <math>a</math> en base a una experiencia previa==
Ahora vamos a intentar calibrar el párametro <math>a</math> para que el número de infectados tenga el máximo en el quinto día, para ello escribimos el siguiente código:

{{matlab|codigo=
%Aplicamos el Método de Heunn

%Datos
t0=0;
tN=40;
h=0.1;
t=t0:h:tN;
a=0.0005:0.0001:0.002;
b=0.3;
c=0.01;
N1=round((tN-t0)/h);
N2=length(a);
Tmax=zeros(1,N2);

%Vectores iniciales de Población e Infectados
y=zeros(2,N1+1);

%Condiciones iniciales
y(1,1)=1600;
y(2,1)=40;

%Bucle que nos calcula todas las posibles soluciones
for j=1:N2
for i=1:N1
K1S=-a(j)*y(1,i)*y(2,i);
K2S=-a(j)*y(1,i)*y(2,i)+K1S*h;
y(1,i+1)=y(1,i)+h/2*(K1S+K2S);

K1I=y(2,i)*(a(j)*y(1,i)-(b+c));
K2I=y(2,i)*(a(j)*y(1,i)-(b+c))+K1I*h;
y(2,i+1)=y(2,i)+h/2*(K1I+K2I);

if i>1 && y(2,i+1)<y(2,i) && y(2,i)>y(2,i-1)
Tmax(j)=t(i);
end
end
end
distancia=(abs(Tmax-5));
[valormin,posmin]=min(distancia);
a=a(posmin)
}}

Hemos encontrado la solución utilizando un bucle de tipo <math>for</math> dos veces, esto nos permite obtener todas los valores de <math>I</math>, infectados, en el tiempo transcurrido entre [0,40] días para cualquier <math>a</math> en el interválo de [0.0005,0.002]. Una vez tenemos todas las posibilidades aplicamos una condición del tipo <math>if</math> para localizar donde se sitúa el máximo de infectados en los distintos <math>a</math> analizados, después vemos cuál de los máximos esta más próximo al quinto día, para ello generamos un vector <math>distancia</math> donde hacemos el valor absoluto de la resta entre el valor obtenido y el quinto día tras esto conseguimos el valor mas cercano. Ahora analizamos el valor de <math>a</math> para dicha posición obteniendo la solución pedida. En nuestro caso el valor pedida daba '''a = 8.0000e-04''' que era el cuarto <math>a</math> analizado en el bucle.

Modelos epidemiológicos (Grupo 3A)

2016-05-02T19:25:00Z

Trabajocampos: /* Calibración del parámetro a en base a una experiencia previa */

{{ TrabajoED | Modelos Epidemológicos (Grupo 3A) | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED15/16|Curso 2015-16]] |
Ignacio Mollá Carcaño 
Pablo Revuelta Aragón 
David González Hernández 
Jose María García Rodríguez 
Alejandro Martínez Gamonal }}

En el desarrollo de una epidemia, se distinguen dos tipos de individuos. Los que ya han contraído la enfermedad (infectados), que llamaremos <math>I(t)</math> ; y los que son susceptibles de contraerla (sanos), a los que llamaremos <math>S(t)</math>. Donde <math>t</math> es la variable del tiempo. Se presentan dos hipótesis para realizar este estudio:

1. La población de personas infectadas se altera por el fallecimiento o la cura de las mismas. En ambos casos, la tasa de cambio depende del número 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ón entre el número de individuos en ambas clases.

A través del siguiente sistema de ecuaciones diferenciales en función del tiempo observamos las variaciones de ambas poblaciones:

<math>\begin{cases}\frac{dS}{dt}=-aSI\\\frac{dI}{dt}=aSI-bI-cI\end{cases}</math>

Donde <math>a</math>, <math>b</math> y <math>c</math> son parámetros.

* <math>\frac{dS}{dt}</math> : Es la variación de la población susceptible de contraer la enfermedad.

* <math>\frac{dI}{dt}</math> : Es la variación de la población infectada.

* <math>aSI</math> : Es la interacción entre ambas poblaciones que dependen de <math>a</math>.

* <math>bI</math> : Indica los individuos que fallecen. Por lo tanto <math>b</math> podría considerarse como la '''tasa de mortalidad'''.

* <math>cI</math> : Indica los individuos que se han curado. Por lo tanto <math>c</math> podría considerarse como la '''tasa de curación'''.

Para interpretar <math>a</math> observamos el sistema y la segunda hipótesis. Como los individuos que pasan de ser susceptibles a estar infectados son proporcionales a su interacción <math>(SI)</math>, <math>a</math> podría considerarse la '''tasa de contagio''' en la población susceptible.

==Resolución del problema con una sola ecuación diferencial==
Primero vamamos a estudiar el modelo epidemiológico de una población en la que incialmente las 2000 personas que conforman la población se encuentran infectadas. Para hacer el seguimiento de la enfermedad contamos tanto con la constante de proporcionalidad de personas que se curan (b=0.3) y la constante de proporcionalidad del número de personas que fallecen (c=0.01). De esta forma nos quedará el siguiente PVI:

<math> \begin{cases} \frac{dI}{dt} = -0.31 I \\I_{0} = 2000\end{cases} </math>

A continuación vamos a proceder a resolverlo por el método de Euler y trapecio asignandole un tamaño de paso de h=0.1 y veremos cuanto tiempo tardará en reducirse el número de infectados a la cuarta parte.

<math>\Rightarrow</math> Método de Euler
[[Archivo:Ejerciciodoseuler.png|500px|thumb|rigth|Aproximación mediante Euler]]
{{matlab|codigo=

clear all
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
y0=2000;
b=0.3;
c=0.01;

t(1)=t0;
ye(1)=y0; %Vector solución Euler

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de y(i) hasta que y(i)>500

%EULER

while ye(i)>500

ye(i+1)=ye(i)+h*(-b*ye(i)-c*ye(i));
t(i+1)=t(i)+h;
i=i+1;
end

[t',ye'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,ye)
legend('Población infectada(Método de Euler)','Location','best');

}}

<math>\Rightarrow</math> Método del trapecio
[[Archivo:Ejerciciodostrapecio.png|500px|thumb|rigth|Aproximación mediante el método del trapecio]]
{{matlab|codigo=
clear all
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
z0=2000;
b=0.3;
c=0.01;

t(1)=t0;
z(1)=z0; %Vector solución trapecio

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de z(i) hasta que z(i)>500

%TRAPECIO
while z(i)>500
z(i+1)=(z(i)*(1-(h/2)*(b+c)))/(1+h/2*(b+c));%trapecio
t(i+1)=t(i)+h;
i=i+1;
end

[t',z'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,z)
legend('Población infectada(Método del trapecio)','Location','best');

}}

Con todo esto podemos determinar que el número de infectados se reducirá a la cuarta parte en 4.5 días, esto ocurrirá o bien por el fallecimiento o por la curación de dicha enfermedad. Cabe señalar que el número de curados será mayor que el de defunciones, esto se puede ver a simple vista viendo que b>a.
De cara a nuestro estudio de esta enfermedad vamos a ver como afectaría al desarrollo de la misma si introducimos 100 sujetos sanos en la población, tomando una constante de interacción entre personas infectadas y sanas de a=0.003. Procederemos igual que antes a resolverlo por el método de Euler y trapecio con un tamaño de paso de 0.1.
Primero vamos a definir nuestro nuevo PVI:
<math> \begin{cases} \frac{dI}{dt} = -0.01 I \\I_{0} = 2000\end{cases} </math>
<math>\Rightarrow</math> Método de Euler
[[Archivo:Ejerciciotreseuler.png|500px|thumb|rigth|Aproximación mediante Euler para S=100]]
{{matlab|codigo=

%Euler
clear all
%Datos
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%a: Es la constante de proporcionalidad entre los infectados por cada
%interacción
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
y0=2000;
S=100;
a=0.003;
b=0.3;
c=0.01;

t(1)=t0;
y(1)=y0;

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de y(i) hasta que y(i)>500
while y(i)>500

y(i+1)=y(i)+h*(a*S*y(i)-b*y(i)-c*y(i));
t(i+1)=t(i)+h;
i=i+1;
end

[t',y'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,y)
legend('Población infectada para S=100(Método de Euler)','Location','best');

}}

<math>\Rightarrow</math> Método del trapecio
[[Archivo:Ejerciciotrestrapecio.png|500px|thumb|rigth|Aproximación mediante el método del trapecio para S=100]]
{{matlab|codigo=

%Trapecio
clear all
%Datos
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%a: Es la constante de proporcionalidad entre los infectados por cada
%interacción
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
z0=2000;
S=100;
a=0.003;
b=0.3;
c=0.01;

t(1)=t0;
z(1)=z0;

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de z(i) hasta que z(i)>500

%TRAPECIO
while z(i)>500
z(i+1)=(z(i)*(1-(h/2)*(b+c-S*a)))/(1+h/2*(b+c-S*a));
t(i+1)=t(i)+h;
i=i+1;
end

[t',z'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,z)
legend('Población infectada para S=100(Método del trapecio)','Location','best');

}}

Ahora con estos nuevos datos vemos que para que se reduzca a la cuarta parte el número de infectados tendrán que pasar 138.7 días, algo lógico al fijarnos en la ED, ambas son ecuaciones decrecientes pero para el caso de S=100 la pendiente será mucho menor.
Si continuamos haciendo pruebas vemos para valores superiores de s, por ejemplo, S=200 el matlab no nos da ninguna respuesta. Esto ocurre porque al aumentar S la ecuación pasa a ser creciente, por lo que lógicamente si intentas buscar un valor del número de infectados menor al inicial no lo podrás encontrar.

==Resolución del problema completo==
Vamos a resolver el problema completo para ver como evolucionarían las poblaciones en un periodo de 40 días.
Primero supongamos una población infectada inicial de 20 individuos y una población en riesgo de contagio de 800.
Después supondremos el caso de 40 individuos enfermos y 10000 en riesgo de contagio.
Para resolverlo usaremos el método de Euler y más tarde el método de Runge-Kutta de orden cuatro con distintas discretizaciones para ver la influencia de estas.

<math>\Rightarrow</math> Método de Euler
{{matlab|codigo=

%Euler
clear all
%Datos
%Introducimos los valores de las constantes
a=0.003;
b=0.3;
c=0.01;
h=input('Introducir valores de h:');
t0=0;
y0=input('Introducir vector [So,Io]:');
tN=40;
% calculamos los subintervalos
t=t0:h:tN;
N=(tN-t0)/h;
y=zeros(2,length(t));
y(:,1)=y0';
for i= 1:N;
y(:,i+1)=y(:,i)+h*[-a*y(1,i).*y(2,i);a*y(1,i).*y(2,i)-b*y(2,i)-c*y(2,i)];
end

% grafico
plot(t,y)
% vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h
legend('Población sana','Location','best','Población enferma','Location','best');

}}

<math>\Rightarrow</math> Método de Runge-Kutta

{{matlab|codigo=

%Runge-Kutta
clear all

%Datos
S0=input('Introducir valor inicial SO: ');
I0=input('Introducir valor inicial IO: ');
a=0.003;
b=0.3;
c=0.01;
%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector solucion
y0=[S0;I0];
y=zeros(2,N+1);
y(:,1)=y0;
for i=1:N
K1S=-a*y(1,i)*y(2,i);
K2S=-a*(y(1,i)+(1/2)*K1S*h)*(y(2,i)+(1/2)*K1S*h);
K3S=-a*(y(1,i)+(1/2)*K2S*h)*(y(2,i)+(1/2)*K2S*h);
K4S=-a*(y(1,i)+K3S*h)*(y(2,i)+K3S*h);
y(1,i+1)=y(1,i)+(h/6)*(K1S+2*K2S+2*K3S+K4S);
K1I=a*y(1,i)*y(2,i)-b*y(2,i)-c*y(2,i);
K2I=a*(y(1,i)+(1/2)*K1I*h)*(y(2,i)+(1/2)*K1I*h)-b*(y(2,i)+(1/2)*K1I*h)-c*(y(2,i)+(1/2)*K1I*h);
K3I=a*(y(1,i)+(1/2)*K2I*h)*(y(2,i)+(1/2)*K2I*h)-b*(y(2,i)+(1/2)*K2I*h)-c*(y(2,i)+(1/2)*K2I*h);
K4I=a*(y(1,i)+K3I*h)*(y(2,i)+K3I*h)-b*(y(2,i)+K3I*h)-c*(y(2,i)+K3I*h);
y(2,i+1)=y(2,i)+(h/6)*(K1I+2*K2I+2*K3I+K4I);
end

%Tabla de resutados
[t',y(1,:)',y(2,:)']
%Gráfico
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'k')
legend('Poblacion sana','Poblacion enferma','Location','best');
hold off

% vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h

}}

La dificultad que hay en el uso del método del trapecio o cualquier otro método implícito está en que las incógnitas de nuestro sistema <math>(S,I)</math> quedan implícitas en ambas ecuaciones y necesitan que se despejen. Al comenzar a realizar cálculos nos damos cuenta que el despeje es imposible, ya que las dos variables incógnita aparecen multiplicándose entre sí, quedando una variable dependiendo de la otra siempre.

===Situación inicial de <math>(S _{0} , I_{0} )=(800,20)</math>===
En estos programas, sustituyendo las <math>h</math> para 0.1, 0.01, 0.001 y 0.0001 obtendremos los siguientes resultados:

-Para <math>h=0.1</math>

[[Archivo:Tborraré.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.1rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

Como se puede ver la población sana se infecta rápidamente y a partir del día 15 prácticamente toda la población ha fallecido.
* '''Método de Euler''': nuestro programa nos indica que en este caso el momento en el que los enfermos son máximos se produce al '''tercer día (2.8)''' y el '''número de enfermos en ese momento es 517'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce al '''tercer día (2.6)''' y el '''número de enfermos en ese momento es 650'''.

-Para <math>h=0.01</math>
[[Archivo:Tborraré2.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.01rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

En el gráfico apenas se encuentra diferencia, lo cual es lógico porque es el mismo problema, pero si sobrepusiéramos las gráficas, lograríamos ver como se va ajustando a la realidad en función del aumento de h.
* '''Método de Euler''': en este caso el momento en el que los enfermos son máximos es también en el '''tercer día (2.7)''' y el número de '''enfermos es 506'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce al '''tercer día (2.69)''' y el número de '''enfermos en ese momento es 517'''.

-Para <math>h=0.001</math>
[[Archivo:Tborraré3.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.001rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': en este caso el momento en el que los enfermos son máximos se produce el '''tercer día (2.696)''' con un valor de '''505 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''tercer día (2.694)''' con un valor de '''506 enfermos'''.

-Para <math>h=0.0001</math>
[[Archivo:Tborraré4.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.0001rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': al igual que antes, se repite el momento en el que los enfermos son máximos en el '''tercer día (2.695)''' con un valor de '''505 enfermos'''.
* '''Método de Runge-Kutta''': se repite el momento en el que los enfermos son máximos en el '''tercer día (2.6948)''' con un valor de '''505 enfermos'''.

Observamos en las distintas gráficas obtenidas, que según vamos empleando un tamaño de paso <math>h</math> más pequeño, nos vamos acercando cada vez más a la población real que pretendemos reflejar.
Nos damos cuenta de que las gráficas tienen bastante lógica, pues a la vez que desciende el número de población sana va aumentando la población infectada.
También observamos que para el método de Runge-Kutta cuando utilizamos tamaño de paso <math>h=0.1</math> pierde precisión en relación a cuando utilizamos tamaños de paso más pequeños.

===Situación inicial de <math>(S _{0} , I_{0} )=(10000,40)</math>===
En estos programas, sustituyendo las <math>h</math> para 0.1, 0.01, 0.001 y 0.0001 obtendremos los siguientes resultados:

-Para <math>h=0.1</math>

[[Archivo:H0.11000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.5)''' con un valor de '''12600 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.4)''' con un valor de '''2.6e74 enfermos'''.

-Para <math>h=0.01</math>
[[Archivo:H0.0110000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.0110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.35)''' con un valor de '''9521 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.31)''' con un valor de '''13550 enfermos'''.

-Para <math>h=0.001</math>
[[Archivo:H0.00110000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.00110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.342)''' con un valor de '''9470 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.338)''' con un valor de '''9758 enfermos'''.

-Para <math>h=0.0001</math>
[[Archivo:H0.000110000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.000110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.3411)''' con un valor de '''9465 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.3407)''' con un valor de '''9493 enfermos'''.

Para este segundo caso, nos damos cuenta que al utilizar tamaño de paso <math>h=0.1</math> , llegamos a un resultado sin sentido, pues nos da valores negativos y es imposible que una población tenga tales valores.
En las demás gráficas, los resultados obtenidos no parecen cuadrar con la realidad. Esto podría ser debido a que la alta diferencia entre los valores iniciales <math>(S _{0} , I_{0} )=(10000,40)</math> no permiten un buen cálculo de su variación respecto del tiempo.
El número inicial de infectados es muy elevado y al disminuir el tamaño de paso, el programa aproxima mejor la solución.

==Resolución del problema completo suponiendo <math> a </math> una función variable con el tiempo==
Para el siguiente caso supondremos que el coeficiente <math> a </math> (tasa de contagio en población sana) es una función dada por :

<math> a(t)= \frac{0.003}{(1+t)} </math>
Suponiendo además las condiciones iniciales <math> (S_{0},I_{0})=(1640,40) </math> , vamos a dibujar las gráficas de infectados y sanos por el método de Heun, tomando un tamaño de paso <math> h=0.01 </math>
[[Archivo: H0.01heunn.png|1000px|thumb|rigth| Población infectada-Población sana para a=a(t)]]
{{matlab|codigo=
%Heun. Coeficiente variable a(t)

%Datos
S0=input('Introducir valor inicial SO: ');
I0=input('Introducir valor inicial IO: ');
b=0.3;
c=0.01;
%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector solucion
y0=[S0;I0];
y=zeros(2,N+1);
y(:,1)=y0;
for i=1:N
K1S=-(0.003/(1+t(i)))*y(1,i)*y(2,i);
K2S=-(0.003/(1+t(i)+h))*(y(1,i)+K1S*h)*(y(2,i)+K1S*h);
y(1,i+1)=y(1,i)+(h/2)*(K1S+K2S);
K1I=(0.003/(1+t(i)))*y(1,i)*y(2,i)-b*y(2,i)-c*y(2,i);
K2I=(0.003/(1+t(i)+h))*(y(1,i)+K1I*h)*(y(2,i)+K1I*h)-b*(y(2,i)+K1I*h)-c*(y(2,i)+K1I*h);
y(2,i+1)=y(2,i)+(h/2)*(K1I+K2I);
end

%Tabla de resutados
[t',y(1,:)',y(2,:)']
%Graficos
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'k')
title('Método de Heun');
legend('Poblacion sana','Población infectada','Location','best');
hold off

% Vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h

}}

Como se puede ver la población sana se infecta rápidamente y a partir del día 15 casi toda la población ha fallecido. Al segundo día <math>(t=2,15)</math> observamos que el número de infectados es de 994 individuos <math>(I=994,3)</math>, siendo este el máximo que se alcanzará.

A medida que va pasando el tiempo el coeficiente <math> a(t) </math> (tasa de contagio en población sana) va a ir decreciendo poco a poco hasta aproximarse a cero según podemos ver en el siguiente gráfico:
<math> \lim_{t\to \infty } a(t)=0</math>
[[Archivo: A(t).png|450px|thumb|right|a(t) para t=[0,40]]]
{{matlab|codigo=

%tasa de contagio en poblacion sana a(t)

%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector a
a=zeros(1,N+1);
for i=1:N+1
a(i)=0.003/(1+t(i));
end

%Grafica
plot(t,a)

}}

[[Archivo: Comparacionheun.png|border|1500px|center|frame|Poblaciones según la tasa de contagio a(t)]]

Ambas poblaciones parten del mismo inicio (para <math>t=0</math> las poblaciones son iguales). Podemos observar que la población sana (susceptible de contraer la enfermedad) en el caso en el que la tasa de contagio en población sana <math>(a)</math> es constante <math>(a=0.003)</math> se ve reducida antes que en el caso en el que depende del tiempo <math>(a=a(t))</math>.

<math>a(t)</math> disminuye a la vez que transcurre el tiempo, por lo que en general las pendientes tanto de la población sana como de la infectada serán menores para este caso que para el caso en el que <math>a=0.003</math>

==Calibración del parámetro <math>a</math> en base a una experiencia previa==
Ahora vamos a intentar calibrar el párametro <math>a</math> para que el número de infectados tenga el máximo en el quinto día, para ello escribimos el siguiente código:

{{matlab|codigo=
%Aplicamos el Método de Heunn

%Datos
t0=0;
tN=40;
h=0.1;
t=t0:h:tN;
a=0.0005:0.0001:0.002;
b=0.3;
c=0.01;
N1=round((tN-t0)/h);
N2=length(a);
Tmax=zeros(1,N2);

%Vectores iniciales de Población e Infectados
y=zeros(2,N1+1);

%Condiciones iniciales
y(1,1)=1600;
y(2,1)=40;

%Bucle que nos calcula todas las posibles soluciones
for j=1:N2
for i=1:N1
K1S=-a(j)*y(1,i)*y(2,i);
K2S=-a(j)*y(1,i)*y(2,i)+K1S*h;
y(1,i+1)=y(1,i)+h/2*(K1S+K2S);

K1I=y(2,i)*(a(j)*y(1,i)-(b+c));
K2I=y(2,i)*(a(j)*y(1,i)-(b+c))+K1I*h;
y(2,i+1)=y(2,i)+h/2*(K1I+K2I);

if i>1 && y(2,i+1)<y(2,i) && y(2,i)>y(2,i-1)
Tmax(j)=t(i);
end
end
end
distancia=(abs(Tmax-5));
[valormin,posmin]=min(distancia);
a=a(posmin)
}}

Hemos encontrado la solución utilizando un bucle de tipo <math>for</math> dos veces, esto nos permite obtener todas los valores de <math>I</math>, infectados, en el tiempo transcurrido entre [0,40] días para cualquier <math>a</math> en el interválo de [0.0005,0.002]. Una vez tenemos todas las posibilidades aplicamos una condición del tipo <math>if</math> para localizar donde se sitúa el máximo de infectados en los distintos <math>a</math> analizados, después vemos cuál de los máximos esta más próximo al quinto día, para ello generamos un vector <math>distancia</math> donde hacemos el valor absoluto de la resta entre el valor obtenido y el quinto día traesto conseguimos el valor mas próximo. Ahora analizamos el valor de <math>a</math> para dicha posición obteniendo la solución pedida. en nuestro caso el valor pedida daba '''a = 8.0000e-04''' que era el cuarto <math>a</math> analizado en el bucle.

Modelos epidemiológicos (Grupo 3A)

2016-05-02T19:23:36Z

Trabajocampos: /* Calibración del parámetro a en base a una experiencia previa */

{{ TrabajoED | Modelos Epidemológicos (Grupo 3A) | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED15/16|Curso 2015-16]] |
Ignacio Mollá Carcaño 
Pablo Revuelta Aragón 
David González Hernández 
Jose María García Rodríguez 
Alejandro Martínez Gamonal }}

En el desarrollo de una epidemia, se distinguen dos tipos de individuos. Los que ya han contraído la enfermedad (infectados), que llamaremos <math>I(t)</math> ; y los que son susceptibles de contraerla (sanos), a los que llamaremos <math>S(t)</math>. Donde <math>t</math> es la variable del tiempo. Se presentan dos hipótesis para realizar este estudio:

1. La población de personas infectadas se altera por el fallecimiento o la cura de las mismas. En ambos casos, la tasa de cambio depende del número 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ón entre el número de individuos en ambas clases.

A través del siguiente sistema de ecuaciones diferenciales en función del tiempo observamos las variaciones de ambas poblaciones:

<math>\begin{cases}\frac{dS}{dt}=-aSI\\\frac{dI}{dt}=aSI-bI-cI\end{cases}</math>

Donde <math>a</math>, <math>b</math> y <math>c</math> son parámetros.

* <math>\frac{dS}{dt}</math> : Es la variación de la población susceptible de contraer la enfermedad.

* <math>\frac{dI}{dt}</math> : Es la variación de la población infectada.

* <math>aSI</math> : Es la interacción entre ambas poblaciones que dependen de <math>a</math>.

* <math>bI</math> : Indica los individuos que fallecen. Por lo tanto <math>b</math> podría considerarse como la '''tasa de mortalidad'''.

* <math>cI</math> : Indica los individuos que se han curado. Por lo tanto <math>c</math> podría considerarse como la '''tasa de curación'''.

Para interpretar <math>a</math> observamos el sistema y la segunda hipótesis. Como los individuos que pasan de ser susceptibles a estar infectados son proporcionales a su interacción <math>(SI)</math>, <math>a</math> podría considerarse la '''tasa de contagio''' en la población susceptible.

==Resolución del problema con una sola ecuación diferencial==
Primero vamamos a estudiar el modelo epidemiológico de una población en la que incialmente las 2000 personas que conforman la población se encuentran infectadas. Para hacer el seguimiento de la enfermedad contamos tanto con la constante de proporcionalidad de personas que se curan (b=0.3) y la constante de proporcionalidad del número de personas que fallecen (c=0.01). De esta forma nos quedará el siguiente PVI:

<math> \begin{cases} \frac{dI}{dt} = -0.31 I \\I_{0} = 2000\end{cases} </math>

A continuación vamos a proceder a resolverlo por el método de Euler y trapecio asignandole un tamaño de paso de h=0.1 y veremos cuanto tiempo tardará en reducirse el número de infectados a la cuarta parte.

<math>\Rightarrow</math> Método de Euler
[[Archivo:Ejerciciodoseuler.png|500px|thumb|rigth|Aproximación mediante Euler]]
{{matlab|codigo=

clear all
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
y0=2000;
b=0.3;
c=0.01;

t(1)=t0;
ye(1)=y0; %Vector solución Euler

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de y(i) hasta que y(i)>500

%EULER

while ye(i)>500

ye(i+1)=ye(i)+h*(-b*ye(i)-c*ye(i));
t(i+1)=t(i)+h;
i=i+1;
end

[t',ye'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,ye)
legend('Población infectada(Método de Euler)','Location','best');

}}

<math>\Rightarrow</math> Método del trapecio
[[Archivo:Ejerciciodostrapecio.png|500px|thumb|rigth|Aproximación mediante el método del trapecio]]
{{matlab|codigo=
clear all
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
z0=2000;
b=0.3;
c=0.01;

t(1)=t0;
z(1)=z0; %Vector solución trapecio

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de z(i) hasta que z(i)>500

%TRAPECIO
while z(i)>500
z(i+1)=(z(i)*(1-(h/2)*(b+c)))/(1+h/2*(b+c));%trapecio
t(i+1)=t(i)+h;
i=i+1;
end

[t',z'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,z)
legend('Población infectada(Método del trapecio)','Location','best');

}}

Con todo esto podemos determinar que el número de infectados se reducirá a la cuarta parte en 4.5 días, esto ocurrirá o bien por el fallecimiento o por la curación de dicha enfermedad. Cabe señalar que el número de curados será mayor que el de defunciones, esto se puede ver a simple vista viendo que b>a.
De cara a nuestro estudio de esta enfermedad vamos a ver como afectaría al desarrollo de la misma si introducimos 100 sujetos sanos en la población, tomando una constante de interacción entre personas infectadas y sanas de a=0.003. Procederemos igual que antes a resolverlo por el método de Euler y trapecio con un tamaño de paso de 0.1.
Primero vamos a definir nuestro nuevo PVI:
<math> \begin{cases} \frac{dI}{dt} = -0.01 I \\I_{0} = 2000\end{cases} </math>
<math>\Rightarrow</math> Método de Euler
[[Archivo:Ejerciciotreseuler.png|500px|thumb|rigth|Aproximación mediante Euler para S=100]]
{{matlab|codigo=

%Euler
clear all
%Datos
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%a: Es la constante de proporcionalidad entre los infectados por cada
%interacción
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
y0=2000;
S=100;
a=0.003;
b=0.3;
c=0.01;

t(1)=t0;
y(1)=y0;

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de y(i) hasta que y(i)>500
while y(i)>500

y(i+1)=y(i)+h*(a*S*y(i)-b*y(i)-c*y(i));
t(i+1)=t(i)+h;
i=i+1;
end

[t',y'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,y)
legend('Población infectada para S=100(Método de Euler)','Location','best');

}}

<math>\Rightarrow</math> Método del trapecio
[[Archivo:Ejerciciotrestrapecio.png|500px|thumb|rigth|Aproximación mediante el método del trapecio para S=100]]
{{matlab|codigo=

%Trapecio
clear all
%Datos
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%a: Es la constante de proporcionalidad entre los infectados por cada
%interacción
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
z0=2000;
S=100;
a=0.003;
b=0.3;
c=0.01;

t(1)=t0;
z(1)=z0;

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de z(i) hasta que z(i)>500

%TRAPECIO
while z(i)>500
z(i+1)=(z(i)*(1-(h/2)*(b+c-S*a)))/(1+h/2*(b+c-S*a));
t(i+1)=t(i)+h;
i=i+1;
end

[t',z'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,z)
legend('Población infectada para S=100(Método del trapecio)','Location','best');

}}

Ahora con estos nuevos datos vemos que para que se reduzca a la cuarta parte el número de infectados tendrán que pasar 138.7 días, algo lógico al fijarnos en la ED, ambas son ecuaciones decrecientes pero para el caso de S=100 la pendiente será mucho menor.
Si continuamos haciendo pruebas vemos para valores superiores de s, por ejemplo, S=200 el matlab no nos da ninguna respuesta. Esto ocurre porque al aumentar S la ecuación pasa a ser creciente, por lo que lógicamente si intentas buscar un valor del número de infectados menor al inicial no lo podrás encontrar.

==Resolución del problema completo==
Vamos a resolver el problema completo para ver como evolucionarían las poblaciones en un periodo de 40 días.
Primero supongamos una población infectada inicial de 20 individuos y una población en riesgo de contagio de 800.
Después supondremos el caso de 40 individuos enfermos y 10000 en riesgo de contagio.
Para resolverlo usaremos el método de Euler y más tarde el método de Runge-Kutta de orden cuatro con distintas discretizaciones para ver la influencia de estas.

<math>\Rightarrow</math> Método de Euler
{{matlab|codigo=

%Euler
clear all
%Datos
%Introducimos los valores de las constantes
a=0.003;
b=0.3;
c=0.01;
h=input('Introducir valores de h:');
t0=0;
y0=input('Introducir vector [So,Io]:');
tN=40;
% calculamos los subintervalos
t=t0:h:tN;
N=(tN-t0)/h;
y=zeros(2,length(t));
y(:,1)=y0';
for i= 1:N;
y(:,i+1)=y(:,i)+h*[-a*y(1,i).*y(2,i);a*y(1,i).*y(2,i)-b*y(2,i)-c*y(2,i)];
end

% grafico
plot(t,y)
% vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h
legend('Población sana','Location','best','Población enferma','Location','best');

}}

<math>\Rightarrow</math> Método de Runge-Kutta

{{matlab|codigo=

%Runge-Kutta
clear all

%Datos
S0=input('Introducir valor inicial SO: ');
I0=input('Introducir valor inicial IO: ');
a=0.003;
b=0.3;
c=0.01;
%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector solucion
y0=[S0;I0];
y=zeros(2,N+1);
y(:,1)=y0;
for i=1:N
K1S=-a*y(1,i)*y(2,i);
K2S=-a*(y(1,i)+(1/2)*K1S*h)*(y(2,i)+(1/2)*K1S*h);
K3S=-a*(y(1,i)+(1/2)*K2S*h)*(y(2,i)+(1/2)*K2S*h);
K4S=-a*(y(1,i)+K3S*h)*(y(2,i)+K3S*h);
y(1,i+1)=y(1,i)+(h/6)*(K1S+2*K2S+2*K3S+K4S);
K1I=a*y(1,i)*y(2,i)-b*y(2,i)-c*y(2,i);
K2I=a*(y(1,i)+(1/2)*K1I*h)*(y(2,i)+(1/2)*K1I*h)-b*(y(2,i)+(1/2)*K1I*h)-c*(y(2,i)+(1/2)*K1I*h);
K3I=a*(y(1,i)+(1/2)*K2I*h)*(y(2,i)+(1/2)*K2I*h)-b*(y(2,i)+(1/2)*K2I*h)-c*(y(2,i)+(1/2)*K2I*h);
K4I=a*(y(1,i)+K3I*h)*(y(2,i)+K3I*h)-b*(y(2,i)+K3I*h)-c*(y(2,i)+K3I*h);
y(2,i+1)=y(2,i)+(h/6)*(K1I+2*K2I+2*K3I+K4I);
end

%Tabla de resutados
[t',y(1,:)',y(2,:)']
%Gráfico
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'k')
legend('Poblacion sana','Poblacion enferma','Location','best');
hold off

% vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h

}}

La dificultad que hay en el uso del método del trapecio o cualquier otro método implícito está en que las incógnitas de nuestro sistema <math>(S,I)</math> quedan implícitas en ambas ecuaciones y necesitan que se despejen. Al comenzar a realizar cálculos nos damos cuenta que el despeje es imposible, ya que las dos variables incógnita aparecen multiplicándose entre sí, quedando una variable dependiendo de la otra siempre.

===Situación inicial de <math>(S _{0} , I_{0} )=(800,20)</math>===
En estos programas, sustituyendo las <math>h</math> para 0.1, 0.01, 0.001 y 0.0001 obtendremos los siguientes resultados:

-Para <math>h=0.1</math>

[[Archivo:Tborraré.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.1rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

Como se puede ver la población sana se infecta rápidamente y a partir del día 15 prácticamente toda la población ha fallecido.
* '''Método de Euler''': nuestro programa nos indica que en este caso el momento en el que los enfermos son máximos se produce al '''tercer día (2.8)''' y el '''número de enfermos en ese momento es 517'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce al '''tercer día (2.6)''' y el '''número de enfermos en ese momento es 650'''.

-Para <math>h=0.01</math>
[[Archivo:Tborraré2.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.01rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

En el gráfico apenas se encuentra diferencia, lo cual es lógico porque es el mismo problema, pero si sobrepusiéramos las gráficas, lograríamos ver como se va ajustando a la realidad en función del aumento de h.
* '''Método de Euler''': en este caso el momento en el que los enfermos son máximos es también en el '''tercer día (2.7)''' y el número de '''enfermos es 506'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce al '''tercer día (2.69)''' y el número de '''enfermos en ese momento es 517'''.

-Para <math>h=0.001</math>
[[Archivo:Tborraré3.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.001rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': en este caso el momento en el que los enfermos son máximos se produce el '''tercer día (2.696)''' con un valor de '''505 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''tercer día (2.694)''' con un valor de '''506 enfermos'''.

-Para <math>h=0.0001</math>
[[Archivo:Tborraré4.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.0001rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': al igual que antes, se repite el momento en el que los enfermos son máximos en el '''tercer día (2.695)''' con un valor de '''505 enfermos'''.
* '''Método de Runge-Kutta''': se repite el momento en el que los enfermos son máximos en el '''tercer día (2.6948)''' con un valor de '''505 enfermos'''.

Observamos en las distintas gráficas obtenidas, que según vamos empleando un tamaño de paso <math>h</math> más pequeño, nos vamos acercando cada vez más a la población real que pretendemos reflejar.
Nos damos cuenta de que las gráficas tienen bastante lógica, pues a la vez que desciende el número de población sana va aumentando la población infectada.
También observamos que para el método de Runge-Kutta cuando utilizamos tamaño de paso <math>h=0.1</math> pierde precisión en relación a cuando utilizamos tamaños de paso más pequeños.

===Situación inicial de <math>(S _{0} , I_{0} )=(10000,40)</math>===
En estos programas, sustituyendo las <math>h</math> para 0.1, 0.01, 0.001 y 0.0001 obtendremos los siguientes resultados:

-Para <math>h=0.1</math>

[[Archivo:H0.11000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.5)''' con un valor de '''12600 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.4)''' con un valor de '''2.6e74 enfermos'''.

-Para <math>h=0.01</math>
[[Archivo:H0.0110000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.0110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.35)''' con un valor de '''9521 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.31)''' con un valor de '''13550 enfermos'''.

-Para <math>h=0.001</math>
[[Archivo:H0.00110000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.00110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.342)''' con un valor de '''9470 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.338)''' con un valor de '''9758 enfermos'''.

-Para <math>h=0.0001</math>
[[Archivo:H0.000110000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.000110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.3411)''' con un valor de '''9465 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.3407)''' con un valor de '''9493 enfermos'''.

Para este segundo caso, nos damos cuenta que al utilizar tamaño de paso <math>h=0.1</math> , llegamos a un resultado sin sentido, pues nos da valores negativos y es imposible que una población tenga tales valores.
En las demás gráficas, los resultados obtenidos no parecen cuadrar con la realidad. Esto podría ser debido a que la alta diferencia entre los valores iniciales <math>(S _{0} , I_{0} )=(10000,40)</math> no permiten un buen cálculo de su variación respecto del tiempo.
El número inicial de infectados es muy elevado y al disminuir el tamaño de paso, el programa aproxima mejor la solución.

==Resolución del problema completo suponiendo <math> a </math> una función variable con el tiempo==
Para el siguiente caso supondremos que el coeficiente <math> a </math> (tasa de contagio en población sana) es una función dada por :

<math> a(t)= \frac{0.003}{(1+t)} </math>
Suponiendo además las condiciones iniciales <math> (S_{0},I_{0})=(1640,40) </math> , vamos a dibujar las gráficas de infectados y sanos por el método de Heun, tomando un tamaño de paso <math> h=0.01 </math>
[[Archivo: H0.01heunn.png|1000px|thumb|rigth| Población infectada-Población sana para a=a(t)]]
{{matlab|codigo=
%Heun. Coeficiente variable a(t)

%Datos
S0=input('Introducir valor inicial SO: ');
I0=input('Introducir valor inicial IO: ');
b=0.3;
c=0.01;
%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector solucion
y0=[S0;I0];
y=zeros(2,N+1);
y(:,1)=y0;
for i=1:N
K1S=-(0.003/(1+t(i)))*y(1,i)*y(2,i);
K2S=-(0.003/(1+t(i)+h))*(y(1,i)+K1S*h)*(y(2,i)+K1S*h);
y(1,i+1)=y(1,i)+(h/2)*(K1S+K2S);
K1I=(0.003/(1+t(i)))*y(1,i)*y(2,i)-b*y(2,i)-c*y(2,i);
K2I=(0.003/(1+t(i)+h))*(y(1,i)+K1I*h)*(y(2,i)+K1I*h)-b*(y(2,i)+K1I*h)-c*(y(2,i)+K1I*h);
y(2,i+1)=y(2,i)+(h/2)*(K1I+K2I);
end

%Tabla de resutados
[t',y(1,:)',y(2,:)']
%Graficos
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'k')
title('Método de Heun');
legend('Poblacion sana','Población infectada','Location','best');
hold off

% Vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h

}}

Como se puede ver la población sana se infecta rápidamente y a partir del día 15 casi toda la población ha fallecido. Al segundo día <math>(t=2,15)</math> observamos que el número de infectados es de 994 individuos <math>(I=994,3)</math>, siendo este el máximo que se alcanzará.

A medida que va pasando el tiempo el coeficiente <math> a(t) </math> (tasa de contagio en población sana) va a ir decreciendo poco a poco hasta aproximarse a cero según podemos ver en el siguiente gráfico:
<math> \lim_{t\to \infty } a(t)=0</math>
[[Archivo: A(t).png|450px|thumb|right|a(t) para t=[0,40]]]
{{matlab|codigo=

%tasa de contagio en poblacion sana a(t)

%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector a
a=zeros(1,N+1);
for i=1:N+1
a(i)=0.003/(1+t(i));
end

%Grafica
plot(t,a)

}}

[[Archivo: Comparacionheun.png|border|1500px|center|frame|Poblaciones según la tasa de contagio a(t)]]

Ambas poblaciones parten del mismo inicio (para <math>t=0</math> las poblaciones son iguales). Podemos observar que la población sana (susceptible de contraer la enfermedad) en el caso en el que la tasa de contagio en población sana <math>(a)</math> es constante <math>(a=0.003)</math> se ve reducida antes que en el caso en el que depende del tiempo <math>(a=a(t))</math>.

<math>a(t)</math> disminuye a la vez que transcurre el tiempo, por lo que en general las pendientes tanto de la población sana como de la infectada serán menores para este caso que para el caso en el que <math>a=0.003</math>

==Calibración del parámetro <math>a</math> en base a una experiencia previa==
Ahora vamos a intentar calibrar el párametro <math>a</math> para que el número de infectados tenga el máximo en el quinto día, para ello escribimos el siguiente código:

{{matlab|codigo=
%Aplicamos el Método de Heunn

%Datos
t0=0;
tN=40;
h=0.1;
t=t0:h:tN;
a=0.0005:0.0001:0.002;
b=0.3;
c=0.01;
N1=round((tN-t0)/h);
N2=length(a);
Tmax=zeros(1,N2);

%Vectores iniciales de Población e Infectados
y=zeros(2,N1+1);

%Condiciones iniciales
y(1,1)=1600;
y(2,1)=40;

%Bucle que nos calcula todas las posibles soluciones
for j=1:N2
for i=1:N1
K1S=-a(j)*y(1,i)*y(2,i);
K2S=-a(j)*y(1,i)*y(2,i)+K1S*h;
y(1,i+1)=y(1,i)+h/2*(K1S+K2S);

K1I=y(2,i)*(a(j)*y(1,i)-(b+c));
K2I=y(2,i)*(a(j)*y(1,i)-(b+c))+K1I*h;
y(2,i+1)=y(2,i)+h/2*(K1I+K2I);

if i>1 && y(2,i+1)<y(2,i) && y(2,i)>y(2,i-1)
Tmax(j)=t(i);
end
end
end
distancia=(abs(Tmax-5));
[valormin,posmin]=min(distancia);
a=a(posmin)
}}

Hemos encontrado la solución utilizando un bucle de tipo <math>for</math> dos veces, esto nos permite obtener todas los valores de <math>I</math>, infectados, en el tiempo transcurrido entre [0,40] días para cualquier <math>a</math> en el interválo de [0.0005,0.002]. Una vez tenemos todas las posibilidades aplicamos una condición del tipo <math>if</math> para localizar donde se sitúa el máximo de infectados de los distintos <math>a</math>, después vemos cuál de los máximos esta más próximo al quinto día, para ello generamos un vector <math>distancia</math> donde hacemos el valor absoluto de la resta entre el valor obtenido y el quinto día traesto conseguimos el valor mas próximo. Ahora analizamos el valor de <math>a</math> para dicha posición obteniendo la solución pedida. en nuestro caso el valor pedida daba '''a = 8.0000e-04''' que era el cuarto <math>a</math> analizado en el bucle.

Modelos epidemiológicos (Grupo 3A)

2016-05-02T18:41:53Z

Trabajocampos: /* Calibración del parámetro a en base a una experiencia previa */

{{ TrabajoED | Modelos Epidemológicos (Grupo 3A) | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED15/16|Curso 2015-16]] |
Ignacio Mollá Carcaño 
Pablo Revuelta Aragón 
David González Hernández 
Jose María García Rodríguez 
Alejandro Martínez Gamonal }}

En el desarrollo de una epidemia, se distinguen dos tipos de individuos. Los que ya han contraído la enfermedad (infectados), que llamaremos <math>I(t)</math> ; y los que son susceptibles de contraerla (sanos), a los que llamaremos <math>S(t)</math>. Donde <math>t</math> es la variable del tiempo. Se presentan dos hipótesis para realizar este estudio:

1. La población de personas infectadas se altera por el fallecimiento o la cura de las mismas. En ambos casos, la tasa de cambio depende del número 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ón entre el número de individuos en ambas clases.

A través del siguiente sistema de ecuaciones diferenciales en función del tiempo observamos las variaciones de ambas poblaciones:

<math>\begin{cases}\frac{dS}{dt}=-aSI\\\frac{dI}{dt}=aSI-bI-cI\end{cases}</math>

Donde <math>a</math>, <math>b</math> y <math>c</math> son parámetros.

* <math>\frac{dS}{dt}</math> : Es la variación de la población susceptible de contraer la enfermedad.

* <math>\frac{dI}{dt}</math> : Es la variación de la población infectada.

* <math>aSI</math> : Es la interacción entre ambas poblaciones que dependen de <math>a</math>.

* <math>bI</math> : Indica los individuos que fallecen. Por lo tanto <math>b</math> podría considerarse como la '''tasa de mortalidad'''.

* <math>cI</math> : Indica los individuos que se han curado. Por lo tanto <math>c</math> podría considerarse como la '''tasa de curación'''.

Para interpretar <math>a</math> observamos el sistema y la segunda hipótesis. Como los individuos que pasan de ser susceptibles a estar infectados son proporcionales a su interacción <math>(SI)</math>, <math>a</math> podría considerarse la '''tasa de contagio''' en la población susceptible.

==Resolución del problema con una sola ecuación diferencial==
Primero vamamos a estudiar el modelo epidemiológico de una población en la que incialmente las 2000 personas que conforman la población se encuentran infectadas. Para hacer el seguimiento de la enfermedad contamos tanto con la constante de proporcionalidad de personas que se curan (b=0.3) y la constante de proporcionalidad del número de personas que fallecen (c=0.01). De esta forma nos quedará el siguiente PVI:

<math> \begin{cases} \frac{dI}{dt} = -0.31 I \\I_{0} = 2000\end{cases} </math>

A continuación vamos a proceder a resolverlo por el método de Euler y trapecio asignandole un tamaño de paso de h=0.1 y veremos cuanto tiempo tardará en reducirse el número de infectados a la cuarta parte.

<math>\Rightarrow</math> Método de Euler
[[Archivo:Ejerciciodoseuler.png|500px|thumb|rigth|Aproximación mediante Euler]]
{{matlab|codigo=

clear all
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
y0=2000;
b=0.3;
c=0.01;

t(1)=t0;
ye(1)=y0; %Vector solución Euler

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de y(i) hasta que y(i)>500

%EULER

while ye(i)>500

ye(i+1)=ye(i)+h*(-b*ye(i)-c*ye(i));
t(i+1)=t(i)+h;
i=i+1;
end

[t',ye'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,ye)
legend('Población infectada(Método de Euler)','Location','best');

}}

<math>\Rightarrow</math> Método del trapecio
[[Archivo:Ejerciciodostrapecio.png|500px|thumb|rigth|Aproximación mediante el método del trapecio]]
{{matlab|codigo=
clear all
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
z0=2000;
b=0.3;
c=0.01;

t(1)=t0;
z(1)=z0; %Vector solución trapecio

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de z(i) hasta que z(i)>500

%TRAPECIO
while z(i)>500
z(i+1)=(z(i)*(1-(h/2)*(b+c)))/(1+h/2*(b+c));%trapecio
t(i+1)=t(i)+h;
i=i+1;
end

[t',z'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,z)
legend('Población infectada(Método del trapecio)','Location','best');

}}

Con todo esto podemos determinar que el número de infectados se reducirá a la cuarta parte en 4.5 días, esto ocurrirá o bien por el fallecimiento o por la curación de dicha enfermedad. Cabe señalar que el número de curados será mayor que el de defunciones, esto se puede ver a simple vista viendo que b>a.
De cara a nuestro estudio de esta enfermedad vamos a ver como afectaría al desarrollo de la misma si introducimos 100 sujetos sanos en la población, tomando una constante de interacción entre personas infectadas y sanas de a=0.003. Procederemos igual que antes a resolverlo por el método de Euler y trapecio con un tamaño de paso de 0.1.
Primero vamos a definir nuestro nuevo PVI:
<math> \begin{cases} \frac{dI}{dt} = -0.01 I \\I_{0} = 2000\end{cases} </math>
<math>\Rightarrow</math> Método de Euler
[[Archivo:Ejerciciotreseuler.png|500px|thumb|rigth|Aproximación mediante Euler para S=100]]
{{matlab|codigo=

%Euler
clear all
%Datos
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%a: Es la constante de proporcionalidad entre los infectados por cada
%interacción
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
y0=2000;
S=100;
a=0.003;
b=0.3;
c=0.01;

t(1)=t0;
y(1)=y0;

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de y(i) hasta que y(i)>500
while y(i)>500

y(i+1)=y(i)+h*(a*S*y(i)-b*y(i)-c*y(i));
t(i+1)=t(i)+h;
i=i+1;
end

[t',y'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,y)
legend('Población infectada para S=100(Método de Euler)','Location','best');

}}

<math>\Rightarrow</math> Método del trapecio
[[Archivo:Ejerciciotrestrapecio.png|500px|thumb|rigth|Aproximación mediante el método del trapecio para S=100]]
{{matlab|codigo=

%Trapecio
clear all
%Datos
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%a: Es la constante de proporcionalidad entre los infectados por cada
%interacción
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
z0=2000;
S=100;
a=0.003;
b=0.3;
c=0.01;

t(1)=t0;
z(1)=z0;

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de z(i) hasta que z(i)>500

%TRAPECIO
while z(i)>500
z(i+1)=(z(i)*(1-(h/2)*(b+c-S*a)))/(1+h/2*(b+c-S*a));
t(i+1)=t(i)+h;
i=i+1;
end

[t',z'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,z)
legend('Población infectada para S=100(Método del trapecio)','Location','best');

}}

Ahora con estos nuevos datos vemos que para que se reduzca a la cuarta parte el número de infectados tendrán que pasar 138.7 días, algo lógico al fijarnos en la ED, ambas son ecuaciones decrecientes pero para el caso de S=100 la pendiente será mucho menor.
Si continuamos haciendo pruebas vemos para valores superiores de s, por ejemplo, S=200 el matlab no nos da ninguna respuesta. Esto ocurre porque al aumentar S la ecuación pasa a ser creciente, por lo que lógicamente si intentas buscar un valor del número de infectados menor al inicial no lo podrás encontrar.

==Resolución del problema completo==
Vamos a resolver el problema completo para ver como evolucionarían las poblaciones en un periodo de 40 días.
Primero supongamos una población infectada inicial de 20 individuos y una población en riesgo de contagio de 800.
Después supondremos el caso de 40 individuos enfermos y 10000 en riesgo de contagio.
Para resolverlo usaremos el método de Euler y más tarde el método de Runge-Kutta de orden cuatro con distintas discretizaciones para ver la influencia de estas.

<math>\Rightarrow</math> Método de Euler
{{matlab|codigo=

%Euler
clear all
%Datos
%Introducimos los valores de las constantes
a=0.003;
b=0.3;
c=0.01;
h=input('Introducir valores de h:');
t0=0;
y0=input('Introducir vector [So,Io]:');
tN=40;
% calculamos los subintervalos
t=t0:h:tN;
N=(tN-t0)/h;
y=zeros(2,length(t));
y(:,1)=y0';
for i= 1:N;
y(:,i+1)=y(:,i)+h*[-a*y(1,i).*y(2,i);a*y(1,i).*y(2,i)-b*y(2,i)-c*y(2,i)];
end

% grafico
plot(t,y)
% vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h
legend('Población sana','Location','best','Población enferma','Location','best');

}}

<math>\Rightarrow</math> Método de Runge-Kutta

{{matlab|codigo=

%Runge-Kutta
clear all

%Datos
S0=input('Introducir valor inicial SO: ');
I0=input('Introducir valor inicial IO: ');
a=0.003;
b=0.3;
c=0.01;
%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector solucion
y0=[S0;I0];
y=zeros(2,N+1);
y(:,1)=y0;
for i=1:N
K1S=-a*y(1,i)*y(2,i);
K2S=-a*(y(1,i)+(1/2)*K1S*h)*(y(2,i)+(1/2)*K1S*h);
K3S=-a*(y(1,i)+(1/2)*K2S*h)*(y(2,i)+(1/2)*K2S*h);
K4S=-a*(y(1,i)+K3S*h)*(y(2,i)+K3S*h);
y(1,i+1)=y(1,i)+(h/6)*(K1S+2*K2S+2*K3S+K4S);
K1I=a*y(1,i)*y(2,i)-b*y(2,i)-c*y(2,i);
K2I=a*(y(1,i)+(1/2)*K1I*h)*(y(2,i)+(1/2)*K1I*h)-b*(y(2,i)+(1/2)*K1I*h)-c*(y(2,i)+(1/2)*K1I*h);
K3I=a*(y(1,i)+(1/2)*K2I*h)*(y(2,i)+(1/2)*K2I*h)-b*(y(2,i)+(1/2)*K2I*h)-c*(y(2,i)+(1/2)*K2I*h);
K4I=a*(y(1,i)+K3I*h)*(y(2,i)+K3I*h)-b*(y(2,i)+K3I*h)-c*(y(2,i)+K3I*h);
y(2,i+1)=y(2,i)+(h/6)*(K1I+2*K2I+2*K3I+K4I);
end

%Tabla de resutados
[t',y(1,:)',y(2,:)']
%Gráfico
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'k')
legend('Poblacion sana','Poblacion enferma','Location','best');
hold off

% vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h

}}

La dificultad que hay en el uso del método del trapecio o cualquier otro método implícito está en que las incógnitas de nuestro sistema <math>(S,I)</math> quedan implícitas en ambas ecuaciones y necesitan que se despejen. Al comenzar a realizar cálculos nos damos cuenta que el despeje es imposible, ya que las dos variables incógnita aparecen multiplicándose entre sí, quedando una variable dependiendo de la otra siempre.

===Situación inicial de <math>(S _{0} , I_{0} )=(800,20)</math>===
En estos programas, sustituyendo las <math>h</math> para 0.1, 0.01, 0.001 y 0.0001 obtendremos los siguientes resultados:

-Para <math>h=0.1</math>

[[Archivo:Tborraré.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.1rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

Como se puede ver la población sana se infecta rápidamente y a partir del día 15 prácticamente toda la población ha fallecido.
* '''Método de Euler''': nuestro programa nos indica que en este caso el momento en el que los enfermos son máximos se produce al '''tercer día (2.8)''' y el '''número de enfermos en ese momento es 517'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce al '''tercer día (2.6)''' y el '''número de enfermos en ese momento es 650'''.

-Para <math>h=0.01</math>
[[Archivo:Tborraré2.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.01rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

En el gráfico apenas se encuentra diferencia, lo cual es lógico porque es el mismo problema, pero si sobrepusiéramos las gráficas, lograríamos ver como se va ajustando a la realidad en función del aumento de h.
* '''Método de Euler''': en este caso el momento en el que los enfermos son máximos es también en el '''tercer día (2.7)''' y el número de '''enfermos es 506'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce al '''tercer día (2.69)''' y el número de '''enfermos en ese momento es 517'''.

-Para <math>h=0.001</math>
[[Archivo:Tborraré3.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.001rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': en este caso el momento en el que los enfermos son máximos se produce el '''tercer día (2.696)''' con un valor de '''505 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''tercer día (2.694)''' con un valor de '''506 enfermos'''.

-Para <math>h=0.0001</math>
[[Archivo:Tborraré4.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.0001rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': al igual que antes, se repite el momento en el que los enfermos son máximos en el '''tercer día (2.695)''' con un valor de '''505 enfermos'''.
* '''Método de Runge-Kutta''': se repite el momento en el que los enfermos son máximos en el '''tercer día (2.6948)''' con un valor de '''505 enfermos'''.

Observamos en las distintas gráficas obtenidas, que según vamos empleando un tamaño de paso <math>h</math> más pequeño, nos vamos acercando cada vez más a la población real que pretendemos reflejar.
Nos damos cuenta de que las gráficas tienen bastante lógica, pues a la vez que desciende el número de población sana va aumentando la población infectada.
También observamos que para el método de Runge-Kutta cuando utilizamos tamaño de paso <math>h=0.1</math> pierde precisión en relación a cuando utilizamos tamaños de paso más pequeños.

===Situación inicial de <math>(S _{0} , I_{0} )=(10000,40)</math>===
En estos programas, sustituyendo las <math>h</math> para 0.1, 0.01, 0.001 y 0.0001 obtendremos los siguientes resultados:

-Para <math>h=0.1</math>

[[Archivo:H0.11000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.5)''' con un valor de '''12600 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.4)''' con un valor de '''2.6e74 enfermos'''.

-Para <math>h=0.01</math>
[[Archivo:H0.0110000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.0110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.35)''' con un valor de '''9521 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.31)''' con un valor de '''13550 enfermos'''.

-Para <math>h=0.001</math>
[[Archivo:H0.00110000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.00110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.342)''' con un valor de '''9470 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.338)''' con un valor de '''9758 enfermos'''.

-Para <math>h=0.0001</math>
[[Archivo:H0.000110000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.000110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.3411)''' con un valor de '''9465 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.3407)''' con un valor de '''9493 enfermos'''.

Para este segundo caso, nos damos cuenta que al utilizar tamaño de paso <math>h=0.1</math> , llegamos a un resultado sin sentido, pues nos da valores negativos y es imposible que una población tenga tales valores.
En las demás gráficas, los resultados obtenidos no parecen cuadrar con la realidad. Esto podría ser debido a que la alta diferencia entre los valores iniciales <math>(S _{0} , I_{0} )=(10000,40)</math> no permiten un buen cálculo de su variación respecto del tiempo.
El número inicial de infectados es muy elevado y al disminuir el tamaño de paso, el programa aproxima mejor la solución.

==Resolución del problema completo suponiendo <math> a </math> una función variable con el tiempo==
Para el siguiente caso supondremos que el coeficiente <math> a </math> (tasa de contagio en población sana) es una función dada por :

<math> a(t)= \frac{0.003}{(1+t)} </math>
Suponiendo además las condiciones iniciales <math> (S_{0},I_{0})=(1640,40) </math> , vamos a dibujar las gráficas de infectados y sanos por el método de Heun, tomando un tamaño de paso <math> h=0.01 </math>
[[Archivo: H0.01heunn.png|1000px|thumb|rigth| Población infectada-Población sana para a=a(t)]]
{{matlab|codigo=
%Heun. Coeficiente variable a(t)

%Datos
S0=input('Introducir valor inicial SO: ');
I0=input('Introducir valor inicial IO: ');
b=0.3;
c=0.01;
%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector solucion
y0=[S0;I0];
y=zeros(2,N+1);
y(:,1)=y0;
for i=1:N
K1S=-(0.003/(1+t(i)))*y(1,i)*y(2,i);
K2S=-(0.003/(1+t(i)+h))*(y(1,i)+K1S*h)*(y(2,i)+K1S*h);
y(1,i+1)=y(1,i)+(h/2)*(K1S+K2S);
K1I=(0.003/(1+t(i)))*y(1,i)*y(2,i)-b*y(2,i)-c*y(2,i);
K2I=(0.003/(1+t(i)+h))*(y(1,i)+K1I*h)*(y(2,i)+K1I*h)-b*(y(2,i)+K1I*h)-c*(y(2,i)+K1I*h);
y(2,i+1)=y(2,i)+(h/2)*(K1I+K2I);
end

%Tabla de resutados
[t',y(1,:)',y(2,:)']
%Graficos
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'k')
title('Método de Heun');
legend('Poblacion sana','Población infectada','Location','best');
hold off

% Vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h

}}

Como se puede ver la población sana se infecta rápidamente y a partir del día 15 casi toda la población ha fallecido. Al segundo día <math>(t=2,15)</math> observamos que el número de infectados es de 994 individuos <math>(I=994,3)</math>, siendo este el máximo que se alcanzará.

A medida que va pasando el tiempo el coeficiente <math> a(t) </math> (tasa de contagio en población sana) va a ir decreciendo poco a poco hasta aproximarse a cero según podemos ver en el siguiente gráfico:
<math> \lim_{t\to \infty } a(t)=0</math>
[[Archivo: A(t).png|450px|thumb|right|a(t) para t=[0,40]]]
{{matlab|codigo=

%tasa de contagio en poblacion sana a(t)

%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector a
a=zeros(1,N+1);
for i=1:N+1
a(i)=0.003/(1+t(i));
end

%Grafica
plot(t,a)

}}

[[Archivo: Comparacionheun.png|border|1500px|center|frame|Poblaciones según la tasa de contagio a(t)]]

Ambas poblaciones parten del mismo inicio (para <math>t=0</math> las poblaciones son iguales). Podemos observar que la población sana (susceptible de contraer la enfermedad) en el caso en el que la tasa de contagio en población sana <math>(a)</math> es constante <math>(a=0.003)</math> se ve reducida antes que en el caso en el que depende del tiempo <math>(a=a(t))</math>.

<math>a(t)</math> disminuye a la vez que transcurre el tiempo, por lo que en general las pendientes tanto de la población sana como de la infectada serán menores para este caso que para el caso en el que <math>a=0.003</math>

==Calibración del parámetro <math>a</math> en base a una experiencia previa==
Ahora vamos a intentar calibrar el párametro <math>a</math> para que el número de infectados tenga el máximo en el quinto día, para ello escribimos el siguiente código:

{{matlab|codigo=
%Aplicamos el Método de Heunn

%Datos
t0=0;
tN=40;
h=0.1;
t=t0:h:tN;
a=0.0005:0.0001:0.002;
b=0.3;
c=0.01;
N1=round((tN-t0)/h);
N2=length(a);
Tmax=zeros(1,N2);

%Vectores iniciales de Población e Infectados
y=zeros(2,N1+1);

%Condiciones iniciales
y(1,1)=1600;
y(2,1)=40;

%Bucle que nos calcula todas las posibles soluciones
for j=1:N2
for i=1:N1
K1S=-a(j)*y(1,i)*y(2,i);
K2S=-a(j)*y(1,i)*y(2,i)+K1S*h;
y(1,i+1)=y(1,i)+h/2*(K1S+K2S);

K1I=y(2,i)*(a(j)*y(1,i)-(b+c));
K2I=y(2,i)*(a(j)*y(1,i)-(b+c))+K1I*h;
y(2,i+1)=y(2,i)+h/2*(K1I+K2I);

if i>1 && y(2,i+1)<y(2,i) && y(2,i)>y(2,i-1)
Tmax(j)=t(i);
end
end
end
distancia=(abs(Tmax-5));
[valormin,posmin]=min(distancia);
a=a(posmin)
}}

Con este código conseguimos

Modelos epidemiológicos (Grupo 3A)

2016-05-02T18:32:35Z

Trabajocampos: /* Calibración del parámetro a en base a una experiencia previa */

{{ TrabajoED | Modelos Epidemológicos (Grupo 3A) | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED15/16|Curso 2015-16]] |
Ignacio Mollá Carcaño 
Pablo Revuelta Aragón 
David González Hernández 
Jose María García Rodríguez 
Alejandro Martínez Gamonal }}

En el desarrollo de una epidemia, se distinguen dos tipos de individuos. Los que ya han contraído la enfermedad (infectados), que llamaremos <math>I(t)</math> ; y los que son susceptibles de contraerla (sanos), a los que llamaremos <math>S(t)</math>. Donde <math>t</math> es la variable del tiempo. Se presentan dos hipótesis para realizar este estudio:

1. La población de personas infectadas se altera por el fallecimiento o la cura de las mismas. En ambos casos, la tasa de cambio depende del número 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ón entre el número de individuos en ambas clases.

A través del siguiente sistema de ecuaciones diferenciales en función del tiempo observamos las variaciones de ambas poblaciones:

<math>\begin{cases}\frac{dS}{dt}=-aSI\\\frac{dI}{dt}=aSI-bI-cI\end{cases}</math>

Donde <math>a</math>, <math>b</math> y <math>c</math> son parámetros.

* <math>\frac{dS}{dt}</math> : Es la variación de la población susceptible de contraer la enfermedad.

* <math>\frac{dI}{dt}</math> : Es la variación de la población infectada.

* <math>aSI</math> : Es la interacción entre ambas poblaciones que dependen de <math>a</math>.

* <math>bI</math> : Indica los individuos que fallecen. Por lo tanto <math>b</math> podría considerarse como la '''tasa de mortalidad'''.

* <math>cI</math> : Indica los individuos que se han curado. Por lo tanto <math>c</math> podría considerarse como la '''tasa de curación'''.

Para interpretar <math>a</math> observamos el sistema y la segunda hipótesis. Como los individuos que pasan de ser susceptibles a estar infectados son proporcionales a su interacción <math>(SI)</math>, <math>a</math> podría considerarse la '''tasa de contagio''' en la población susceptible.

==Resolución del problema con una sola ecuación diferencial==
Primero vamamos a estudiar el modelo epidemiológico de una población en la que incialmente las 2000 personas que conforman la población se encuentran infectadas. Para hacer el seguimiento de la enfermedad contamos tanto con la constante de proporcionalidad de personas que se curan (b=0.3) y la constante de proporcionalidad del número de personas que fallecen (c=0.01). De esta forma nos quedará el siguiente PVI:

<math> \begin{cases} \frac{dI}{dt} = -0.31 I \\I_{0} = 2000\end{cases} </math>

A continuación vamos a proceder a resolverlo por el método de Euler y trapecio asignandole un tamaño de paso de h=0.1 y veremos cuanto tiempo tardará en reducirse el número de infectados a la cuarta parte.

<math>\Rightarrow</math> Método de Euler
[[Archivo:Ejerciciodoseuler.png|500px|thumb|rigth|Aproximación mediante Euler]]
{{matlab|codigo=

clear all
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
y0=2000;
b=0.3;
c=0.01;

t(1)=t0;
ye(1)=y0; %Vector solución Euler

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de y(i) hasta que y(i)>500

%EULER

while ye(i)>500

ye(i+1)=ye(i)+h*(-b*ye(i)-c*ye(i));
t(i+1)=t(i)+h;
i=i+1;
end

[t',ye'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,ye)
legend('Población infectada(Método de Euler)','Location','best');

}}

<math>\Rightarrow</math> Método del trapecio
[[Archivo:Ejerciciodostrapecio.png|500px|thumb|rigth|Aproximación mediante el método del trapecio]]
{{matlab|codigo=
clear all
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
z0=2000;
b=0.3;
c=0.01;

t(1)=t0;
z(1)=z0; %Vector solución trapecio

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de z(i) hasta que z(i)>500

%TRAPECIO
while z(i)>500
z(i+1)=(z(i)*(1-(h/2)*(b+c)))/(1+h/2*(b+c));%trapecio
t(i+1)=t(i)+h;
i=i+1;
end

[t',z'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,z)
legend('Población infectada(Método del trapecio)','Location','best');

}}

Con todo esto podemos determinar que el número de infectados se reducirá a la cuarta parte en 4.5 días, esto ocurrirá o bien por el fallecimiento o por la curación de dicha enfermedad. Cabe señalar que el número de curados será mayor que el de defunciones, esto se puede ver a simple vista viendo que b>a.
De cara a nuestro estudio de esta enfermedad vamos a ver como afectaría al desarrollo de la misma si introducimos 100 sujetos sanos en la población, tomando una constante de interacción entre personas infectadas y sanas de a=0.003. Procederemos igual que antes a resolverlo por el método de Euler y trapecio con un tamaño de paso de 0.1.
Primero vamos a definir nuestro nuevo PVI:
<math> \begin{cases} \frac{dI}{dt} = -0.01 I \\I_{0} = 2000\end{cases} </math>
<math>\Rightarrow</math> Método de Euler
[[Archivo:Ejerciciotreseuler.png|500px|thumb|rigth|Aproximación mediante Euler para S=100]]
{{matlab|codigo=

%Euler
clear all
%Datos
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%a: Es la constante de proporcionalidad entre los infectados por cada
%interacción
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
y0=2000;
S=100;
a=0.003;
b=0.3;
c=0.01;

t(1)=t0;
y(1)=y0;

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de y(i) hasta que y(i)>500
while y(i)>500

y(i+1)=y(i)+h*(a*S*y(i)-b*y(i)-c*y(i));
t(i+1)=t(i)+h;
i=i+1;
end

[t',y'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,y)
legend('Población infectada para S=100(Método de Euler)','Location','best');

}}

<math>\Rightarrow</math> Método del trapecio
[[Archivo:Ejerciciotrestrapecio.png|500px|thumb|rigth|Aproximación mediante el método del trapecio para S=100]]
{{matlab|codigo=

%Trapecio
clear all
%Datos
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%a: Es la constante de proporcionalidad entre los infectados por cada
%interacción
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
z0=2000;
S=100;
a=0.003;
b=0.3;
c=0.01;

t(1)=t0;
z(1)=z0;

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de z(i) hasta que z(i)>500

%TRAPECIO
while z(i)>500
z(i+1)=(z(i)*(1-(h/2)*(b+c-S*a)))/(1+h/2*(b+c-S*a));
t(i+1)=t(i)+h;
i=i+1;
end

[t',z'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,z)
legend('Población infectada para S=100(Método del trapecio)','Location','best');

}}

Ahora con estos nuevos datos vemos que para que se reduzca a la cuarta parte el número de infectados tendrán que pasar 138.7 días, algo lógico al fijarnos en la ED, ambas son ecuaciones decrecientes pero para el caso de S=100 la pendiente será mucho menor.
Si continuamos haciendo pruebas vemos para valores superiores de s, por ejemplo, S=200 el matlab no nos da ninguna respuesta. Esto ocurre porque al aumentar S la ecuación pasa a ser creciente, por lo que lógicamente si intentas buscar un valor del número de infectados menor al inicial no lo podrás encontrar.

==Resolución del problema completo==
Vamos a resolver el problema completo para ver como evolucionarían las poblaciones en un periodo de 40 días.
Primero supongamos una población infectada inicial de 20 individuos y una población en riesgo de contagio de 800.
Después supondremos el caso de 40 individuos enfermos y 10000 en riesgo de contagio.
Para resolverlo usaremos el método de Euler y más tarde el método de Runge-Kutta de orden cuatro con distintas discretizaciones para ver la influencia de estas.

<math>\Rightarrow</math> Método de Euler
{{matlab|codigo=

%Euler
clear all
%Datos
%Introducimos los valores de las constantes
a=0.003;
b=0.3;
c=0.01;
h=input('Introducir valores de h:');
t0=0;
y0=input('Introducir vector [So,Io]:');
tN=40;
% calculamos los subintervalos
t=t0:h:tN;
N=(tN-t0)/h;
y=zeros(2,length(t));
y(:,1)=y0';
for i= 1:N;
y(:,i+1)=y(:,i)+h*[-a*y(1,i).*y(2,i);a*y(1,i).*y(2,i)-b*y(2,i)-c*y(2,i)];
end

% grafico
plot(t,y)
% vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h
legend('Población sana','Location','best','Población enferma','Location','best');

}}

<math>\Rightarrow</math> Método de Runge-Kutta

{{matlab|codigo=

%Runge-Kutta
clear all

%Datos
S0=input('Introducir valor inicial SO: ');
I0=input('Introducir valor inicial IO: ');
a=0.003;
b=0.3;
c=0.01;
%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector solucion
y0=[S0;I0];
y=zeros(2,N+1);
y(:,1)=y0;
for i=1:N
K1S=-a*y(1,i)*y(2,i);
K2S=-a*(y(1,i)+(1/2)*K1S*h)*(y(2,i)+(1/2)*K1S*h);
K3S=-a*(y(1,i)+(1/2)*K2S*h)*(y(2,i)+(1/2)*K2S*h);
K4S=-a*(y(1,i)+K3S*h)*(y(2,i)+K3S*h);
y(1,i+1)=y(1,i)+(h/6)*(K1S+2*K2S+2*K3S+K4S);
K1I=a*y(1,i)*y(2,i)-b*y(2,i)-c*y(2,i);
K2I=a*(y(1,i)+(1/2)*K1I*h)*(y(2,i)+(1/2)*K1I*h)-b*(y(2,i)+(1/2)*K1I*h)-c*(y(2,i)+(1/2)*K1I*h);
K3I=a*(y(1,i)+(1/2)*K2I*h)*(y(2,i)+(1/2)*K2I*h)-b*(y(2,i)+(1/2)*K2I*h)-c*(y(2,i)+(1/2)*K2I*h);
K4I=a*(y(1,i)+K3I*h)*(y(2,i)+K3I*h)-b*(y(2,i)+K3I*h)-c*(y(2,i)+K3I*h);
y(2,i+1)=y(2,i)+(h/6)*(K1I+2*K2I+2*K3I+K4I);
end

%Tabla de resutados
[t',y(1,:)',y(2,:)']
%Gráfico
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'k')
legend('Poblacion sana','Poblacion enferma','Location','best');
hold off

% vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h

}}

La dificultad que hay en el uso del método del trapecio o cualquier otro método implícito está en que las incógnitas de nuestro sistema <math>(S,I)</math> quedan implícitas en ambas ecuaciones y necesitan que se despejen. Al comenzar a realizar cálculos nos damos cuenta que el despeje es imposible, ya que las dos variables incógnita aparecen multiplicándose entre sí, quedando una variable dependiendo de la otra siempre.

===Situación inicial de <math>(S _{0} , I_{0} )=(800,20)</math>===
En estos programas, sustituyendo las <math>h</math> para 0.1, 0.01, 0.001 y 0.0001 obtendremos los siguientes resultados:

-Para <math>h=0.1</math>

[[Archivo:Tborraré.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.1rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

Como se puede ver la población sana se infecta rápidamente y a partir del día 15 prácticamente toda la población ha fallecido.
* '''Método de Euler''': nuestro programa nos indica que en este caso el momento en el que los enfermos son máximos se produce al '''tercer día (2.8)''' y el '''número de enfermos en ese momento es 517'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce al '''tercer día (2.6)''' y el '''número de enfermos en ese momento es 650'''.

-Para <math>h=0.01</math>
[[Archivo:Tborraré2.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.01rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

En el gráfico apenas se encuentra diferencia, lo cual es lógico porque es el mismo problema, pero si sobrepusiéramos las gráficas, lograríamos ver como se va ajustando a la realidad en función del aumento de h.
* '''Método de Euler''': en este caso el momento en el que los enfermos son máximos es también en el '''tercer día (2.7)''' y el número de '''enfermos es 506'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce al '''tercer día (2.69)''' y el número de '''enfermos en ese momento es 517'''.

-Para <math>h=0.001</math>
[[Archivo:Tborraré3.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.001rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': en este caso el momento en el que los enfermos son máximos se produce el '''tercer día (2.696)''' con un valor de '''505 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''tercer día (2.694)''' con un valor de '''506 enfermos'''.

-Para <math>h=0.0001</math>
[[Archivo:Tborraré4.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.0001rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': al igual que antes, se repite el momento en el que los enfermos son máximos en el '''tercer día (2.695)''' con un valor de '''505 enfermos'''.
* '''Método de Runge-Kutta''': se repite el momento en el que los enfermos son máximos en el '''tercer día (2.6948)''' con un valor de '''505 enfermos'''.

Observamos en las distintas gráficas obtenidas, que según vamos empleando un tamaño de paso <math>h</math> más pequeño, nos vamos acercando cada vez más a la población real que pretendemos reflejar.
Nos damos cuenta de que las gráficas tienen bastante lógica, pues a la vez que desciende el número de población sana va aumentando la población infectada.
También observamos que para el método de Runge-Kutta cuando utilizamos tamaño de paso <math>h=0.1</math> pierde precisión en relación a cuando utilizamos tamaños de paso más pequeños.

===Situación inicial de <math>(S _{0} , I_{0} )=(10000,40)</math>===
En estos programas, sustituyendo las <math>h</math> para 0.1, 0.01, 0.001 y 0.0001 obtendremos los siguientes resultados:

-Para <math>h=0.1</math>

[[Archivo:H0.11000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.5)''' con un valor de '''12600 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.4)''' con un valor de '''2.6e74 enfermos'''.

-Para <math>h=0.01</math>
[[Archivo:H0.0110000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.0110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.35)''' con un valor de '''9521 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.31)''' con un valor de '''13550 enfermos'''.

-Para <math>h=0.001</math>
[[Archivo:H0.00110000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.00110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.342)''' con un valor de '''9470 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.338)''' con un valor de '''9758 enfermos'''.

-Para <math>h=0.0001</math>
[[Archivo:H0.000110000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.000110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.3411)''' con un valor de '''9465 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.3407)''' con un valor de '''9493 enfermos'''.

Para este segundo caso, nos damos cuenta que al utilizar tamaño de paso <math>h=0.1</math> , llegamos a un resultado sin sentido, pues nos da valores negativos y es imposible que una población tenga tales valores.
En las demás gráficas, los resultados obtenidos no parecen cuadrar con la realidad. Esto podría ser debido a que la alta diferencia entre los valores iniciales <math>(S _{0} , I_{0} )=(10000,40)</math> no permiten un buen cálculo de su variación respecto del tiempo.
El número inicial de infectados es muy elevado y al disminuir el tamaño de paso, el programa aproxima mejor la solución.

==Resolución del problema completo suponiendo <math> a </math> una función variable con el tiempo==
Para el siguiente caso supondremos que el coeficiente <math> a </math> (tasa de contagio en población sana) es una función dada por :

<math> a(t)= \frac{0.003}{(1+t)} </math>
Suponiendo además las condiciones iniciales <math> (S_{0},I_{0})=(1640,40) </math> , vamos a dibujar las gráficas de infectados y sanos por el método de Heun, tomando un tamaño de paso <math> h=0.01 </math>
[[Archivo: H0.01heunn.png|1000px|thumb|rigth| Población infectada-Población sana para a=a(t)]]
{{matlab|codigo=
%Heun. Coeficiente variable a(t)

%Datos
S0=input('Introducir valor inicial SO: ');
I0=input('Introducir valor inicial IO: ');
b=0.3;
c=0.01;
%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector solucion
y0=[S0;I0];
y=zeros(2,N+1);
y(:,1)=y0;
for i=1:N
K1S=-(0.003/(1+t(i)))*y(1,i)*y(2,i);
K2S=-(0.003/(1+t(i)+h))*(y(1,i)+K1S*h)*(y(2,i)+K1S*h);
y(1,i+1)=y(1,i)+(h/2)*(K1S+K2S);
K1I=(0.003/(1+t(i)))*y(1,i)*y(2,i)-b*y(2,i)-c*y(2,i);
K2I=(0.003/(1+t(i)+h))*(y(1,i)+K1I*h)*(y(2,i)+K1I*h)-b*(y(2,i)+K1I*h)-c*(y(2,i)+K1I*h);
y(2,i+1)=y(2,i)+(h/2)*(K1I+K2I);
end

%Tabla de resutados
[t',y(1,:)',y(2,:)']
%Graficos
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'k')
title('Método de Heun');
legend('Poblacion sana','Población infectada','Location','best');
hold off

% Vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h

}}

Como se puede ver la población sana se infecta rápidamente y a partir del día 15 casi toda la población ha fallecido. Al segundo día <math>(t=2,15)</math> observamos que el número de infectados es de 994 individuos <math>(I=994,3)</math>, siendo este el máximo que se alcanzará.

A medida que va pasando el tiempo el coeficiente <math> a(t) </math> (tasa de contagio en población sana) va a ir decreciendo poco a poco hasta aproximarse a cero según podemos ver en el siguiente gráfico:
<math> \lim_{t\to \infty } a(t)=0</math>
[[Archivo: A(t).png|450px|thumb|right|a(t) para t=[0,40]]]
{{matlab|codigo=

%tasa de contagio en poblacion sana a(t)

%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector a
a=zeros(1,N+1);
for i=1:N+1
a(i)=0.003/(1+t(i));
end

%Grafica
plot(t,a)

}}

[[Archivo: Comparacionheun.png|border|1500px|center|frame|Poblaciones según la tasa de contagio a(t)]]

Ambas poblaciones parten del mismo inicio (para <math>t=0</math> las poblaciones son iguales). Podemos observar que la población sana (susceptible de contraer la enfermedad) en el caso en el que la tasa de contagio en población sana <math>(a)</math> es constante <math>(a=0.003)</math> se ve reducida antes que en el caso en el que depende del tiempo <math>(a=a(t))</math>.

<math>a(t)</math> disminuye a la vez que transcurre el tiempo, por lo que en general las pendientes tanto de la población sana como de la infectada serán menores para este caso que para el caso en el que <math>a=0.003</math>

==Calibración del parámetro <math>a</math> en base a una experiencia previa==
Ahora vamos a intentar calibrar el párametro <math>a</math> para que el número de infectados tenga el máximo en el quinto día, para ello escribimos el siguiente código:

{{matlab|codigo=
%Aplicamos el Método de Heunn

%Datos
t0=0;
tN=40;
h=0.1;
t=t0:h:tN;
a=0.0005:0.0001:0.002;
b=0.3;
c=0.01;
N1=round((tN-t0)/h);
N2=length(a);
Tmax=zeros(1,N2);

%Vectores iniciales de Población e Infectados
S=zeros(1,N1+1);
I=zeros(1,N1+1);

%Condiciones iniciales
S(1)=1600;
I(1)=40;

%Bucle que nos calcula todas las posibles soluciones
for j=1:N2
for i=1:N1
K1S=-a(j)*S(i)*I(i);
K2S=-a(j)*S(i)*I(i)+K1S*h;
S(i+1)=S(i)+h/2*(K1S+K2S);

K1I=I(i)*(a(j)*S(i)-(b+c));
K2I=I(i)*(a(j)*S(i)-(b+c))+K1I*h;
I(i+1)=I(i)+h/2*(K1I+K2I);

if i>1 && I(i+1)<I(i) && I(i)>I(i-1)
Tmax(j)=t(i);
end
end
end
distancia=(abs(Tmax-5));
[valormin,posmin]=min(distancia);
a=a(posmin)
}}

Modelos epidemiológicos (Grupo 3A)

2016-05-02T18:30:50Z

Trabajocampos: /* Calibración del parámetro a en base a una experiencia previa */

{{ TrabajoED | Modelos Epidemológicos (Grupo 3A) | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED15/16|Curso 2015-16]] |
Ignacio Mollá Carcaño 
Pablo Revuelta Aragón 
David González Hernández 
Jose María García Rodríguez 
Alejandro Martínez Gamonal }}

En el desarrollo de una epidemia, se distinguen dos tipos de individuos. Los que ya han contraído la enfermedad (infectados), que llamaremos <math>I(t)</math> ; y los que son susceptibles de contraerla (sanos), a los que llamaremos <math>S(t)</math>. Donde <math>t</math> es la variable del tiempo. Se presentan dos hipótesis para realizar este estudio:

1. La población de personas infectadas se altera por el fallecimiento o la cura de las mismas. En ambos casos, la tasa de cambio depende del número 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ón entre el número de individuos en ambas clases.

A través del siguiente sistema de ecuaciones diferenciales en función del tiempo observamos las variaciones de ambas poblaciones:

<math>\begin{cases}\frac{dS}{dt}=-aSI\\\frac{dI}{dt}=aSI-bI-cI\end{cases}</math>

Donde <math>a</math>, <math>b</math> y <math>c</math> son parámetros.

* <math>\frac{dS}{dt}</math> : Es la variación de la población susceptible de contraer la enfermedad.

* <math>\frac{dI}{dt}</math> : Es la variación de la población infectada.

* <math>aSI</math> : Es la interacción entre ambas poblaciones que dependen de <math>a</math>.

* <math>bI</math> : Indica los individuos que fallecen. Por lo tanto <math>b</math> podría considerarse como la '''tasa de mortalidad'''.

* <math>cI</math> : Indica los individuos que se han curado. Por lo tanto <math>c</math> podría considerarse como la '''tasa de curación'''.

Para interpretar <math>a</math> observamos el sistema y la segunda hipótesis. Como los individuos que pasan de ser susceptibles a estar infectados son proporcionales a su interacción <math>(SI)</math>, <math>a</math> podría considerarse la '''tasa de contagio''' en la población susceptible.

==Resolución del problema con una sola ecuación diferencial==
Primero vamamos a estudiar el modelo epidemiológico de una población en la que incialmente las 2000 personas que conforman la población se encuentran infectadas. Para hacer el seguimiento de la enfermedad contamos tanto con la constante de proporcionalidad de personas que se curan (b=0.3) y la constante de proporcionalidad del número de personas que fallecen (c=0.01). De esta forma nos quedará el siguiente PVI:

<math> \begin{cases} \frac{dI}{dt} = -0.31 I \\I_{0} = 2000\end{cases} </math>

A continuación vamos a proceder a resolverlo por el método de Euler y trapecio asignandole un tamaño de paso de h=0.1 y veremos cuanto tiempo tardará en reducirse el número de infectados a la cuarta parte.

<math>\Rightarrow</math> Método de Euler
[[Archivo:Ejerciciodoseuler.png|500px|thumb|rigth|Aproximación mediante Euler]]
{{matlab|codigo=

clear all
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
y0=2000;
b=0.3;
c=0.01;

t(1)=t0;
ye(1)=y0; %Vector solución Euler

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de y(i) hasta que y(i)>500

%EULER

while ye(i)>500

ye(i+1)=ye(i)+h*(-b*ye(i)-c*ye(i));
t(i+1)=t(i)+h;
i=i+1;
end

[t',ye'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,ye)
legend('Población infectada(Método de Euler)','Location','best');

}}

<math>\Rightarrow</math> Método del trapecio
[[Archivo:Ejerciciodostrapecio.png|500px|thumb|rigth|Aproximación mediante el método del trapecio]]
{{matlab|codigo=
clear all
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
z0=2000;
b=0.3;
c=0.01;

t(1)=t0;
z(1)=z0; %Vector solución trapecio

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de z(i) hasta que z(i)>500

%TRAPECIO
while z(i)>500
z(i+1)=(z(i)*(1-(h/2)*(b+c)))/(1+h/2*(b+c));%trapecio
t(i+1)=t(i)+h;
i=i+1;
end

[t',z'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,z)
legend('Población infectada(Método del trapecio)','Location','best');

}}

Con todo esto podemos determinar que el número de infectados se reducirá a la cuarta parte en 4.5 días, esto ocurrirá o bien por el fallecimiento o por la curación de dicha enfermedad. Cabe señalar que el número de curados será mayor que el de defunciones, esto se puede ver a simple vista viendo que b>a.
De cara a nuestro estudio de esta enfermedad vamos a ver como afectaría al desarrollo de la misma si introducimos 100 sujetos sanos en la población, tomando una constante de interacción entre personas infectadas y sanas de a=0.003. Procederemos igual que antes a resolverlo por el método de Euler y trapecio con un tamaño de paso de 0.1.
Primero vamos a definir nuestro nuevo PVI:
<math> \begin{cases} \frac{dI}{dt} = -0.01 I \\I_{0} = 2000\end{cases} </math>
<math>\Rightarrow</math> Método de Euler
[[Archivo:Ejerciciotreseuler.png|500px|thumb|rigth|Aproximación mediante Euler para S=100]]
{{matlab|codigo=

%Euler
clear all
%Datos
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%a: Es la constante de proporcionalidad entre los infectados por cada
%interacción
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
y0=2000;
S=100;
a=0.003;
b=0.3;
c=0.01;

t(1)=t0;
y(1)=y0;

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de y(i) hasta que y(i)>500
while y(i)>500

y(i+1)=y(i)+h*(a*S*y(i)-b*y(i)-c*y(i));
t(i+1)=t(i)+h;
i=i+1;
end

[t',y'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,y)
legend('Población infectada para S=100(Método de Euler)','Location','best');

}}

<math>\Rightarrow</math> Método del trapecio
[[Archivo:Ejerciciotrestrapecio.png|500px|thumb|rigth|Aproximación mediante el método del trapecio para S=100]]
{{matlab|codigo=

%Trapecio
clear all
%Datos
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%a: Es la constante de proporcionalidad entre los infectados por cada
%interacción
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
z0=2000;
S=100;
a=0.003;
b=0.3;
c=0.01;

t(1)=t0;
z(1)=z0;

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de z(i) hasta que z(i)>500

%TRAPECIO
while z(i)>500
z(i+1)=(z(i)*(1-(h/2)*(b+c-S*a)))/(1+h/2*(b+c-S*a));
t(i+1)=t(i)+h;
i=i+1;
end

[t',z'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,z)
legend('Población infectada para S=100(Método del trapecio)','Location','best');

}}

Ahora con estos nuevos datos vemos que para que se reduzca a la cuarta parte el número de infectados tendrán que pasar 138.7 días, algo lógico al fijarnos en la ED, ambas son ecuaciones decrecientes pero para el caso de S=100 la pendiente será mucho menor.
Si continuamos haciendo pruebas vemos para valores superiores de s, por ejemplo, S=200 el matlab no nos da ninguna respuesta. Esto ocurre porque al aumentar S la ecuación pasa a ser creciente, por lo que lógicamente si intentas buscar un valor del número de infectados menor al inicial no lo podrás encontrar.

==Resolución del problema completo==
Vamos a resolver el problema completo para ver como evolucionarían las poblaciones en un periodo de 40 días.
Primero supongamos una población infectada inicial de 20 individuos y una población en riesgo de contagio de 800.
Después supondremos el caso de 40 individuos enfermos y 10000 en riesgo de contagio.
Para resolverlo usaremos el método de Euler y más tarde el método de Runge-Kutta de orden cuatro con distintas discretizaciones para ver la influencia de estas.

<math>\Rightarrow</math> Método de Euler
{{matlab|codigo=

%Euler
clear all
%Datos
%Introducimos los valores de las constantes
a=0.003;
b=0.3;
c=0.01;
h=input('Introducir valores de h:');
t0=0;
y0=input('Introducir vector [So,Io]:');
tN=40;
% calculamos los subintervalos
t=t0:h:tN;
N=(tN-t0)/h;
y=zeros(2,length(t));
y(:,1)=y0';
for i= 1:N;
y(:,i+1)=y(:,i)+h*[-a*y(1,i).*y(2,i);a*y(1,i).*y(2,i)-b*y(2,i)-c*y(2,i)];
end

% grafico
plot(t,y)
% vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h
legend('Población sana','Location','best','Población enferma','Location','best');

}}

<math>\Rightarrow</math> Método de Runge-Kutta

{{matlab|codigo=

%Runge-Kutta
clear all

%Datos
S0=input('Introducir valor inicial SO: ');
I0=input('Introducir valor inicial IO: ');
a=0.003;
b=0.3;
c=0.01;
%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector solucion
y0=[S0;I0];
y=zeros(2,N+1);
y(:,1)=y0;
for i=1:N
K1S=-a*y(1,i)*y(2,i);
K2S=-a*(y(1,i)+(1/2)*K1S*h)*(y(2,i)+(1/2)*K1S*h);
K3S=-a*(y(1,i)+(1/2)*K2S*h)*(y(2,i)+(1/2)*K2S*h);
K4S=-a*(y(1,i)+K3S*h)*(y(2,i)+K3S*h);
y(1,i+1)=y(1,i)+(h/6)*(K1S+2*K2S+2*K3S+K4S);
K1I=a*y(1,i)*y(2,i)-b*y(2,i)-c*y(2,i);
K2I=a*(y(1,i)+(1/2)*K1I*h)*(y(2,i)+(1/2)*K1I*h)-b*(y(2,i)+(1/2)*K1I*h)-c*(y(2,i)+(1/2)*K1I*h);
K3I=a*(y(1,i)+(1/2)*K2I*h)*(y(2,i)+(1/2)*K2I*h)-b*(y(2,i)+(1/2)*K2I*h)-c*(y(2,i)+(1/2)*K2I*h);
K4I=a*(y(1,i)+K3I*h)*(y(2,i)+K3I*h)-b*(y(2,i)+K3I*h)-c*(y(2,i)+K3I*h);
y(2,i+1)=y(2,i)+(h/6)*(K1I+2*K2I+2*K3I+K4I);
end

%Tabla de resutados
[t',y(1,:)',y(2,:)']
%Gráfico
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'k')
legend('Poblacion sana','Poblacion enferma','Location','best');
hold off

% vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h

}}

La dificultad que hay en el uso del método del trapecio o cualquier otro método implícito está en que las incógnitas de nuestro sistema <math>(S,I)</math> quedan implícitas en ambas ecuaciones y necesitan que se despejen. Al comenzar a realizar cálculos nos damos cuenta que el despeje es imposible, ya que las dos variables incógnita aparecen multiplicándose entre sí, quedando una variable dependiendo de la otra siempre.

===Situación inicial de <math>(S _{0} , I_{0} )=(800,20)</math>===
En estos programas, sustituyendo las <math>h</math> para 0.1, 0.01, 0.001 y 0.0001 obtendremos los siguientes resultados:

-Para <math>h=0.1</math>

[[Archivo:Tborraré.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.1rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

Como se puede ver la población sana se infecta rápidamente y a partir del día 15 prácticamente toda la población ha fallecido.
* '''Método de Euler''': nuestro programa nos indica que en este caso el momento en el que los enfermos son máximos se produce al '''tercer día (2.8)''' y el '''número de enfermos en ese momento es 517'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce al '''tercer día (2.6)''' y el '''número de enfermos en ese momento es 650'''.

-Para <math>h=0.01</math>
[[Archivo:Tborraré2.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.01rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

En el gráfico apenas se encuentra diferencia, lo cual es lógico porque es el mismo problema, pero si sobrepusiéramos las gráficas, lograríamos ver como se va ajustando a la realidad en función del aumento de h.
* '''Método de Euler''': en este caso el momento en el que los enfermos son máximos es también en el '''tercer día (2.7)''' y el número de '''enfermos es 506'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce al '''tercer día (2.69)''' y el número de '''enfermos en ese momento es 517'''.

-Para <math>h=0.001</math>
[[Archivo:Tborraré3.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.001rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': en este caso el momento en el que los enfermos son máximos se produce el '''tercer día (2.696)''' con un valor de '''505 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''tercer día (2.694)''' con un valor de '''506 enfermos'''.

-Para <math>h=0.0001</math>
[[Archivo:Tborraré4.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.0001rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': al igual que antes, se repite el momento en el que los enfermos son máximos en el '''tercer día (2.695)''' con un valor de '''505 enfermos'''.
* '''Método de Runge-Kutta''': se repite el momento en el que los enfermos son máximos en el '''tercer día (2.6948)''' con un valor de '''505 enfermos'''.

Observamos en las distintas gráficas obtenidas, que según vamos empleando un tamaño de paso <math>h</math> más pequeño, nos vamos acercando cada vez más a la población real que pretendemos reflejar.
Nos damos cuenta de que las gráficas tienen bastante lógica, pues a la vez que desciende el número de población sana va aumentando la población infectada.
También observamos que para el método de Runge-Kutta cuando utilizamos tamaño de paso <math>h=0.1</math> pierde precisión en relación a cuando utilizamos tamaños de paso más pequeños.

===Situación inicial de <math>(S _{0} , I_{0} )=(10000,40)</math>===
En estos programas, sustituyendo las <math>h</math> para 0.1, 0.01, 0.001 y 0.0001 obtendremos los siguientes resultados:

-Para <math>h=0.1</math>

[[Archivo:H0.11000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.5)''' con un valor de '''12600 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.4)''' con un valor de '''2.6e74 enfermos'''.

-Para <math>h=0.01</math>
[[Archivo:H0.0110000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.0110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.35)''' con un valor de '''9521 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.31)''' con un valor de '''13550 enfermos'''.

-Para <math>h=0.001</math>
[[Archivo:H0.00110000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.00110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.342)''' con un valor de '''9470 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.338)''' con un valor de '''9758 enfermos'''.

-Para <math>h=0.0001</math>
[[Archivo:H0.000110000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.000110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.3411)''' con un valor de '''9465 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.3407)''' con un valor de '''9493 enfermos'''.

Para este segundo caso, nos damos cuenta que al utilizar tamaño de paso <math>h=0.1</math> , llegamos a un resultado sin sentido, pues nos da valores negativos y es imposible que una población tenga tales valores.
En las demás gráficas, los resultados obtenidos no parecen cuadrar con la realidad. Esto podría ser debido a que la alta diferencia entre los valores iniciales <math>(S _{0} , I_{0} )=(10000,40)</math> no permiten un buen cálculo de su variación respecto del tiempo.
El número inicial de infectados es muy elevado y al disminuir el tamaño de paso, el programa aproxima mejor la solución.

==Resolución del problema completo suponiendo <math> a </math> una función variable con el tiempo==
Para el siguiente caso supondremos que el coeficiente <math> a </math> (tasa de contagio en población sana) es una función dada por :

<math> a(t)= \frac{0.003}{(1+t)} </math>
Suponiendo además las condiciones iniciales <math> (S_{0},I_{0})=(1640,40) </math> , vamos a dibujar las gráficas de infectados y sanos por el método de Heun, tomando un tamaño de paso <math> h=0.01 </math>
[[Archivo: H0.01heunn.png|1000px|thumb|rigth| Población infectada-Población sana para a=a(t)]]
{{matlab|codigo=
%Heun. Coeficiente variable a(t)

%Datos
S0=input('Introducir valor inicial SO: ');
I0=input('Introducir valor inicial IO: ');
b=0.3;
c=0.01;
%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector solucion
y0=[S0;I0];
y=zeros(2,N+1);
y(:,1)=y0;
for i=1:N
K1S=-(0.003/(1+t(i)))*y(1,i)*y(2,i);
K2S=-(0.003/(1+t(i)+h))*(y(1,i)+K1S*h)*(y(2,i)+K1S*h);
y(1,i+1)=y(1,i)+(h/2)*(K1S+K2S);
K1I=(0.003/(1+t(i)))*y(1,i)*y(2,i)-b*y(2,i)-c*y(2,i);
K2I=(0.003/(1+t(i)+h))*(y(1,i)+K1I*h)*(y(2,i)+K1I*h)-b*(y(2,i)+K1I*h)-c*(y(2,i)+K1I*h);
y(2,i+1)=y(2,i)+(h/2)*(K1I+K2I);
end

%Tabla de resutados
[t',y(1,:)',y(2,:)']
%Graficos
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'k')
title('Método de Heun');
legend('Poblacion sana','Población infectada','Location','best');
hold off

% Vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h

}}

Como se puede ver la población sana se infecta rápidamente y a partir del día 15 casi toda la población ha fallecido. Al segundo día <math>(t=2,15)</math> observamos que el número de infectados es de 994 individuos <math>(I=994,3)</math>, siendo este el máximo que se alcanzará.

A medida que va pasando el tiempo el coeficiente <math> a(t) </math> (tasa de contagio en población sana) va a ir decreciendo poco a poco hasta aproximarse a cero según podemos ver en el siguiente gráfico:
<math> \lim_{t\to \infty } a(t)=0</math>
[[Archivo: A(t).png|450px|thumb|right|a(t) para t=[0,40]]]
{{matlab|codigo=

%tasa de contagio en poblacion sana a(t)

%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector a
a=zeros(1,N+1);
for i=1:N+1
a(i)=0.003/(1+t(i));
end

%Grafica
plot(t,a)

}}

[[Archivo: Comparacionheun.png|border|1500px|center|frame|Poblaciones según la tasa de contagio a(t)]]

Ambas poblaciones parten del mismo inicio (para <math>t=0</math> las poblaciones son iguales). Podemos observar que la población sana (susceptible de contraer la enfermedad) en el caso en el que la tasa de contagio en población sana <math>(a)</math> es constante <math>(a=0.003)</math> se ve reducida antes que en el caso en el que depende del tiempo <math>(a=a(t))</math>.

<math>a(t)</math> disminuye a la vez que transcurre el tiempo, por lo que en general las pendientes tanto de la población sana como de la infectada serán menores para este caso que para el caso en el que <math>a=0.003</math>

==Calibración del parámetro <math>a</math> en base a una experiencia previa==
Ahora vamos a intentar calibrar el párametro <math>a</math> para que el número de infectados tenga el máximo en el quinto día, para ello escribimos el siguiente código:

%Aplicamos el Método de Heunn

%Datos
t0=0;
tN=40;
h=0.1;
t=t0:h:tN;
a=0.0005:0.0001:0.002;
b=0.3;
c=0.01;
N1=round((tN-t0)/h);
N2=length(a);
Tmax=zeros(1,N2);

%Vectores iniciales de Población e Infectados
S=zeros(1,N1+1);
I=zeros(1,N1+1);

%Condiciones iniciales
S(1)=1600;
I(1)=40;

%Bucle que nos calcula todas las posibles soluciones
for j=1:N2
for i=1:N1
K1S=-a(j)*S(i)*I(i);
K2S=-a(j)*S(i)*I(i)+K1S*h;
S(i+1)=S(i)+h/2*(K1S+K2S);

K1I=I(i)*(a(j)*S(i)-(b+c));
K2I=I(i)*(a(j)*S(i)-(b+c))+K1I*h;
I(i+1)=I(i)+h/2*(K1I+K2I);

if i>1 && I(i+1)<I(i) && I(i)>I(i-1)
Tmax(j)=t(i);
end
end
end
distancia=(abs(Tmax-5));
[valormin,posmin]=min(distancia);
a=a(posmin)

Modelos epidemiológicos (Grupo 3A)

2016-05-02T18:27:03Z

Trabajocampos: /* Calibración del parámetro a en base a una experiencia previa */

{{ TrabajoED | Modelos Epidemológicos (Grupo 3A) | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED15/16|Curso 2015-16]] |
Ignacio Mollá Carcaño 
Pablo Revuelta Aragón 
David González Hernández 
Jose María García Rodríguez 
Alejandro Martínez Gamonal }}

En el desarrollo de una epidemia, se distinguen dos tipos de individuos. Los que ya han contraído la enfermedad (infectados), que llamaremos <math>I(t)</math> ; y los que son susceptibles de contraerla (sanos), a los que llamaremos <math>S(t)</math>. Donde <math>t</math> es la variable del tiempo. Se presentan dos hipótesis para realizar este estudio:

1. La población de personas infectadas se altera por el fallecimiento o la cura de las mismas. En ambos casos, la tasa de cambio depende del número 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ón entre el número de individuos en ambas clases.

A través del siguiente sistema de ecuaciones diferenciales en función del tiempo observamos las variaciones de ambas poblaciones:

<math>\begin{cases}\frac{dS}{dt}=-aSI\\\frac{dI}{dt}=aSI-bI-cI\end{cases}</math>

Donde <math>a</math>, <math>b</math> y <math>c</math> son parámetros.

* <math>\frac{dS}{dt}</math> : Es la variación de la población susceptible de contraer la enfermedad.

* <math>\frac{dI}{dt}</math> : Es la variación de la población infectada.

* <math>aSI</math> : Es la interacción entre ambas poblaciones que dependen de <math>a</math>.

* <math>bI</math> : Indica los individuos que fallecen. Por lo tanto <math>b</math> podría considerarse como la '''tasa de mortalidad'''.

* <math>cI</math> : Indica los individuos que se han curado. Por lo tanto <math>c</math> podría considerarse como la '''tasa de curación'''.

Para interpretar <math>a</math> observamos el sistema y la segunda hipótesis. Como los individuos que pasan de ser susceptibles a estar infectados son proporcionales a su interacción <math>(SI)</math>, <math>a</math> podría considerarse la '''tasa de contagio''' en la población susceptible.

==Resolución del problema con una sola ecuación diferencial==
Primero vamamos a estudiar el modelo epidemiológico de una población en la que incialmente las 2000 personas que conforman la población se encuentran infectadas. Para hacer el seguimiento de la enfermedad contamos tanto con la constante de proporcionalidad de personas que se curan (b=0.3) y la constante de proporcionalidad del número de personas que fallecen (c=0.01). De esta forma nos quedará el siguiente PVI:

<math> \begin{cases} \frac{dI}{dt} = -0.31 I \\I_{0} = 2000\end{cases} </math>

A continuación vamos a proceder a resolverlo por el método de Euler y trapecio asignandole un tamaño de paso de h=0.1 y veremos cuanto tiempo tardará en reducirse el número de infectados a la cuarta parte.

<math>\Rightarrow</math> Método de Euler
[[Archivo:Ejerciciodoseuler.png|500px|thumb|rigth|Aproximación mediante Euler]]
{{matlab|codigo=

clear all
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
y0=2000;
b=0.3;
c=0.01;

t(1)=t0;
ye(1)=y0; %Vector solución Euler

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de y(i) hasta que y(i)>500

%EULER

while ye(i)>500

ye(i+1)=ye(i)+h*(-b*ye(i)-c*ye(i));
t(i+1)=t(i)+h;
i=i+1;
end

[t',ye'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,ye)
legend('Población infectada(Método de Euler)','Location','best');

}}

<math>\Rightarrow</math> Método del trapecio
[[Archivo:Ejerciciodostrapecio.png|500px|thumb|rigth|Aproximación mediante el método del trapecio]]
{{matlab|codigo=
clear all
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
z0=2000;
b=0.3;
c=0.01;

t(1)=t0;
z(1)=z0; %Vector solución trapecio

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de z(i) hasta que z(i)>500

%TRAPECIO
while z(i)>500
z(i+1)=(z(i)*(1-(h/2)*(b+c)))/(1+h/2*(b+c));%trapecio
t(i+1)=t(i)+h;
i=i+1;
end

[t',z'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,z)
legend('Población infectada(Método del trapecio)','Location','best');

}}

Con todo esto podemos determinar que el número de infectados se reducirá a la cuarta parte en 4.5 días, esto ocurrirá o bien por el fallecimiento o por la curación de dicha enfermedad. Cabe señalar que el número de curados será mayor que el de defunciones, esto se puede ver a simple vista viendo que b>a.
De cara a nuestro estudio de esta enfermedad vamos a ver como afectaría al desarrollo de la misma si introducimos 100 sujetos sanos en la población, tomando una constante de interacción entre personas infectadas y sanas de a=0.003. Procederemos igual que antes a resolverlo por el método de Euler y trapecio con un tamaño de paso de 0.1.
Primero vamos a definir nuestro nuevo PVI:
<math> \begin{cases} \frac{dI}{dt} = -0.01 I \\I_{0} = 2000\end{cases} </math>
<math>\Rightarrow</math> Método de Euler
[[Archivo:Ejerciciotreseuler.png|500px|thumb|rigth|Aproximación mediante Euler para S=100]]
{{matlab|codigo=

%Euler
clear all
%Datos
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%a: Es la constante de proporcionalidad entre los infectados por cada
%interacción
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
y0=2000;
S=100;
a=0.003;
b=0.3;
c=0.01;

t(1)=t0;
y(1)=y0;

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de y(i) hasta que y(i)>500
while y(i)>500

y(i+1)=y(i)+h*(a*S*y(i)-b*y(i)-c*y(i));
t(i+1)=t(i)+h;
i=i+1;
end

[t',y'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,y)
legend('Población infectada para S=100(Método de Euler)','Location','best');

}}

<math>\Rightarrow</math> Método del trapecio
[[Archivo:Ejerciciotrestrapecio.png|500px|thumb|rigth|Aproximación mediante el método del trapecio para S=100]]
{{matlab|codigo=

%Trapecio
clear all
%Datos
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%a: Es la constante de proporcionalidad entre los infectados por cada
%interacción
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
z0=2000;
S=100;
a=0.003;
b=0.3;
c=0.01;

t(1)=t0;
z(1)=z0;

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de z(i) hasta que z(i)>500

%TRAPECIO
while z(i)>500
z(i+1)=(z(i)*(1-(h/2)*(b+c-S*a)))/(1+h/2*(b+c-S*a));
t(i+1)=t(i)+h;
i=i+1;
end

[t',z'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,z)
legend('Población infectada para S=100(Método del trapecio)','Location','best');

}}

Ahora con estos nuevos datos vemos que para que se reduzca a la cuarta parte el número de infectados tendrán que pasar 138.7 días, algo lógico al fijarnos en la ED, ambas son ecuaciones decrecientes pero para el caso de S=100 la pendiente será mucho menor.
Si continuamos haciendo pruebas vemos para valores superiores de s, por ejemplo, S=200 el matlab no nos da ninguna respuesta. Esto ocurre porque al aumentar S la ecuación pasa a ser creciente, por lo que lógicamente si intentas buscar un valor del número de infectados menor al inicial no lo podrás encontrar.

==Resolución del problema completo==
Vamos a resolver el problema completo para ver como evolucionarían las poblaciones en un periodo de 40 días.
Primero supongamos una población infectada inicial de 20 individuos y una población en riesgo de contagio de 800.
Después supondremos el caso de 40 individuos enfermos y 10000 en riesgo de contagio.
Para resolverlo usaremos el método de Euler y más tarde el método de Runge-Kutta de orden cuatro con distintas discretizaciones para ver la influencia de estas.

<math>\Rightarrow</math> Método de Euler
{{matlab|codigo=

%Euler
clear all
%Datos
%Introducimos los valores de las constantes
a=0.003;
b=0.3;
c=0.01;
h=input('Introducir valores de h:');
t0=0;
y0=input('Introducir vector [So,Io]:');
tN=40;
% calculamos los subintervalos
t=t0:h:tN;
N=(tN-t0)/h;
y=zeros(2,length(t));
y(:,1)=y0';
for i= 1:N;
y(:,i+1)=y(:,i)+h*[-a*y(1,i).*y(2,i);a*y(1,i).*y(2,i)-b*y(2,i)-c*y(2,i)];
end

% grafico
plot(t,y)
% vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h
legend('Población sana','Location','best','Población enferma','Location','best');

}}

<math>\Rightarrow</math> Método de Runge-Kutta

{{matlab|codigo=

%Runge-Kutta
clear all

%Datos
S0=input('Introducir valor inicial SO: ');
I0=input('Introducir valor inicial IO: ');
a=0.003;
b=0.3;
c=0.01;
%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector solucion
y0=[S0;I0];
y=zeros(2,N+1);
y(:,1)=y0;
for i=1:N
K1S=-a*y(1,i)*y(2,i);
K2S=-a*(y(1,i)+(1/2)*K1S*h)*(y(2,i)+(1/2)*K1S*h);
K3S=-a*(y(1,i)+(1/2)*K2S*h)*(y(2,i)+(1/2)*K2S*h);
K4S=-a*(y(1,i)+K3S*h)*(y(2,i)+K3S*h);
y(1,i+1)=y(1,i)+(h/6)*(K1S+2*K2S+2*K3S+K4S);
K1I=a*y(1,i)*y(2,i)-b*y(2,i)-c*y(2,i);
K2I=a*(y(1,i)+(1/2)*K1I*h)*(y(2,i)+(1/2)*K1I*h)-b*(y(2,i)+(1/2)*K1I*h)-c*(y(2,i)+(1/2)*K1I*h);
K3I=a*(y(1,i)+(1/2)*K2I*h)*(y(2,i)+(1/2)*K2I*h)-b*(y(2,i)+(1/2)*K2I*h)-c*(y(2,i)+(1/2)*K2I*h);
K4I=a*(y(1,i)+K3I*h)*(y(2,i)+K3I*h)-b*(y(2,i)+K3I*h)-c*(y(2,i)+K3I*h);
y(2,i+1)=y(2,i)+(h/6)*(K1I+2*K2I+2*K3I+K4I);
end

%Tabla de resutados
[t',y(1,:)',y(2,:)']
%Gráfico
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'k')
legend('Poblacion sana','Poblacion enferma','Location','best');
hold off

% vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h

}}

La dificultad que hay en el uso del método del trapecio o cualquier otro método implícito está en que las incógnitas de nuestro sistema <math>(S,I)</math> quedan implícitas en ambas ecuaciones y necesitan que se despejen. Al comenzar a realizar cálculos nos damos cuenta que el despeje es imposible, ya que las dos variables incógnita aparecen multiplicándose entre sí, quedando una variable dependiendo de la otra siempre.

===Situación inicial de <math>(S _{0} , I_{0} )=(800,20)</math>===
En estos programas, sustituyendo las <math>h</math> para 0.1, 0.01, 0.001 y 0.0001 obtendremos los siguientes resultados:

-Para <math>h=0.1</math>

[[Archivo:Tborraré.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.1rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

Como se puede ver la población sana se infecta rápidamente y a partir del día 15 prácticamente toda la población ha fallecido.
* '''Método de Euler''': nuestro programa nos indica que en este caso el momento en el que los enfermos son máximos se produce al '''tercer día (2.8)''' y el '''número de enfermos en ese momento es 517'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce al '''tercer día (2.6)''' y el '''número de enfermos en ese momento es 650'''.

-Para <math>h=0.01</math>
[[Archivo:Tborraré2.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.01rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

En el gráfico apenas se encuentra diferencia, lo cual es lógico porque es el mismo problema, pero si sobrepusiéramos las gráficas, lograríamos ver como se va ajustando a la realidad en función del aumento de h.
* '''Método de Euler''': en este caso el momento en el que los enfermos son máximos es también en el '''tercer día (2.7)''' y el número de '''enfermos es 506'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce al '''tercer día (2.69)''' y el número de '''enfermos en ese momento es 517'''.

-Para <math>h=0.001</math>
[[Archivo:Tborraré3.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.001rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': en este caso el momento en el que los enfermos son máximos se produce el '''tercer día (2.696)''' con un valor de '''505 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''tercer día (2.694)''' con un valor de '''506 enfermos'''.

-Para <math>h=0.0001</math>
[[Archivo:Tborraré4.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.0001rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': al igual que antes, se repite el momento en el que los enfermos son máximos en el '''tercer día (2.695)''' con un valor de '''505 enfermos'''.
* '''Método de Runge-Kutta''': se repite el momento en el que los enfermos son máximos en el '''tercer día (2.6948)''' con un valor de '''505 enfermos'''.

Observamos en las distintas gráficas obtenidas, que según vamos empleando un tamaño de paso <math>h</math> más pequeño, nos vamos acercando cada vez más a la población real que pretendemos reflejar.
Nos damos cuenta de que las gráficas tienen bastante lógica, pues a la vez que desciende el número de población sana va aumentando la población infectada.
También observamos que para el método de Runge-Kutta cuando utilizamos tamaño de paso <math>h=0.1</math> pierde precisión en relación a cuando utilizamos tamaños de paso más pequeños.

===Situación inicial de <math>(S _{0} , I_{0} )=(10000,40)</math>===
En estos programas, sustituyendo las <math>h</math> para 0.1, 0.01, 0.001 y 0.0001 obtendremos los siguientes resultados:

-Para <math>h=0.1</math>

[[Archivo:H0.11000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.5)''' con un valor de '''12600 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.4)''' con un valor de '''2.6e74 enfermos'''.

-Para <math>h=0.01</math>
[[Archivo:H0.0110000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.0110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.35)''' con un valor de '''9521 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.31)''' con un valor de '''13550 enfermos'''.

-Para <math>h=0.001</math>
[[Archivo:H0.00110000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.00110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.342)''' con un valor de '''9470 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.338)''' con un valor de '''9758 enfermos'''.

-Para <math>h=0.0001</math>
[[Archivo:H0.000110000e.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.000110000r.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* '''Método de Euler''': el momento en el que los enfermos son máximos se produce el '''primer día (0.3411)''' con un valor de '''9465 enfermos'''.
* '''Método de Runge-Kutta''': el momento en el que los enfermos son máximos se produce el '''primer día (0.3407)''' con un valor de '''9493 enfermos'''.

Para este segundo caso, nos damos cuenta que al utilizar tamaño de paso <math>h=0.1</math> , llegamos a un resultado sin sentido, pues nos da valores negativos y es imposible que una población tenga tales valores.
En las demás gráficas, los resultados obtenidos no parecen cuadrar con la realidad. Esto podría ser debido a que la alta diferencia entre los valores iniciales <math>(S _{0} , I_{0} )=(10000,40)</math> no permiten un buen cálculo de su variación respecto del tiempo.
El número inicial de infectados es muy elevado y al disminuir el tamaño de paso, el programa aproxima mejor la solución.

==Resolución del problema completo suponiendo <math> a </math> una función variable con el tiempo==
Para el siguiente caso supondremos que el coeficiente <math> a </math> (tasa de contagio en población sana) es una función dada por :

<math> a(t)= \frac{0.003}{(1+t)} </math>
Suponiendo además las condiciones iniciales <math> (S_{0},I_{0})=(1640,40) </math> , vamos a dibujar las gráficas de infectados y sanos por el método de Heun, tomando un tamaño de paso <math> h=0.01 </math>
[[Archivo: H0.01heunn.png|1000px|thumb|rigth| Población infectada-Población sana para a=a(t)]]
{{matlab|codigo=
%Heun. Coeficiente variable a(t)

%Datos
S0=input('Introducir valor inicial SO: ');
I0=input('Introducir valor inicial IO: ');
b=0.3;
c=0.01;
%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector solucion
y0=[S0;I0];
y=zeros(2,N+1);
y(:,1)=y0;
for i=1:N
K1S=-(0.003/(1+t(i)))*y(1,i)*y(2,i);
K2S=-(0.003/(1+t(i)+h))*(y(1,i)+K1S*h)*(y(2,i)+K1S*h);
y(1,i+1)=y(1,i)+(h/2)*(K1S+K2S);
K1I=(0.003/(1+t(i)))*y(1,i)*y(2,i)-b*y(2,i)-c*y(2,i);
K2I=(0.003/(1+t(i)+h))*(y(1,i)+K1I*h)*(y(2,i)+K1I*h)-b*(y(2,i)+K1I*h)-c*(y(2,i)+K1I*h);
y(2,i+1)=y(2,i)+(h/2)*(K1I+K2I);
end

%Tabla de resutados
[t',y(1,:)',y(2,:)']
%Graficos
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'k')
title('Método de Heun');
legend('Poblacion sana','Población infectada','Location','best');
hold off

% Vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h

}}

Como se puede ver la población sana se infecta rápidamente y a partir del día 15 casi toda la población ha fallecido. Al segundo día <math>(t=2,15)</math> observamos que el número de infectados es de 994 individuos <math>(I=994,3)</math>, siendo este el máximo que se alcanzará.

A medida que va pasando el tiempo el coeficiente <math> a(t) </math> (tasa de contagio en población sana) va a ir decreciendo poco a poco hasta aproximarse a cero según podemos ver en el siguiente gráfico:
<math> \lim_{t\to \infty } a(t)=0</math>
[[Archivo: A(t).png|450px|thumb|right|a(t) para t=[0,40]]]
{{matlab|codigo=

%tasa de contagio en poblacion sana a(t)

%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector a
a=zeros(1,N+1);
for i=1:N+1
a(i)=0.003/(1+t(i));
end

%Grafica
plot(t,a)

}}

[[Archivo: Comparacionheun.png|border|1500px|center|frame|Poblaciones según la tasa de contagio a(t)]]

Ambas poblaciones parten del mismo inicio (para <math>t=0</math> las poblaciones son iguales). Podemos observar que la población sana (susceptible de contraer la enfermedad) en el caso en el que la tasa de contagio en población sana <math>(a)</math> es constante <math>(a=0.003)</math> se ve reducida antes que en el caso en el que depende del tiempo <math>(a=a(t))</math>.

<math>a(t)</math> disminuye a la vez que transcurre el tiempo, por lo que en general las pendientes tanto de la población sana como de la infectada serán menores para este caso que para el caso en el que <math>a=0.003</math>

==Calibración del parámetro <math>a</math> en base a una experiencia previa==
Ahora vamos a intentar calibrar el párametro <math>a</math> para que el número de infectados tenga el máximo en el quinto día, para ello escribimos el siguiente código:

%Datos
t0=0;
tN=40;
h=0.1;
t=t0:h:tN;
a=0.0005:0.0001:0.002;
b=0.3;
c=0.01;
N1=round((tN-t0)/h);
N2=length(a);
Tmax=zeros(1,N2);

%Vectores iniciales de Población e Infectados
S=zeros(1,N1+1);
I=zeros(1,N1+1);

%Condiciones iniciales
S(1)=1600;
I(1)=40;

%Bucle que nos calcula todas las posibles soluciones
for j=1:N2
for i=1:N1
K1S=-a(j)*S(i)*I(i);
K2S=-a(j)*S(i)*I(i)+K1S*h;
S(i+1)=S(i)+h/2*(K1S+K2S);

K1I=I(i)*(a(j)*S(i)-(b+c));
K2I=I(i)*(a(j)*S(i)-(b+c))+K1I*h;
I(i+1)=I(i)+h/2*(K1I+K2I);

if i>1 && I(i+1)<I(i) && I(i)>I(i-1)
Tmax(j)=t(i);
end
end
end
distancia=(abs(Tmax-5));
[valormin,posmin]=min(distancia);
a=a(posmin)

Modelos epidemiológicos (Grupo 3A)

2016-05-02T09:14:23Z

Trabajocampos: /* Introducción y planteamiento del problema */

{{ TrabajoED | Modelos Epidemológicos (Grupo 3A) | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED15/16|Curso 2015-16]] |
Ignacio Mollá Carcaño 
Pablo Revuelta Aragón 
David González Hernández 
Jose María García Rodríguez 
Alejandro Martínez Gamonal }}

En el desarrollo de una epidemia, se distinguen dos tipos de individuos. Los que ya han contraído la enfermedad (infectados), que llamaremos I(t); y los que son susceptibles de contraerla, a los que llamaremos S(t). Donde t es la variable del tiempo. Se presentan dos hipótesis para realizar este estudio:

1. La población de personas infectadas se altera por el fallecimiento o la cura de las mismas. En ambos casos, la tasa de cambio depende del número 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ón entre el número de individuos en ambas clases.

A través del siguiente sistema de ecuaciones diferenciales en función del tiempo observamos las variaciones de ambas poblaciones:

dS/dt =-aSI
dI/dt =aSI-bI-cI
Donde a, b y c son parámetros.

• dS/dt : Es la variación de la población susceptible de contraer la enfermedad.
• dI/dt : Es la variación de la población infectada.
• aSI : Es la interacción entre ambas poblaciones que dependen de a.
• bI : Indica los individuos que fallecen. Por lo tanto 'b' podría considerarse como la '''tasa de mortalidad'''.
• cI : Indica los individuos que se han curado. Por lo tanto 'c' podría considerarse como la '''tasa de curación'''.

Para interpretar 'a' observamos el sistema y la segunda hipótesis. Como los individuos que pasan de ser susceptibles a estar infectados son proporcionales a su interacción(SI), 'a' podría considerarse la '''tasa de contagio''' en la población susceptible.

==Resolución del problema con una sola ecuación diferencial==
Primero vamamos a estudiar el modelo epidemiológico de una población en la que incialmente las 2000 personas que conforman la población se encuentran infectadas. Para hacer el seguimiento de la enfermedad contamos tanto con la constante de proporcionalidad de personas que se curan (b=0.3) y la constante de proporcionalidad del número de personas que fallecen (c=0.01). De esta forma nos quedará el siguiente PVI:

<math> \begin{cases} \frac{dI}{dt} = -0.31 I \\I_{0} = 2000\end{cases} </math>

A continuación vamos a proceder a resolverlo por el método de Euler y trapecio asignandole un tamaño de paso de h=0.1 y veremos cuanto tiempo tardará en reducirse el número de infectados a la cuarta parte.

<math>\Rightarrow</math> Método de Euler
[[Archivo:Ejerciciodoseuler.png|500px|thumb|rigth|Aproximación mediante Euler]]
{{matlab|codigo=

clear all
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
y0=2000;
b=0.3;
c=0.01;

t(1)=t0;
ye(1)=y0; %Vector solución Euler

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de y(i) hasta que y(i)>500

%EULER

while ye(i)>500

ye(i+1)=ye(i)+h*(-b*ye(i)-c*ye(i));
t(i+1)=t(i)+h;
i=i+1;
end

[t',ye'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,ye)
legend('Población infectada(Método de Euler)','Location','best');

}}

<math>\Rightarrow</math> Método del trapecio
[[Archivo:Ejerciciodostrapecio.png|500px|thumb|rigth|Aproximación mediante el método del trapecio]]
{{matlab|codigo=
clear all
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
z0=2000;
b=0.3;
c=0.01;

t(1)=t0;
z(1)=z0; %Vector solución trapecio

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de z(i) hasta que z(i)>500

%TRAPECIO
while z(i)>500
z(i+1)=(z(i)*(1-(h/2)*(b+c)))/(1+h/2*(b+c));%trapecio
t(i+1)=t(i)+h;
i=i+1;
end

[t',z'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,z)
legend('Población infectada(Método del trapecio)','Location','best');

}}

Con todo esto podemos determinar que el número de infectados se reducirá a la cuarta parte en 4.5 días, esto ocurrirá o bien por el fallecimiento o por la curación de dicha enfermedad. Cabe señalar que el número de curados será mayor que el de defunciones, esto se puede ver a simple vista viendo que b>a.
De cara a nuestro estudio de esta enfermedad vamos a ver como afectaría al desarrollo de la misma si introducimos 100 sujetos sanos en la población, tomando una constante de interacción entre personas infectadas y sanas de a=0.003. Procederemos igual que antes a resolverlo por el método de Euler y trapecio con un tamaño de paso de 0.1.
Primero vamos a definir nuestro nuevo PVI:
<math> \begin{cases} \frac{dI}{dt} = -0.01 I \\I_{0} = 2000\end{cases} </math>
<math>\Rightarrow</math> Método de Euler
[[Archivo:Ejerciciotreseuler.png|500px|thumb|rigth|Aproximación mediante Euler para S=100]]
{{matlab|codigo=

%Euler
clear all
%Datos
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%a: Es la constante de proporcionalidad entre los infectados por cada
%interacción
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
y0=2000;
S=100;
a=0.003;
b=0.3;
c=0.01;

t(1)=t0;
y(1)=y0;

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de y(i) hasta que y(i)>500
while y(i)>500

y(i+1)=y(i)+h*(a*S*y(i)-b*y(i)-c*y(i));
t(i+1)=t(i)+h;
i=i+1;
end

[t',y'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,y)
legend('Población infectada para S=100(Método de Euler)','Location','best');

}}

<math>\Rightarrow</math> Método del trapecio
[[Archivo:Ejerciciotrestrapecio.png|500px|thumb|rigth|Aproximación mediante el método del trapecio para S=100]]
{{matlab|codigo=

%Trapecio
clear all
%Datos
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%a: Es la constante de proporcionalidad entre los infectados por cada
%interacción
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
z0=2000;
S=100;
a=0.003;
b=0.3;
c=0.01;

t(1)=t0;
z(1)=z0;

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de z(i) hasta que z(i)>500

%TRAPECIO
while z(i)>500
z(i+1)=(z(i)*(1-(h/2)*(b+c-S*a)))/(1+h/2*(b+c-S*a));
t(i+1)=t(i)+h;
i=i+1;
end

[t',z'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,z)
legend('Población infectada para S=100(Método del trapecio)','Location','best');

}}

Ahora con estos nuevos datos vemos que para que se reduzca a la cuarta parte el número de infectados tendrán que pasar 138.7 días, algo lógico al fijarnos en la ED, ambas son ecuaciones decrecientes pero para el caso de S=100 la pendiente será mucho menor.
Si continuamos haciendo pruebas vemos para valores superiores de s, por ejemplo, S=200 el matlab no nos da ninguna respuesta. Esto ocurre porque al aumentar S la ecuación pasa a ser creciente, por lo que lógicamente si intentas buscar un valor del número de infectados menor al inicial no lo podrás encontrar.

==Resolución del problema completo==
Vamos a resolver el problema completo para ver como evolucionarían las poblaciones en un periodo de 40 días.
Primero supongamos una población infectada inicial de 20 individuos y una población en riesgo de contagio de 800.
Después supondremos el caso de 40 individuos enfermos y 10000 en riesgo de contagio.
Para resolverlo usaremos el método de Euler y más tarde el método de Runge-Kutta de orden cuatro con distintas discretizaciones para ver la influencia de estas.

<math>\Rightarrow</math> Método de Euler
{{matlab|codigo=

%Euler
clear all
%Datos
%Introducimos los valores de las constantes
a=0.003;
b=0.3;
c=0.01;
h=input('Introducir valores de h:');
t0=0;
y0=input('Introducir vector [So,Io]:');
tN=40;
% calculamos los subintervalos
t=t0:h:tN;
N=(tN-t0)/h;
y=zeros(2,length(t));
y(:,1)=y0';
for i= 1:N;
y(:,i+1)=y(:,i)+h*[-a*y(1,i).*y(2,i);a*y(1,i).*y(2,i)-b*y(2,i)-c*y(2,i)];
end

% grafico
plot(t,y)
% vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h
legend('Población sana','Location','best','Población enferma','Location','best');

}}

<math>\Rightarrow</math> Método de Runge-Kutta

{{matlab|codigo=

%Runge-Kutta
clear all

%Datos
S0=input('Introducir valor inicial SO: ');
I0=input('Introducir valor inicial IO: ');
a=0.003;
b=0.3;
c=0.01;
%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector solucion
y0=[S0;I0];
y=zeros(2,N+1);
y(:,1)=y0;
for i=1:N
K1S=-a*y(1,i)*y(2,i);
K2S=-a*(y(1,i)+(1/2)*K1S*h)*(y(2,i)+(1/2)*K1S*h);
K3S=-a*(y(1,i)+(1/2)*K2S*h)*(y(2,i)+(1/2)*K2S*h);
K4S=-a*(y(1,i)+K3S*h)*(y(2,i)+K3S*h);
y(1,i+1)=y(1,i)+(h/6)*(K1S+2*K2S+2*K3S+K4S);
K1I=a*y(1,i)*y(2,i)-b*y(2,i)-c*y(2,i);
K2I=a*(y(1,i)+(1/2)*K1I*h)*(y(2,i)+(1/2)*K1I*h)-b*(y(2,i)+(1/2)*K1I*h)-c*(y(2,i)+(1/2)*K1I*h);
K3I=a*(y(1,i)+(1/2)*K2I*h)*(y(2,i)+(1/2)*K2I*h)-b*(y(2,i)+(1/2)*K2I*h)-c*(y(2,i)+(1/2)*K2I*h);
K4I=a*(y(1,i)+K3I*h)*(y(2,i)+K3I*h)-b*(y(2,i)+K3I*h)-c*(y(2,i)+K3I*h);
y(2,i+1)=y(2,i)+(h/6)*(K1I+2*K2I+2*K3I+K4I);
end

%Tabla de resutados
[t',y(1,:)',y(2,:)']
%Gráfico
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'k')
legend('Poblacion sana','Poblacion enferma','Location','best');
hold off

% vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h

}}

La dificultad que hay en el uso del método del trapecio o cualquier otro método implícito está en que las incógnitas de nuestro sistema <math>(S,I)</math> quedan implícitas en ambas ecuaciones y necesitan que se despejen. Al comenzar a realizar cálculos nos damos cuenta que el despeje es imposible, ya que las dos variables incógnita aparecen multiplicándose entre sí, quedando una variable dependiendo de la otra siempre.

===Situación inicial de <math>(S _{0} , I_{0} )=(800,20)</math>===
En estos programas, sustituyendo las <math>h</math> para 0.1, 0.01, 0.001 y 0.0001 obtendremos los siguientes resultados:

-Para <math>h=0.1</math>

[[Archivo:Tborraré.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.1rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

Como se puede ver la población sana se infecta rápidamente y a partir del día 15 prácticamente toda la población ha fallecido.
* Método de Euler: nuestro programa nos indica que en este caso el momento en el que los enfermos son máximos se produce al tercer día (Día 2.8) y el número de enfermos en ese momento es 517.
* Método de Runge-Kutta: el momento en el que los enfermos son máximos se produce al tercer día (Día 2.6) y el número de enfermos en ese momento es 650.

-Para <math>h=0.01</math>
[[Archivo:Tborraré2.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.01rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

En el gráfico apenas se encuentra diferencia, lo cual es lógico porque es el mismo problema, pero si sobrepusiéramos las gráficas, lograríamos ver como se va ajustando a la realidad en función del aumento de h.
* Método de Euler: en este caso el momento en el que los enfermos son máximos es también en el tercer día (2.7) y el número de enfermos es 506.
* Método de Runge-Kutta: el momento en el que los enfermos son máximos se produce al tercer día (Día 2.69) y el número de enfermos en ese momento es 517.

-Para <math>h=0.001</math>
[[Archivo:Tborraré3.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.001rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* Método de Euler: en este caso el momento en el que los enfermos son máximos se produce el tercer día (2.696) con un valor de 505 enfermos.
* Método de Runge-Kutta: el momento en el que los enfermos son máximos se produce el tercer día (2.694) con un valor de 506 enfermos.

-Para <math>h=0.0001</math>
[[Archivo:Tborraré4.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.0001rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* Método de Euler: al igual que antes, se repite el momento en el que los enfermos son máximos en el tercer día (2.695) con un valor de 505 enfermos.
* Método de Runge-Kutta: se repite el momento en el que los enfermos son máximos en el tercer día (2.6948) con un valor de 505 enfermos.

===Situación inicial de <math>(S _{0} , I_{0} )=(10000,40)</math>===
En estos programas, sustituyendo las <math>h</math> para 0.1, 0.01, 0.001 y 0.0001 obtendremos los siguientes resultados:

-Para <math>h=0.1</math>

[[Archivo:|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* Método de Euler:0.5 1.26e4
* Método de Runge-Kutta: 0.4 2.596e74

-Para <math>h=0.01</math>
[[Archivo:|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* Método de Euler:0.35 9.5211e3
* Método de Runge-Kutta: 0.3100 1.3546e4

-Para <math>h=0.001</math>
[[Archivo:|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* Método de Euler: 0.3420 9.4696e3
* Método de Runge-Kutta: 0.3380 9.5781e3
9
-Para <math>h=0.0001</math>
[[Archivo:|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* Método de Euler: 0.3411 9.4647e3
* Método de Runge-Kutta: 0.347 9.4298e3

==Resolución del problema completo suponiendo <math> a </math> una función variable con el tiempo==
Para el siguiente caso supondremos que el coeficiente <math> a </math> (tasa de contagio en población sana) es una función dada por :

<math> a(t)= \frac{0.003}{(1+t)} </math>
Suponiendo además las condiciones iniciales <math> (S_{0},I_{0})=(1640,40) </math> , vamos a dibujar las gráficas de infectados y sanos por el método de Heun, tomando un tamaño de paso <math> h=0.01 </math>
[[Archivo: H0.01heunn.png|1000px|thumb|rigth| Población infectada-Población sana para a=a(t)]]
{{matlab|codigo=
%Heun. Coeficiente variable a(t)

%Datos
S0=input('Introducir valor inicial SO: ');
I0=input('Introducir valor inicial IO: ');
b=0.3;
c=0.01;
%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector solucion
y0=[S0;I0];
y=zeros(2,N+1);
y(:,1)=y0;
for i=1:N
K1S=-(0.003/(1+t(i)))*y(1,i)*y(2,i);
K2S=-(0.003/(1+t(i)+h))*(y(1,i)+K1S*h)*(y(2,i)+K1S*h);
y(1,i+1)=y(1,i)+(h/2)*(K1S+K2S);
K1I=(0.003/(1+t(i)))*y(1,i)*y(2,i)-b*y(2,i)-c*y(2,i);
K2I=(0.003/(1+t(i)+h))*(y(1,i)+K1I*h)*(y(2,i)+K1I*h)-b*(y(2,i)+K1I*h)-c*(y(2,i)+K1I*h);
y(2,i+1)=y(2,i)+(h/2)*(K1I+K2I);
end

%Tabla de resutados
[t',y(1,:)',y(2,:)']
%Graficos
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'k')
title('Método de Heun');
legend('Poblacion sana','Población infectada','Location','best');
hold off

% Vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h

}}

Como se puede ver la población sana se infecta rápidamente y a partir del día 15 casi toda la población ha fallecido. Al segundo día <math>(t=2,15)</math> observamos que el número de infectados es de 994 individuos <math>(I=994,3)</math>, siendo este el máximo que se alcanzará.

A medida que va pasando el tiempo el coeficiente <math> a(t) </math> (tasa de contagio en población sana) va a ir decreciendo poco a poco hasta aproximarse a cero según podemos ver en el siguiente gráfico:
<math> \lim_{t\to \infty } a(t)=0</math>
[[Archivo: A(t).png|450px|thumb|right|a(t) para t=[0,40]]]
{{matlab|codigo=

%tasa de contagio en poblacion sana a(t)

%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector a
a=zeros(1,N+1);
for i=1:N+1
a(i)=0.003/(1+t(i));
end

%Grafica
plot(t,a)

}}

[[Archivo: Comparacionheun.png|border|1500px|center|frame|Poblaciones según la tasa de contagio a(t)]]

Ambas poblaciones parten del mismo inicio (para <math>t=0</math> las poblaciones son iguales). Podemos observar que la población sana (susceptible de contraer la enfermedad) en el caso en el que la tasa de contagio en población sana <math>(a)</math> es constante <math>(a=0.003)</math> se ve reducida antes que en el caso en el que depende del tiempo <math>(a=a(t))</math>.

<math>a(t)</math> disminuye a la vez que transcurre el tiempo, por lo que en general las pendientes tanto de la población sana como de la infectada serán menores para este caso que para el caso en el que <math>a=0.003</math>

==Calibración del parámetro <math>a</math> en base a una experiencia previa==

Modelos epidemiológicos (Grupo 3A)

2016-05-01T21:35:49Z

Trabajocampos: /* Situación inicial de (S _{0} , I_{0} )=(10000,40) */

{{ TrabajoED | Modelos Epidemológicos (Grupo 3A) | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED15/16|Curso 2015-16]] |
Ignacio Mollá Carcaño 
Pablo Revuelta Aragón 
David González Hernández 
Jose María García Rodríguez 
Alejandro Martínez Gamonal }}

==Introducción y planteamiento del problema==

==Resolución del problema con una sola ecuación diferencial==
Primero vamamos a estudiar el modelo epidemiológico de una población en la que incialmente las 2000 personas que conforman la población se encuentran infectadas. Para hacer el seguimiento de la enfermedad contamos tanto con la constante de proporcionalidad de personas que se curan (b=0.3) y la constante de proporcionalidad del número de personas que fallecen (c=0.01). De esta forma nos quedará el siguiente PVI:

<math> \begin{cases} \frac{dI}{dt} = -0.31 I \\I_{0} = 2000\end{cases} </math>

A continuación vamos a proceder a resolverlo por el método de Euler y trapecio asignandole un tamaño de paso de h=0.1 y veremos cuanto tiempo tardará en reducirse el número de infectados a la cuarta parte.

<math>\Rightarrow</math> Método de Euler
[[Archivo:Ejerciciodoseuler.png|500px|thumb|rigth|Aproximación mediante Euler]]
{{matlab|codigo=

clear all
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
y0=2000;
b=0.3;
c=0.01;

t(1)=t0;
ye(1)=y0; %Vector solución Euler

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de y(i) hasta que y(i)>500

%EULER

while ye(i)>500

ye(i+1)=ye(i)+h*(-b*ye(i)-c*ye(i));
t(i+1)=t(i)+h;
i=i+1;
end

[t',ye'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,ye)
legend('Población infectada(Método de Euler)','Location','best');

}}

<math>\Rightarrow</math> Método del trapecio
[[Archivo:Ejerciciodostrapecio.png|500px|thumb|rigth|Aproximación mediante el método del trapecio]]
{{matlab|codigo=
clear all
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
z0=2000;
b=0.3;
c=0.01;

t(1)=t0;
z(1)=z0; %Vector solución trapecio

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de z(i) hasta que z(i)>500

%TRAPECIO
while z(i)>500
z(i+1)=(z(i)*(1-(h/2)*(b+c)))/(1+h/2*(b+c));%trapecio
t(i+1)=t(i)+h;
i=i+1;
end

[t',z'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,z)
legend('Población infectada(Método del trapecio)','Location','best');

}}

Con todo esto podemos determinar que el número de infectados se reducirá a la cuarta parte en 4.5 días, esto ocurrirá o bien por el fallecimiento o por la curación de dicha enfermedad. Cabe señalar que el número de curados será mayor que el de defunciones, esto se puede ver a simple vista viendo que b>a.
De cara a nuestro estudio de esta enfermedad vamos a ver como afectaría al desarrollo de la misma si introducimos 100 sujetos sanos en la población, tomando una constante de interacción entre personas infectadas y sanas de a=0.003. Procederemos igual que antes a resolverlo por el método de Euler y trapecio con un tamaño de paso de 0.1.
Primero vamos a definir nuestro nuevo PVI:
<math> \begin{cases} \frac{dI}{dt} = -0.01 I \\I_{0} = 2000\end{cases} </math>
<math>\Rightarrow</math> Método de Euler
[[Archivo:Ejerciciotreseuler.png|500px|thumb|rigth|Aproximación mediante Euler para S=100]]
{{matlab|codigo=

%Euler
clear all
%Datos
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%a: Es la constante de proporcionalidad entre los infectados por cada
%interacción
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
y0=2000;
S=100;
a=0.003;
b=0.3;
c=0.01;

t(1)=t0;
y(1)=y0;

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de y(i) hasta que y(i)>500
while y(i)>500

y(i+1)=y(i)+h*(a*S*y(i)-b*y(i)-c*y(i));
t(i+1)=t(i)+h;
i=i+1;
end

[t',y'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,y)
legend('Población infectada para S=100(Método de Euler)','Location','best');

}}

<math>\Rightarrow</math> Método del trapecio
[[Archivo:Ejerciciotrestrapecio.png|500px|thumb|rigth|Aproximación mediante el método del trapecio para S=100]]
{{matlab|codigo=

%Trapecio
clear all
%Datos
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%a: Es la constante de proporcionalidad entre los infectados por cada
%interacción
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
z0=2000;
S=100;
a=0.003;
b=0.3;
c=0.01;

t(1)=t0;
z(1)=z0;

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de z(i) hasta que z(i)>500

%TRAPECIO
while z(i)>500
z(i+1)=(z(i)*(1-(h/2)*(b+c-S*a)))/(1+h/2*(b+c-S*a));
t(i+1)=t(i)+h;
i=i+1;
end

[t',z'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,z)
legend('Población infectada para S=100(Método del trapecio)','Location','best');

}}

Ahora con estos nuevos datos vemos que para que se reduzca a la cuarta parte el número de infectados tendrán que pasar 138.7 días, algo lógico al fijarnos en la ED, ambas son ecuaciones decrecientes pero para el caso de S=100 la pendiente será mucho menor.
Si continuamos haciendo pruebas vemos para valores superiores de s, por ejemplo, S=200 el matlab no nos da ninguna respuesta. Esto ocurre porque al aumentar S la ecuación pasa a ser creciente, por lo que lógicamente si intentas buscar un valor del número de infectados menor al inicial no lo podrás encontrar.

==Resolución del problema completo==
Vamos a resolver el problema completo para ver como evolucionarían las poblaciones en un periodo de 40 días.
Primero supongamos una población infectada inicial de 20 individuos y una población en riesgo de contagio de 800.
Después supondremos el caso de 40 individuos enfermos y 10000 en riesgo de contagio.
Para resolverlo usaremos el método de Euler y más tarde el método de Runge-Kutta de orden cuatro con distintas discretizaciones para ver la influencia de estas.

<math>\Rightarrow</math> Método de Euler
{{matlab|codigo=

%Euler
clear all
%Datos
%Introducimos los valores de las constantes
a=0.003;
b=0.3;
c=0.01;
h=input('Introducir valores de h:');
t0=0;
y0=input('Introducir vector [So,Io]:');
tN=40;
% calculamos los subintervalos
t=t0:h:tN;
N=(tN-t0)/h;
y=zeros(2,length(t));
y(:,1)=y0';
for i= 1:N;
y(:,i+1)=y(:,i)+h*[-a*y(1,i).*y(2,i);a*y(1,i).*y(2,i)-b*y(2,i)-c*y(2,i)];
end

% grafico
plot(t,y)
% vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h
legend('Población sana','Location','best','Población enferma','Location','best');

}}

<math>\Rightarrow</math> Método de Runge-Kutta

{{matlab|codigo=

%Runge-Kutta
clear all

%Datos
S0=input('Introducir valor inicial SO: ');
I0=input('Introducir valor inicial IO: ');
a=0.003;
b=0.3;
c=0.01;
%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector solucion
y0=[S0;I0];
y=zeros(2,N+1);
y(:,1)=y0;
for i=1:N
K1S=-a*y(1,i)*y(2,i);
K2S=-a*(y(1,i)+(1/2)*K1S*h)*(y(2,i)+(1/2)*K1S*h);
K3S=-a*(y(1,i)+(1/2)*K2S*h)*(y(2,i)+(1/2)*K2S*h);
K4S=-a*(y(1,i)+K3S*h)*(y(2,i)+K3S*h);
y(1,i+1)=y(1,i)+(h/6)*(K1S+2*K2S+2*K3S+K4S);
K1I=a*y(1,i)*y(2,i)-b*y(2,i)-c*y(2,i);
K2I=a*(y(1,i)+(1/2)*K1I*h)*(y(2,i)+(1/2)*K1I*h)-b*(y(2,i)+(1/2)*K1I*h)-c*(y(2,i)+(1/2)*K1I*h);
K3I=a*(y(1,i)+(1/2)*K2I*h)*(y(2,i)+(1/2)*K2I*h)-b*(y(2,i)+(1/2)*K2I*h)-c*(y(2,i)+(1/2)*K2I*h);
K4I=a*(y(1,i)+K3I*h)*(y(2,i)+K3I*h)-b*(y(2,i)+K3I*h)-c*(y(2,i)+K3I*h);
y(2,i+1)=y(2,i)+(h/6)*(K1I+2*K2I+2*K3I+K4I);
end

%Tabla de resutados
[t',y(1,:)',y(2,:)']
%Gráfico
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'k')
legend('Poblacion sana','Poblacion enferma','Location','best');
hold off

% vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h

}}

La dificultad que hay en el uso del método del trapecio o cualquier otro método implícito está en que las incógnitas de nuestro sistema <math>(S,I)</math> quedan implícitas en ambas ecuaciones y necesitan que se despejen. Al comenzar a realizar cálculos nos damos cuenta que el despeje es imposible, ya que las dos variables incógnita aparecen multiplicándose entre sí, quedando una variable dependiendo de la otra siempre.

===Situación inicial de <math>(S _{0} , I_{0} )=(800,20)</math>===
En estos programas, sustituyendo las <math>h</math> para 0.1, 0.01, 0.001 y 0.0001 obtendremos los siguientes resultados:

-Para <math>h=0.1</math>

[[Archivo:Tborraré.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.1rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

Como se puede ver la población sana se infecta rápidamente y a partir del día 15 prácticamente toda la población ha fallecido.
* Método de Euler: nuestro programa nos indica que en este caso el momento en el que los enfermos son máximos se produce al tercer día (Día 2.8) y el número de enfermos en ese momento es 517.
* Método de Runge-Kutta: el momento en el que los enfermos son máximos se produce al tercer día (Día 2.6) y el número de enfermos en ese momento es 650.

-Para <math>h=0.01</math>
[[Archivo:Tborraré2.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.01rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

En el gráfico apenas se encuentra diferencia, lo cual es lógico porque es el mismo problema, pero si sobrepusiéramos las gráficas, lograríamos ver como se va ajustando a la realidad en función del aumento de h.
* Método de Euler: en este caso el momento en el que los enfermos son máximos es también en el tercer día (2.7) y el número de enfermos es 506.
* Método de Runge-Kutta: el momento en el que los enfermos son máximos se produce al tercer día (Día 2.69) y el número de enfermos en ese momento es 517.

-Para <math>h=0.001</math>
[[Archivo:Tborraré3.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.001rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* Método de Euler: en este caso el momento en el que los enfermos son máximos se produce el tercer día (2.696) con un valor de 505 enfermos.
* Método de Runge-Kutta: el momento en el que los enfermos son máximos se produce el tercer día (2.694) con un valor de 506 enfermos.

-Para <math>h=0.0001</math>
[[Archivo:Tborraré4.png|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:H0.0001rkaa.png|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* Método de Euler: al igual que antes, se repite el momento en el que los enfermos son máximos en el tercer día (2.695) con un valor de 505 enfermos.
* Método de Runge-Kutta: se repite el momento en el que los enfermos son máximos en el tercer día (2.6948) con un valor de 505 enfermos.

===Situación inicial de <math>(S _{0} , I_{0} )=(10000,40)</math>===
En estos programas, sustituyendo las <math>h</math> para 0.1, 0.01, 0.001 y 0.0001 obtendremos los siguientes resultados:

-Para <math>h=0.1</math>

[[Archivo:|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* Método de Euler:0.5 1.26e4
* Método de Runge-Kutta: 0.4 2.596e74

-Para <math>h=0.01</math>
[[Archivo:|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* Método de Euler:0.35 9.5211e3
* Método de Runge-Kutta: 0.3100 1.3546e4

-Para <math>h=0.001</math>
[[Archivo:|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* Método de Euler: 0.3420 9.4696e3
* Método de Runge-Kutta: 0.3380 9.5781e3
9
-Para <math>h=0.0001</math>
[[Archivo:|400px|miniaturadeimagen|left|Método de Euler]] [[Archivo:|400px|miniaturadeimagen|centro|Método de Runge-Kutta]]

* Método de Euler: 0.3411 9.4647e3
* Método de Runge-Kutta: 0.347 9.4298e3

==Resolución del problema completo suponiendo <math> a </math> una función variable con el tiempo==
Para el siguiente caso supondremos que el coeficiente <math> a </math> (tasa de contagio en población sana) es una función dada por :

<math> a(t)= \frac{0.003}{(1+t)} </math>
Suponiendo además las condiciones iniciales <math> (S_{0},I_{0})=(1640,40) </math> , vamos a dibujar las gráficas de infectados y sanos por el método de Heun, tomando un tamaño de paso <math> h=0.01 </math>
[[Archivo: H0.01heunn.png|1000px|thumb|rigth| Población infectada-Población sana para a=a(t)]]
{{matlab|codigo=
%Heun. Coeficiente variable a(t)

%Datos
S0=input('Introducir valor inicial SO: ');
I0=input('Introducir valor inicial IO: ');
b=0.3;
c=0.01;
%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector solucion
y0=[S0;I0];
y=zeros(2,N+1);
y(:,1)=y0;
for i=1:N
K1S=-(0.003/(1+t(i)))*y(1,i)*y(2,i);
K2S=-(0.003/(1+t(i)+h))*(y(1,i)+K1S*h)*(y(2,i)+K1S*h);
y(1,i+1)=y(1,i)+(h/2)*(K1S+K2S);
K1I=(0.003/(1+t(i)))*y(1,i)*y(2,i)-b*y(2,i)-c*y(2,i);
K2I=(0.003/(1+t(i)+h))*(y(1,i)+K1I*h)*(y(2,i)+K1I*h)-b*(y(2,i)+K1I*h)-c*(y(2,i)+K1I*h);
y(2,i+1)=y(2,i)+(h/2)*(K1I+K2I);
end

%Tabla de resutados
[t',y(1,:)',y(2,:)']
%Graficos
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'k')
title('Método de Heun');
legend('Poblacion sana','Población infectada','Location','best');
hold off

% Vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h

}}

Como se puede ver la población sana se infecta rápidamente y a partir del día 15 casi toda la población ha fallecido. Al segundo día <math>(t=2,15)</math> observamos que el número de infectados es de 994 individuos <math>(I=994,3)</math>, siendo este el máximo que se alcanzará.

A medida que va pasando el tiempo el coeficiente <math> a(t) </math> (tasa de contagio en población sana) va a ir decreciendo poco a poco hasta aproximarse a cero según podemos ver en el siguiente gráfico:
<math> \lim_{t\to \infty } a(t)=0</math>
[[Archivo: A(t).png|450px|thumb|right|a(t) para t=[0,40]]]
{{matlab|codigo=

%tasa de contagio en poblacion sana a(t)

%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector a
a=zeros(1,N+1);
for i=1:N+1
a(i)=0.003/(1+t(i));
end

%Grafica
plot(t,a)

}}

[[Archivo: Comparacionheun.png|border|1500px|center|frame|Poblaciones según la tasa de contagio a(t)]]

Ambas poblaciones parten del mismo inicio (para <math>t=0</math> las poblaciones son iguales). Podemos observar que la población sana (susceptible de contraer la enfermedad) en el caso en el que la tasa de contagio en población sana <math>(a)</math> es constante <math>(a=0.003)</math> se ve reducida antes que en el caso en el que depende del tiempo <math>(a=a(t))</math>.

<math>a(t)</math> disminuye a la vez que transcurre el tiempo, por lo que en general las pendientes tanto de la población sana como de la infectada serán menores para este caso que para el caso en el que <math>a=0.003</math>

==Calibración del parámetro <math>a</math> en base a una experiencia previa==

Modelos epidemiológicos (Grupo 3A)

2016-05-01T15:23:55Z

Trabajocampos: /* Resolución del problema con una sola ecuación diferencial */

Modelos epidemiológicos (Grupo 3A)

2016-05-01T15:01:54Z

Trabajocampos: /* Resolución del problema con una sola ecuación diferencial */

Archivo:Ejerciciotrestrapecio.png

2016-05-01T15:00:41Z

Trabajocampos:

Modelos epidemiológicos (Grupo 3A)

2016-05-01T14:42:47Z

Trabajocampos: /* Resolución del problema completo */

Modelos epidemiológicos (Grupo 3A)

2016-05-01T14:38:16Z

Trabajocampos:

Modelos epidemiológicos (Grupo 3A)

2016-05-01T13:03:53Z

Trabajocampos: /* Resolución del problema con una sola ecuación diferencial */

Modelos epidemiológicos (Grupo 3A)

2016-05-01T13:02:37Z

Trabajocampos: /* Resolución del problema con una sola ecuación diferencial */

Archivo:Ejerciciotreseuler.png

2016-05-01T13:02:28Z

Trabajocampos:

Modelos epidemiológicos (Grupo 3A)

2016-05-01T12:57:55Z

Trabajocampos: /* Resolución del problema con una sola ecuación diferencial */

Modelos epidemiológicos (Grupo 3A)

2016-05-01T12:54:16Z

Trabajocampos: /* Resolución del problema con una sola ecuación diferencial */

Modelos epidemiológicos (Grupo 3A)

2016-05-01T12:43:02Z

Trabajocampos: /* Resolución del problema con una sola ecuación diferencial */

Modelos epidemiológicos (Grupo 3A)

2016-05-01T12:26:07Z

Trabajocampos: /* Resolución del problema con una sola ecuación diferencial */

Archivo:Ejerciciodostrapecio.png

2016-05-01T12:23:19Z

Trabajocampos:

Modelos epidemiológicos (Grupo 3A)

2016-05-01T12:19:16Z

Trabajocampos: /* Resolución del problema con una sola ecuación diferencial */

Modelos epidemiológicos (Grupo 3A)

2016-05-01T11:45:49Z

Trabajocampos: /* Resolución del problema con una sola ecuación diferencial */

Archivo:Ejerciciodoseuler.png

2016-05-01T11:43:58Z

Trabajocampos:

Modelos epidemiológicos (Grupo 3A)

2016-05-01T11:39:21Z

Trabajocampos: /* Resolución del problema con una sola ecuación diferencial */

Modelos epidemiológicos (Grupo 3A)

2016-05-01T11:38:26Z

Trabajocampos:

Modelos epidemiológicos (Grupo 3A)

2016-04-30T17:59:39Z

Trabajocampos: /* Resolución del problema con una sola ecuación diferencial */

Modelos epidemiológicos (Grupo 3A)

2016-04-30T17:59:22Z

Trabajocampos: /* Resolución del problema con una sola ecuación diferencial */