Diferencia entre revisiones de «Modelos Epidemiológicos Grupo 3A»
(→Resolución del modelo completo mediante Euler) |
(→Resolución del PVI mediante los métodos de Euler y Trapecio) |
||
| (No se muestran 42 ediciones intermedias de 4 usuarios) | |||
| Línea 1: | Línea 1: | ||
| + | {{ TrabajoED | Modelos Epidemiológicos Grupo 3A | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]| | ||
| + | |||
| + | [[:Categoría:ED14/15|2014-15]] | Daniel Bello Bárcenas, Luis de la Fuente Puig, Ignacio Cobo Muñoz, Germán Sampedro Feito, Mercedes Cid Ruiz , Alfonso Martín de Soto Aláez }} | ||
| + | [[Categoría:Ecuaciones Diferentiales]] | ||
| + | [[Categoría:ED14/15]] | ||
| + | [[Categoría:Trabajos 2014-15]] | ||
| + | |||
| + | |||
=== Introducción e hipótesis iniciales === | === Introducción e hipótesis iniciales === | ||
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 temporal. | 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 temporal. | ||
| Línea 17: | Línea 25: | ||
donde a,b y c son parámetros. Los interpretamos como: | donde a,b y c son parámetros. Los interpretamos como: | ||
| − | + | ||
| + | (dI/dt) es la variación de la población infectada y (dS/dt) la variación de la población susceptible . | ||
| + | La población infectada se ve alterada por: | ||
* aSI ( con signo positivo), donde SI es la interacción entre ambas poblaciones, que depende de a | * aSI ( con signo positivo), donde SI es la interacción entre ambas poblaciones, que depende de a | ||
* -bI (con signo negativo) que son los que fallecen. Por lo tanto, 'b' es la '''tasa de fallecimiento'''. | * -bI (con signo negativo) que son los que fallecen. Por lo tanto, 'b' es la '''tasa de fallecimiento'''. | ||
| Línea 31: | Línea 41: | ||
Antes de pasar a la resolución completa del sistema, lo simplificaremos asignando el valor "cero" a la variable que representa a los '''sujetos susceptibles de contraer la enfermedad''' ('''S=0'''). Con ello, reducimos el sistema a una única ecuación: | Antes de pasar a la resolución completa del sistema, lo simplificaremos asignando el valor "cero" a la variable que representa a los '''sujetos susceptibles de contraer la enfermedad''' ('''S=0'''). Con ello, reducimos el sistema a una única ecuación: | ||
| − | \frac{dI}{dt}=-bI-cI | + | |
| + | <math> \frac{dI}{dt}=-bI-cI</math> | ||
Mediante sentencias Matlab y la asignación a la variable "'''I'''" ('''sujetos infectados'''), y a los coeficientes "'''b'''" y "'''c'''" de los valores deseados, resolvemos la ecuación diferencial simplificada aplicando los métodos de '''Euler''' y '''Trapecio''' (para ello se ha utilizado, en ambos métodos, una espaciación de '''0.1 unidades'''). Procedemos, en primer lugar, a la aproximación de nuestra función por el método de Euler: | Mediante sentencias Matlab y la asignación a la variable "'''I'''" ('''sujetos infectados'''), y a los coeficientes "'''b'''" y "'''c'''" de los valores deseados, resolvemos la ecuación diferencial simplificada aplicando los métodos de '''Euler''' y '''Trapecio''' (para ello se ha utilizado, en ambos métodos, una espaciación de '''0.1 unidades'''). Procedemos, en primer lugar, a la aproximación de nuestra función por el método de Euler: | ||
| Línea 101: | Línea 112: | ||
Es observable que en ambas gráficas se ha empleado el mismo valor inicial de sujetos infectados ('''I'''=2000) y que dicha representación de la gráfica finaliza en un valor para el tiempo de '''t'''=4.5, el cual se identifica con una reducción de la población inicial de sujetos infectados a '''la cuarta parte''' ('''I'''=500), ya sea por defunción o por cura de la patología (el carácter descendente de la función quedó precisado mediante la definición de los coeficientes "b" y "c" al inicio del artículo). | Es observable que en ambas gráficas se ha empleado el mismo valor inicial de sujetos infectados ('''I'''=2000) y que dicha representación de la gráfica finaliza en un valor para el tiempo de '''t'''=4.5, el cual se identifica con una reducción de la población inicial de sujetos infectados a '''la cuarta parte''' ('''I'''=500), ya sea por defunción o por cura de la patología (el carácter descendente de la función quedó precisado mediante la definición de los coeficientes "b" y "c" al inicio del artículo). | ||
| − | + | De forma análoga a la situación anterior, vamos a proceder a tomar como valor "'''S'''" de los '''sujetos susceptibles de ser infectados''' un valor constante de "cien" ('''S=100'''). La resolución, al igual que en caso anterior, se ejecutará mediante los métodos de Euler y del Trapecio respectivamente: | |
[[File:Euler3.jpg|500px|thumb|rigth|Aproximación del sistema por el método de Euler]] | [[File:Euler3.jpg|500px|thumb|rigth|Aproximación del sistema por el método de Euler]] | ||
| Línea 174: | Línea 185: | ||
}} | }} | ||
| − | + | No obstante, el sistema anteriormente resuelto, "esconde" una peculiaridad que quizás no resulta tan trivial como las cuestiones anteriormente expuestas: la asignación a la variable "'''S'''" (sujetos susceptibles de contraer la enfermedad) del valor "'''200'''". | |
| − | + | Si asignamos dicho valor a la sentencia MATLAB "input" del programa anterior (independientemente de qué método se haya aplicado), este último se colapsará. La respuesta por tanto, a esta situación que nos plantea "'''S=200'''", reside en que es a partir de este valor cuando la función '''abandona su carácter decreciente''' y comienza a crecer de forma exponencial tal y como marcaría la solución trivial de la segunda ecuación del sistema (una exponencial). De esta forma, el valor de "'''S=200'''", en combinación con la ya definida '''tasa de contagio''' "'''a'''", hace que la enfermedad sea imparable hasta el punto de que el número de infectados supere al de defunciones y curas; '''toda la población termina por infectarse'''. | |
| − | + | ||
| − | + | ||
| − | + | ||
=== Resolución del modelo completo mediante Euler === | === Resolución del modelo completo mediante Euler === | ||
| Línea 222: | Línea 230: | ||
}} | }} | ||
| − | Como hemos dicho | + | [[File:figura1trabajo4.jpg|border|550px|left|Aproximación para los valores: h=0.1 y (S<sub>0</sub>, I<sub>0</sub>)=(800, 20)]] |
| + | [[File:figura2trabajo4.jpg|border|550px|right|Aproximación para los valores: h=0.1 y (S<sub>0</sub>, I<sub>0</sub>)=(10000, 40)]] | ||
| + | (izquierda: ''h=0.1 y (S<sub>0</sub>, I<sub>0</sub>)=(800, 20)'', derecha:''h=0.1 y (S<sub>0</sub>, I<sub>0</sub>)=(10000, 40)'') | ||
| + | [[File:figura3trabajo4.jpg|border|550px|left|Aproximación para los valores: h=0.01 y (S<sub>0</sub>, I<sub>0</sub>)=(800, 20)]] | ||
| + | [[File:figura4trabajo4.jpg|border|550px|right|Aproximación para los valores: h=0.01 y (S<sub>0</sub>, I<sub>0</sub>)=(1000, 40)]] | ||
| + | (izquierda: ''h=0.01 y (S<sub>0</sub>, I<sub>0</sub>)=(800, 20)'', derecha:''h=0.01 y (S<sub>0</sub>, I<sub>0</sub>)=(10000, 40)'') | ||
| + | [[File:figura5trabajo4.jpg|border|550px|left|Aproximación para los valores: h=0.001 y (S<sub>0</sub>, I<sub>0</sub>)=(800, 20)]] | ||
| + | [[File:figura6trabajo4.jpg|border|550px|right|Aproximación para los valores: h=0.001 y (S<sub>0</sub>, I<sub>0</sub>)=(1000, 40)]] | ||
| + | (izquierda: ''h=0.001 y (S<sub>0</sub>, I<sub>0</sub>)=(800, 20)'', derecha:''h=0.001 y (S<sub>0</sub>, I<sub>0</sub>)=(10000, 40)'') | ||
| + | [[File:figura7trabajo4.jpg|border|550px|left|Aproximación para los valores: h=0.0001 y (S<sub>0</sub>, I<sub>0</sub>)=(800, 20)]] | ||
| + | [[File:figura8trabajo4.jpg|border|550px|right|Aproximación para los valores: h=0.0001 y (S<sub>0</sub>, I<sub>0</sub>)=(1000, 40)]] | ||
| + | (izquierda: ''h=0.0001 y (S<sub>0</sub>, I<sub>0</sub>)=(800, 20)'', derecha:''h=0.0001 y (S<sub>0</sub>, I<sub>0</sub>)=(10000, 40)'') | ||
| + | |||
| + | Como hemos dicho, aquí están las ocho gráficas referentes a los distintos valores. Podemos observar que en las gráficas que tienen como valores iniciales (S<sub>0</sub>, I<sub>0</sub>)=(1000, 40) no apreciamos cómo descienden y ascienden los susceptibles y los infectados. Todo los contrario pasa con las gráficas que tienen como valores iniciales (S<sub>0</sub>, I<sub>0</sub>)=(800, 20). En ellas sí que podemos observar que cuando aumenta el número de infectados, disminuye el número de susceptibles. Otra observación que hacemos, es que el paso de discretización personal no nos afecta a la hora de visualización de las gráficas, todas tienen el mismo aspecto, cada uno referido a sus valores iniciales. | ||
Por lo tanto, a la hora de analizar las gráficas, lo podemos hacer metiendo en un grupo a las gráficas que tengan como valores iniciales (S<sub>0</sub>, I<sub>0</sub>)=(800, 20) y un segundo grupo donde los valores iniciales son (S<sub>0</sub>, I<sub>0</sub>)=(1000, 40). | Por lo tanto, a la hora de analizar las gráficas, lo podemos hacer metiendo en un grupo a las gráficas que tengan como valores iniciales (S<sub>0</sub>, I<sub>0</sub>)=(800, 20) y un segundo grupo donde los valores iniciales son (S<sub>0</sub>, I<sub>0</sub>)=(1000, 40). | ||
| Línea 229: | Línea 250: | ||
En las segundas podemos decir que la curva de los infectados es descendente, pero no podemos decir nada acerca de la curva de los susceptibles, ya que no la apreciamos. En ese segundo grupo no podemos hablar de los máximos y mínimos que tienen las curvas, solo que ambas tienden a cero. Podemos destacar la primera gráfica con estos datos iniciales, donde podemos observar como en los primeros cinco días aproximadamente, las dos curvas suben y bajan sin constancia. Observamos en este con gran claridad cual es el punto mínimo de los susceptibles y el máximo de los infectados; pero por el contrario, no podemos ver cual es el máximo de los susceptibles y el mínimo de los infectados ya que una curva se monta en la otra. Como en las otras gráficas, al final, las dos curvas tienden a cero. Una importante apreciación de esta gráfica es que nos salen personas negativas, lo que es físicamente imposible. Esto nos indica que el método de Euler no es un buen método para la aproximación de estas ecuaciones con los valores iniciales (S<sub>0</sub>, I<sub>0</sub>)=(1000, 40) y h=0.1. | En las segundas podemos decir que la curva de los infectados es descendente, pero no podemos decir nada acerca de la curva de los susceptibles, ya que no la apreciamos. En ese segundo grupo no podemos hablar de los máximos y mínimos que tienen las curvas, solo que ambas tienden a cero. Podemos destacar la primera gráfica con estos datos iniciales, donde podemos observar como en los primeros cinco días aproximadamente, las dos curvas suben y bajan sin constancia. Observamos en este con gran claridad cual es el punto mínimo de los susceptibles y el máximo de los infectados; pero por el contrario, no podemos ver cual es el máximo de los susceptibles y el mínimo de los infectados ya que una curva se monta en la otra. Como en las otras gráficas, al final, las dos curvas tienden a cero. Una importante apreciación de esta gráfica es que nos salen personas negativas, lo que es físicamente imposible. Esto nos indica que el método de Euler no es un buen método para la aproximación de estas ecuaciones con los valores iniciales (S<sub>0</sub>, I<sub>0</sub>)=(1000, 40) y h=0.1. | ||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
=== Resolución del modelo completo mediante Runge-Kutta de orden 4 === | === Resolución del modelo completo mediante Runge-Kutta de orden 4 === | ||
| Línea 292: | Línea 303: | ||
A continuación, comparamos las gráficas obtenidas con (S<sub>0</sub>, I<sub>0</sub>)=(800, 20) entre Runge-Kutta y Euler: | A continuación, comparamos las gráficas obtenidas con (S<sub>0</sub>, I<sub>0</sub>)=(800, 20) entre Runge-Kutta y Euler: | ||
| − | [[File:Rungekutta3A15.jpg| | + | [[File:Rungekutta3A15.jpg|border|550px|left|Aproximación del sistema mediante RK-4]] |
| − | [[File:Euler3A15.jpg| | + | [[File:Euler3A15.jpg|border|550px|right|Aproximación del sistema mediante Euler]] |
| + | (izquierda:''Aproximación mediante Runge-Kutta'',derecha:''Aproximación mediante Euler'') | ||
| − | Ambas gráficas son prácticamente idénticas, se observa que el total de personas susceptibles baja drásticamente en los primeros 5 días hasta ser prácticamente nulo, consecuentemente se observa que el número de infectados aumenta rápidamente en ese periodo hasta alcanzar su máximo(650 infectados aprox.)para luego descender. Finalmente, cuando el tiempo es de 20 días ambas poblaciones son muy próximas a cero. | + | Ambas gráficas son prácticamente idénticas(varía algo el máximo de infectados), se observa que el total de personas susceptibles baja drásticamente en los primeros 5 días hasta ser prácticamente nulo, consecuentemente se observa que el número de infectados aumenta rápidamente en ese periodo hasta alcanzar su máximo(650 infectados aprox.)para luego descender. Finalmente, cuando el tiempo es de 20 días ambas poblaciones son muy próximas a cero. |
Al trabajar con los valores (S<sub>0</sub>, I<sub>0</sub>)=(10000, 40) las aproximaciones mediante Euler y RK4 muestran gráficas incongruentes con la evolución de personas susceptibles e infectadas observada hasta ahora, en las aproximaciones anteriores hemos obtenido gráficas de tipo exponencial con un máximo de infectados y un descenso gradual de ambas poblaciones cuando el tiempo toma valores altos, sin embargo, en este caso mediante Euler hemos obtenido gráficas con una sucesión de máximos y mínimos en intervalos de tiempo muy pequeños( en cuestión de 1 día las cifras de infectados oscilan entre pocos cientos y más de 10000 alternativamente) y además las personas susceptibles adquieren valores negativos(Físicamente imposible) lo cual nos indica que '''el método de Euler es inadecuado''' para la aproximación requerida, seguramente debido a la inestabilidad que presenta en la resolución de sistemas.Con el método de Runge-Kutta se obtiene de la misma manera un resultado desconcertante con respecto a las gráficas expuestas anteriormente, ofreciendo aproximaciones que no son equiparables a una evolución epidemiológica real. Otro posible motivo del error en las gráficas es la incompatibilidad de la pareja de valores iniciales con la aproximación de una epidemia expuesta en este trabajo, es decir, para valores muy altos de S(t) la modelización de una epidemia no puede ser aproximada por un sistema de ecuaciones diferenciales como el proporcionado para este trabajo(Obtenemos valores que no concuerdan con la realidad). | Al trabajar con los valores (S<sub>0</sub>, I<sub>0</sub>)=(10000, 40) las aproximaciones mediante Euler y RK4 muestran gráficas incongruentes con la evolución de personas susceptibles e infectadas observada hasta ahora, en las aproximaciones anteriores hemos obtenido gráficas de tipo exponencial con un máximo de infectados y un descenso gradual de ambas poblaciones cuando el tiempo toma valores altos, sin embargo, en este caso mediante Euler hemos obtenido gráficas con una sucesión de máximos y mínimos en intervalos de tiempo muy pequeños( en cuestión de 1 día las cifras de infectados oscilan entre pocos cientos y más de 10000 alternativamente) y además las personas susceptibles adquieren valores negativos(Físicamente imposible) lo cual nos indica que '''el método de Euler es inadecuado''' para la aproximación requerida, seguramente debido a la inestabilidad que presenta en la resolución de sistemas.Con el método de Runge-Kutta se obtiene de la misma manera un resultado desconcertante con respecto a las gráficas expuestas anteriormente, ofreciendo aproximaciones que no son equiparables a una evolución epidemiológica real. Otro posible motivo del error en las gráficas es la incompatibilidad de la pareja de valores iniciales con la aproximación de una epidemia expuesta en este trabajo, es decir, para valores muy altos de S(t) la modelización de una epidemia no puede ser aproximada por un sistema de ecuaciones diferenciales como el proporcionado para este trabajo(Obtenemos valores que no concuerdan con la realidad). | ||
| Línea 301: | Línea 313: | ||
A modo de aclaración(y curiosidad) se adjuntan dichas gráficas: | A modo de aclaración(y curiosidad) se adjuntan dichas gráficas: | ||
| − | [[File:RK3A15.jpg| | + | [[File:RK3A15.jpg|border|550px|left|Aproximación del sistema mediante RK-4]] |
| − | [[File:Figura2trabajo4.jpg| | + | [[File:Figura2trabajo4.jpg|border|550px|right|Aproximación del sistema mediante Euler]] |
| − | + | (izquierda:''Aproximación mediante Runge-Kutta'',derecha:''Aproximación mediante Euler'') | |
| Línea 371: | Línea 383: | ||
| − | |||
| − | [[Archivo:Figurafunciona.jpg| | + | |
| + | [[Archivo:Figurafunciona.jpg|400px|frame|Visualización de a(t) en t=[0,40]]] | ||
a(t) se trata de una función no definida en t=-1. Esto no nos importa, ya que t es el tiempo, que siempre es positivo y concretamente toma valores entre 0 y 40. Como la t>=0 siempre: | a(t) se trata de una función no definida en t=-1. Esto no nos importa, ya que t es el tiempo, que siempre es positivo y concretamente toma valores entre 0 y 40. Como la t>=0 siempre: | ||
| Línea 384: | Línea 396: | ||
| − | En la gráfica | + | En la gráfica inferior se muestra cómo varia la misma población ( 40 infectados en 1640 individuos) con diferentes coeficientes a(t): |
* Uno constante e igual a 0.003 como en apartados anteriores. | * Uno constante e igual a 0.003 como en apartados anteriores. | ||
| Línea 390: | Línea 402: | ||
* Otro variable variable que es <math> a(t)=\frac{0.003}{1+t}</math> | * Otro variable variable que es <math> a(t)=\frac{0.003}{1+t}</math> | ||
| − | [[File:Figuracomparativa5.jpg|border| | + | [[File:Figuracomparativa5.jpg|border|400px|center|frame|Comparación de poblaciones según la tasa de contagio 'a']] |
| Línea 400: | Línea 412: | ||
=== Calibración del parámetro 'a' en base a una experiencia previa=== | === Calibración del parámetro 'a' en base a una experiencia previa=== | ||
| − | En este apartado cogemos los mismos datos del apartado anterior, con la diferencia que el parámetro a deja de ser una función y se convierte en un vector, de donde iremos cogiendo valores equidistribuidos hasta agotarlos. El objetivo final es encontrar de todos los valores que ha ido cogiendo el parámetro a, el que más infectados tenga estando | + | En este apartado cogemos los mismos datos del apartado anterior, con la diferencia que el parámetro a deja de ser una función y se convierte en un vector comprendido entre 0.0005 y 0.002, de donde iremos cogiendo valores equidistribuidos hasta agotarlos. El objetivo final es encontrar de todos los valores que ha ido cogiendo el parámetro a, el que más infectados tenga estando lo más cerca del quinto día. La aproximación la hemos hecho de nuevo mediante el '''método de Heun'''. |
{{matlab|codigo= | {{matlab|codigo= | ||
| Línea 413: | Línea 425: | ||
b=0.3; | b=0.3; | ||
c=0.01; | c=0.01; | ||
| − | h=0.1;% | + | h=0.1;%Tamaño de paso para t |
%Intervalo de tiempos | %Intervalo de tiempos | ||
t0=0; | t0=0; | ||
| Línea 421: | Línea 433: | ||
%Creo un vector t | %Creo un vector t | ||
t=t0:h:tN; | t=t0:h:tN; | ||
| − | % | + | %Preparo vector S e I para guardar las soluciones |
S=zeros(1,N+1);%vector incognita 1 | S=zeros(1,N+1);%vector incognita 1 | ||
I=zeros(1,N+1);%vector incognita 2 | I=zeros(1,N+1);%vector incognita 2 | ||
S(1)=S0; | S(1)=S0; | ||
I(1)=I0; | I(1)=I0; | ||
| − | % | + | %Preparo vectores m1(máximos de I) m2(tiempos de los máximos) |
m1=zeros(1,length(a)); | m1=zeros(1,length(a)); | ||
m2=zeros(1,length(a)); | m2=zeros(1,length(a)); | ||
| Línea 460: | Línea 472: | ||
}} | }} | ||
| − | En la siguiente gráfica se muestran los máximos valores de infectados y los tiempos correspondientes para los que se produce el máximo según los parámetros a equidistribuidos del vector inicial. El objetivo es encontrar el más próximo al tiempo 5. | + | En la siguiente gráfica se muestran los máximos valores de infectados y los tiempos correspondientes para los que se produce el máximo según los parámetros a equidistribuidos del vector inicial.El objetivo es encontrar el más próximo al tiempo 5. |
| − | [[File:Graficamatlab.jpg| | + | [[File:Graficamatlab.jpg|550px|left]] |
| + | |||
| + | |||
| + | [[File:disp3A15.jpg|550px|right]] | ||
| + | |||
| + | |||
| + | |||
| + | |||
| + | |||
| + | |||
| + | |||
| + | |||
| + | |||
| + | |||
| + | |||
| + | |||
| + | |||
| + | |||
| + | |||
| + | |||
| + | Observación: para el valor de a mas pequeño le corresponde el nº menor de infectados y que a medida que va incrementando el valor del parámetro, también crece el valor máximo de infectados, lo que es lógico. En la gráfica hemos representado los valores máximos para cada parámetro y son independientes uno de otro, digo esto porque podría llevarnos al error de pensar que para el valor más pequeño de a, hay un número mayor de infectados. En nuestro caso, hay un valor de 'a' para el cual el tiempo del máximo es exactamente 5, siendo 1192 el máximo de infectados como queda reflejado en el gráfico y en el recuadro. | ||
Revisión actual del 00:12 7 mar 2015
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Modelos Epidemiológicos Grupo 3A |
| Asignatura | Ecuaciones Diferenciales |
| Curso | |
| Autores | Daniel Bello Bárcenas, Luis de la Fuente Puig, Ignacio Cobo Muñoz, Germán Sampedro Feito, Mercedes Cid Ruiz , Alfonso Martín de Soto Aláez |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
- 1 Introducción e hipótesis iniciales
- 2 Resolución del PVI mediante los métodos de Euler y Trapecio
- 3 Resolución del modelo completo mediante Euler
- 4 Resolución del modelo completo mediante Runge-Kutta de orden 4
- 5 Resolución del modelo completo mediante Heun con a(t) una función dependiente del tiempo
- 6 Calibración del parámetro 'a' en base a una experiencia previa
1 Introducción e hipótesis iniciales
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 temporal. Se dan 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.
Mediante el siguiente sistema de ecuaciones diferenciales se muestran las variaciones de ambas poblaciones respecto al tiempo t:
\begin{matrix}
\frac{dS}{dt} = -aSI \\
\frac{dI}{dt} = aSI-bI-cI
\end{matrix}
donde a,b y c son parámetros. Los interpretamos como:
(dI/dt) es la variación de la población infectada y (dS/dt) la variación de la población susceptible .
La población infectada se ve alterada por:
- aSI ( con signo positivo), donde SI es la interacción entre ambas poblaciones, que depende de a
- -bI (con signo negativo) que son los que fallecen. Por lo tanto, 'b' es la tasa de fallecimiento.
- -cI (con signo negativo) que son los que se han curado. Es decir, 'c' es la tasa de curación.
Para interpretar 'a' observamos el sistema y la segunda hipótesis:
[math] \frac{dI}{dt} = aSI-bI-cI[/math]
Como los individuos que pasan de ser susceptibles a estar infectados son proporcionales a su interacción(SI), 'a' es la tasa de contagio en la población susceptible.
2 Resolución del PVI mediante los métodos de Euler y Trapecio
Antes de pasar a la resolución completa del sistema, lo simplificaremos asignando el valor "cero" a la variable que representa a los sujetos susceptibles de contraer la enfermedad (S=0). Con ello, reducimos el sistema a una única ecuación:
[math] \frac{dI}{dt}=-bI-cI[/math]
Mediante sentencias Matlab y la asignación a la variable "I" (sujetos infectados), y a los coeficientes "b" y "c" de los valores deseados, resolvemos la ecuación diferencial simplificada aplicando los métodos de Euler y Trapecio (para ello se ha utilizado, en ambos métodos, una espaciación de 0.1 unidades). Procedemos, en primer lugar, a la aproximación de nuestra función por el método de Euler:
%I(t): Población de individuos infectados
%S(t): Población susceptible de ser infectada
%los parámetros a=Tasa de contagio entre S e I,b=Tasa de
%fallecimiento de personas infectadas,c=Tasa de cura de personas infectadas
clear all clc
h=0.1;
format long
%Inicializo los vectores t,i e y asignando a la primera posición los valores
%y0 ,t0
t0=0;
y0=2000;
t(1)=t0;
y(1)=y0;
i=1;
%Defino las ctes b y c
b=0.3;
c=0.01;
%Bucle while: calculara valores de y(i) mientras y(i)>500
while y(i)>500
y(i+1)=y(i)+h*(-y(i)*(b+c));%Euler
t(i+1)=t(i)+h;
i=i+1;
end
[t',y']
disp('Tiempo final:')
disp(t(end));
plot(t,y);
legend('Población infectada(Euler)','Location','best');%Optimizar
Resolvemos de nuevo la ecuación simplificada, en este caso, mediante el método del trapecio:
%Método del Trapecio
%I(t): Población de individuos infectados
%S(t: Población susceptible de ser infectada
%los parámetros a=Tasa de contagio,b=Tasa de
%fallecimiento de personas infectadas,c=Tasa de cura de personas infectadas
clear all clc
format long
h=0.1;
%Inicializo los vectores t,i,z asignando a la primera posición los valores
%z0 ,t0
t0=0;
z0=2000;
t(1)=t0;
z(1)=z0;
i=1;
%Defino las ctes.b,c
b=0.3;
c=0.01;
%Bucle while: calculara valores de z(i) mientras z(i)>500
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 final:')
disp(t(end));
plot(t,z,'r');
legend('Población infectada(Trapecio)','Location','best');%OptimizarEs observable que en ambas gráficas se ha empleado el mismo valor inicial de sujetos infectados (I=2000) y que dicha representación de la gráfica finaliza en un valor para el tiempo de t=4.5, el cual se identifica con una reducción de la población inicial de sujetos infectados a la cuarta parte (I=500), ya sea por defunción o por cura de la patología (el carácter descendente de la función quedó precisado mediante la definición de los coeficientes "b" y "c" al inicio del artículo).
De forma análoga a la situación anterior, vamos a proceder a tomar como valor "S" de los sujetos susceptibles de ser infectados un valor constante de "cien" (S=100). La resolución, al igual que en caso anterior, se ejecutará mediante los métodos de Euler y del Trapecio respectivamente:
%Método de Euler
%I(t): Población de individuos infectados
%S(t): Población susceptible de ser infectada
%los parámetros a=Tasa de contagio,b=Tasa de
%fallecimiento de personas infectadas,c=Tasa de cura de personas infectadas
clear all clc
format long
h=0.1;
%Inicializo los vectores t,i e y asignando a la primera posición los valores
%y0 ,t0
i=1;
t0=0;
y0=2000;
t(1)=t0;
y(1)=y0;
%Defino las ctes. S,a,b,c
S=input('Introduce población susceptible de contagio:');
a=0.003;
b=0.3;
c=0.01;
%Bucle while: calculara valores de y(i) mientras y(i)>500
while y(i)>500
y(i+1)=y(i)+h*(-y(i)*(b+c-S*a));%Euler
t(i+1)=t(i)+h;
i=i+1;
end
[t',y']
disp('Tiempo final:')
disp(t(end));
plot(t,y);
Utilizamos ahora el método del trapecio:
%Método del Trapecio
%I(t): Población de individuos infectados
%S(t): Población susceptible de ser infectada
%los parámetros a=Tasa de contagio,b=Tasa de
%fallecimiento de personas infectadas,c=Tasa de cura de personas infectadas
clear all clc
format long
h=0.1;
%Inicializo los vectores t,i y z asignando a la primera posición los valores
%z0 ,t0
t0=0;
z0=2000;
t(1)=t0;
z(1)=z0;
i=1;
%Defino las ctes. S,a,b,c
S=input('Introduce población susceptible de contagio:');
a=0.003;
b=0.3;
c=0.01;
%Bucle while: calculara valores de z(i) mientras z(i)>500
while z(i)>500
z(i+1)=(z(i)*(1-(h/2)*(b+c-S*a)))/(1+h/2*(b+c-S*a));%trapecio
t(i+1)=t(i)+h;
i=i+1;
end
[t',z']
disp('Tiempo final:')
disp(t(end));
plot(t,z,'r');
legend('Población infectada(Trapecio)','Location','best');%Optimizar
No obstante, el sistema anteriormente resuelto, "esconde" una peculiaridad que quizás no resulta tan trivial como las cuestiones anteriormente expuestas: la asignación a la variable "S" (sujetos susceptibles de contraer la enfermedad) del valor "200".
Si asignamos dicho valor a la sentencia MATLAB "input" del programa anterior (independientemente de qué método se haya aplicado), este último se colapsará. La respuesta por tanto, a esta situación que nos plantea "S=200", reside en que es a partir de este valor cuando la función abandona su carácter decreciente y comienza a crecer de forma exponencial tal y como marcaría la solución trivial de la segunda ecuación del sistema (una exponencial). De esta forma, el valor de "S=200", en combinación con la ya definida tasa de contagio "a", hace que la enfermedad sea imparable hasta el punto de que el número de infectados supere al de defunciones y curas; toda la población termina por infectarse.
3 Resolución del modelo completo mediante Euler
En este apartado resolveremos el sistema completo y analizaremos las gráficas para diferentes valores de las variables. Los valores iniciales que tomamos son (S0, I0)=(800, 20) y (S0, I0)=(10000, 40); el tiempo estará en el intervalo [0,40] y como pasos de discretización temporal (h) 10-1, 10-2, 10-3, 10-4. El código de matlab utilizado en este apartado para resolverlo por el método de Euler:
%Metodo de Euler
% Considerando:
%I(t): Población de individuos infectados
%S(t): Población susceptible de ser infectada
%los parámetros a=Tasa de contagio,b=Tasa de
%fallecimiento de personas infectadas,c=Tasa de cura de personas infectadas
clear all clc
format long
h=input('Introduce tamaño de paso:');
S0=input('introduce un valor inicial S0:');
I0=input('introduce un valor inicial I0:');
a=0.003;
b=0.3;
c=0.01;
%Intervalo de tiempos
t0=0;
tN=40;
%Calculamos número de subintervalos
N=(tN-t0)/h;
%Creo un vector t
t=t0:h:tN;
%preparo vectores S,I para guardar las soluciones
S=zeros(1,N+1);%vector incognita 1
I=zeros(1,N+1);%vector incognita 2
S(1)=S0;
I(1)=I0;
for i=1:N
S(i+1)=S(i)+h*(-a*S(i)*I(i));%Euler para S
I(i+1)=I(i)+h*(a*S(i)*I(i)-b*I(i)-c*I(i));%Euler para I
end
[t',S',I']
hold on
plot (t,S)
plot(t,I,'r')
legend('Susceptibles','Infectados','Location','best');%Optimizar
hold off
(izquierda: h=0.1 y (S0, I0)=(800, 20), derecha:h=0.1 y (S0, I0)=(10000, 40))
(izquierda: h=0.01 y (S0, I0)=(800, 20), derecha:h=0.01 y (S0, I0)=(10000, 40))
(izquierda: h=0.001 y (S0, I0)=(800, 20), derecha:h=0.001 y (S0, I0)=(10000, 40))
(izquierda: h=0.0001 y (S0, I0)=(800, 20), derecha:h=0.0001 y (S0, I0)=(10000, 40))
Como hemos dicho, aquí están las ocho gráficas referentes a los distintos valores. Podemos observar que en las gráficas que tienen como valores iniciales (S0, I0)=(1000, 40) no apreciamos cómo descienden y ascienden los susceptibles y los infectados. Todo los contrario pasa con las gráficas que tienen como valores iniciales (S0, I0)=(800, 20). En ellas sí que podemos observar que cuando aumenta el número de infectados, disminuye el número de susceptibles. Otra observación que hacemos, es que el paso de discretización personal no nos afecta a la hora de visualización de las gráficas, todas tienen el mismo aspecto, cada uno referido a sus valores iniciales.
Por lo tanto, a la hora de analizar las gráficas, lo podemos hacer metiendo en un grupo a las gráficas que tengan como valores iniciales (S0, I0)=(800, 20) y un segundo grupo donde los valores iniciales son (S0, I0)=(1000, 40).
En las primeras vemos que la curva de los susceptibles es una curva descendente, la cual desciende con gran rapidez, ya que pasado el día cinco está muy próxima al cero. En cambio la curva de los infectados es ascendente hasta casi el día cinco, donde va a tener un máximo, y luego empieza a descender con gran rapidez. Al descender ambas curvas, se encontrarán aproximadamente sobre el día 20, donde las dos tienen una asíntota en y=0, y por lo tanto ambas tienden a cero.
En las segundas podemos decir que la curva de los infectados es descendente, pero no podemos decir nada acerca de la curva de los susceptibles, ya que no la apreciamos. En ese segundo grupo no podemos hablar de los máximos y mínimos que tienen las curvas, solo que ambas tienden a cero. Podemos destacar la primera gráfica con estos datos iniciales, donde podemos observar como en los primeros cinco días aproximadamente, las dos curvas suben y bajan sin constancia. Observamos en este con gran claridad cual es el punto mínimo de los susceptibles y el máximo de los infectados; pero por el contrario, no podemos ver cual es el máximo de los susceptibles y el mínimo de los infectados ya que una curva se monta en la otra. Como en las otras gráficas, al final, las dos curvas tienden a cero. Una importante apreciación de esta gráfica es que nos salen personas negativas, lo que es físicamente imposible. Esto nos indica que el método de Euler no es un buen método para la aproximación de estas ecuaciones con los valores iniciales (S0, I0)=(1000, 40) y h=0.1.
4 Resolución del modelo completo mediante Runge-Kutta de orden 4
En esta sección se pretende resolver de nuevo el mismo problema que en el apartado anterior, sin embargo, esta vez lo haremos con el método de Runge-Kutta de cuarto orden con el objetivo de compararlo con los resultados obtenidos mediante el método de Euler.El código matlab correspondiente queda expuesto a continuación:
%Trabajo Ecuaciones diferenciales
% Considerando:
%I(t): Población de individuos infectados
%S(t: Población susceptible de ser infectada
%los parámetros a=Tasa de contagio entre S e I,b=Tasa de
%fallecimiento de personas infectadas,c=Tasa de cura de personas infectadas
clear all clc
format long
S0=input('introduce un valor inicial S0:');
I0=input('introduce un valor inicial I0:');
h=input('introduce tamano de paso:');
a=0.003;
b=0.3;
c=0.01;
%Intervalo de tiempos
t0=0;
tN=40;
%Calculamos número de subintervalos
N=(tN-t0)/h;
%Creo un vector t
t=t0:h:tN;
%preparo vectores S,I para guardar las soluciones uso el comando zeros(filas,columnas)
S=zeros(1,N+1);%vector incognita 1
I=zeros(1,N+1);%vector incognita 2
S(1)=S0;
I(1)=I0;
for i=1:N
K1=-a*S(i)*I(i);
K2=-a*(S(i)+(1/2)*K1*h)*(I(i)+(1/2)*K1*h);
K3=-a*(S(i)+(1/2)*K2*h)*(I(i)+(1/2)*K2*h);
K4=-a*(S(i)+K3*h)*(I(i)+K3*h);
S(i+1)=S(i)+(h/6)*(K1+2*K2+2*K3+K4);%Runge-Kutta orden 4 para S
K5=a*S(i)*I(i)-b*I(i)-c*I(i);
K6=a*(S(i)+(1/2)*K5*h)*(I(i)+(1/2)*K5*h)-b*(I(i)+(1/2)*K5*h)-c*(I(i)+(1/2)*K5*h);
K7=a*(S(i)+(1/2)*K6*h)*(I(i)+(1/2)*K6*h)-b*(I(i)+(1/2)*K6*h)-c*(I(i)+(1/2)*K6*h);
K8=a*(S(i)+K7*h)*(I(i)+K7*h)-b*(I(i)+K7*h)-c*(I(i)+K7*h);
I(i+1)=I(i)+(h/6)*(K5+2*K6+2*K7+K8);%Runge-Kutta orden 4 para I
end
[t',S',I']
hold on
plot (t,S)
plot(t,I,'r')
legend('Susceptibles','Infectados','Location','best');%Optimizar
hold offA continuación, comparamos las gráficas obtenidas con (S0, I0)=(800, 20) entre Runge-Kutta y Euler:
(izquierda:Aproximación mediante Runge-Kutta,derecha:Aproximación mediante Euler)
Ambas gráficas son prácticamente idénticas(varía algo el máximo de infectados), se observa que el total de personas susceptibles baja drásticamente en los primeros 5 días hasta ser prácticamente nulo, consecuentemente se observa que el número de infectados aumenta rápidamente en ese periodo hasta alcanzar su máximo(650 infectados aprox.)para luego descender. Finalmente, cuando el tiempo es de 20 días ambas poblaciones son muy próximas a cero.
Al trabajar con los valores (S0, I0)=(10000, 40) las aproximaciones mediante Euler y RK4 muestran gráficas incongruentes con la evolución de personas susceptibles e infectadas observada hasta ahora, en las aproximaciones anteriores hemos obtenido gráficas de tipo exponencial con un máximo de infectados y un descenso gradual de ambas poblaciones cuando el tiempo toma valores altos, sin embargo, en este caso mediante Euler hemos obtenido gráficas con una sucesión de máximos y mínimos en intervalos de tiempo muy pequeños( en cuestión de 1 día las cifras de infectados oscilan entre pocos cientos y más de 10000 alternativamente) y además las personas susceptibles adquieren valores negativos(Físicamente imposible) lo cual nos indica que el método de Euler es inadecuado para la aproximación requerida, seguramente debido a la inestabilidad que presenta en la resolución de sistemas.Con el método de Runge-Kutta se obtiene de la misma manera un resultado desconcertante con respecto a las gráficas expuestas anteriormente, ofreciendo aproximaciones que no son equiparables a una evolución epidemiológica real. Otro posible motivo del error en las gráficas es la incompatibilidad de la pareja de valores iniciales con la aproximación de una epidemia expuesta en este trabajo, es decir, para valores muy altos de S(t) la modelización de una epidemia no puede ser aproximada por un sistema de ecuaciones diferenciales como el proporcionado para este trabajo(Obtenemos valores que no concuerdan con la realidad).
A modo de aclaración(y curiosidad) se adjuntan dichas gráficas:
(izquierda:Aproximación mediante Runge-Kutta,derecha:Aproximación mediante Euler)
Por otro lado, en métodos de aproximación como el trapezoidal o trapecio dado por la expresión:
\begin{matrix}
y_{0} & & & & \\
y_{n+1} = h/2 (f(t_{n},y_{n}) + f(t_{n+1},y_{n+1}))
\end{matrix}
La dificultad a la hora de resolver un sistema de ED como el arriba expuesto radica en que las variables incógnita de nuestro sistema (S e I) quedan implícitas en ambas ecuaciones y requieren despeje. Al comenzar a operar observaríamos que el despeje único es sumamente complicado o imposible, ya que las dos variables incógnita están acopladas simultáneamente en ambas ecuaciones y además multiplicadas entre sí y otros factores, es decir, al despejar completamente una de las dos variables, la otra queda incluida al otro lado de la igualdad y viceversa.
5 Resolución del modelo completo mediante Heun con a(t) una función dependiente del tiempo
Hasta ahora, el coeficiente a ( tasa de contagio en población susceptible) ha sido dada como un coeficiente constante. Ahora 'a' es una función dependiente del tiempo t (coeficiente variable) dada por:
[math] a(t)=\frac{0.003}{1+t}[/math]
%Trabajo Ecuaciones diferenciales
% Considerando:
%I(t): Población de individuos infectados
%S(t: Población susceptible de ser infectada
%los parámetros a=Tasa de contagio,b=Tasa de
%fallecimiento de personas infectadas,c=Tasa de cura de personas infectadas
clear all clc
format long
S0=1600;
I0=40;
h=input('introduce tamano de paso:');
b=0.3;
c=0.01;
%Intervalo de tiempos
t0=0;
tN=40;
%Calculamos número de subintervalos
N=(tN-t0)/h;
%Creo un vector t
t=t0:h:tN;
%preparo vector y para guardar las soluciones uso el comando zeros(filas,columnas)
S=zeros(1,N+1);%vector incognita 1
I=zeros(1,N+1);%vector incognita 2
S(1)=S0;
I(1)=I0;
for i=1:N
K1=-(0.003./(1+t(i)))*(S(i)*I(i));
K2=-(0.003./(1+t(i))).*(S(i)+K1*h)*(I(i)+K1*h);
K3=(0.003./(1+t(i))).*(S(i)*I(i)-b*I(i)-c*I(i));
K4=(0.003./(1+t(i))).*(S(i)+K3*h)*(I(i)+K3*h)-b*(I(i)+K3*h)-c*(I(i)+K3*h);
S(i+1)=S(i)+(h/2)*(K1+K2);
I(i+1)=I(i)+(h/2)*(K3+K4);
end
[t',S',I']
hold on
plot(t,S);
plot(t,I,'r');
legend('Susceptibles','Infectados','Location','best');%Optimizar
hold off
a(t) se trata de una función no definida en t=-1. Esto no nos importa, ya que t es el tiempo, que siempre es positivo y concretamente toma valores entre 0 y 40. Como la t>=0 siempre:
[math] \lim_{t\to \infty } a(t)=0[/math]
. Esto quiere decir que a medida que va pasando el tiempo, el coeficiente a(t) va a ir disminuyendo progresivamente hacia cero. Por tanto, la tasa de contagio disminuirá con el tiempo.
En la gráfica inferior se muestra cómo varia la misma población ( 40 infectados en 1640 individuos) con diferentes coeficientes a(t):
- Uno constante e igual a 0.003 como en apartados anteriores.
- Otro variable variable que es [math] a(t)=\frac{0.003}{1+t}[/math]
La función a(t) alcanza su máximo en el intervalo t= [0,40] en t=0. Coincidiendo en este punto con el coeficiente constante. Esto es, a=0.003. Esto significa que en t=0, ambas poblaciones son iguales.
Teniendo en cuenta que ambas poblaciones parten del mismo inicio, se aprecia claramente como la población de infectados que depende de un coeficiente variable a(t) es en todo momento en el intervalo t=[0,40] inferior a la población de infectados con un coeficiente a=0.003 constante a lo largo del tiempo.
Por otro lado, en lo que se refiere a la población de susceptibles, en la gráfica queda reflejado que cuando hablamos de a=0.003 constante a lo largo del tiempo, la población susceptible se ve reducida antes que cuando hablamos de a(t) variable. Esto tiene lógica porque si a(t) va disminuyendo a lo largo del tiempo, la tasa de contagio va disminuyendo, y menos cantidad de susceptibles se transforman en infectados en un determinado intervalo de tiempo . Por ejemplo, en t=2, la población de susceptibles con a=0.003 constante es prácticamente nula, lo que no ocurre con a variable donde aún queda cantidad de susceptibles apreciable ( no se contagia el mismo número de individuos porque la tasa de contagio es menor).
6 Calibración del parámetro 'a' en base a una experiencia previa
En este apartado cogemos los mismos datos del apartado anterior, con la diferencia que el parámetro a deja de ser una función y se convierte en un vector comprendido entre 0.0005 y 0.002, de donde iremos cogiendo valores equidistribuidos hasta agotarlos. El objetivo final es encontrar de todos los valores que ha ido cogiendo el parámetro a, el que más infectados tenga estando lo más cerca del quinto día. La aproximación la hemos hecho de nuevo mediante el método de Heun.
clear clc
format long
%Defino valores iniciales
S0=1600;
I0=40;
h1=0.0001;%Tamaño de paso para intervalo de a
a=0.0005:h1:0.002;
b=0.3;
c=0.01;
h=0.1;%Tamaño de paso para t
%Intervalo de tiempos
t0=0;
tN=40;
%Calculamos número de subintervalos
N=(tN-t0)/h;
%Creo un vector t
t=t0:h:tN;
%Preparo vector S e I para guardar las soluciones
S=zeros(1,N+1);%vector incognita 1
I=zeros(1,N+1);%vector incognita 2
S(1)=S0;
I(1)=I0;
%Preparo vectores m1(máximos de I) m2(tiempos de los máximos)
m1=zeros(1,length(a));
m2=zeros(1,length(a));
%Bucle for: calcula valores maximos de I para cada valor 'a' del intervalo y los almacena
%junto con su tiempo correspondiente en los vectores m1 y m2
%respectivamente.Metodo numérico de Heun
for j=1:length(a)
for i=1:N
K1=-(a(j))*(S(i)*I(i));
K2=-(a(j)).*(S(i)+K1*h)*(I(i)+K1*h);
K3=(a(j)).*(S(i)*I(i)-b*I(i)-c*I(i));
K4=(a(j)).*(S(i)+K3*h)*(I(i)+K3*h)-b*(I(i)+K3*h)-c*(I(i)+K3*h);
S(i+1)=S(i)+(h/2)*(K1+K2);
I(i+1)=I(i)+(h/2)*(K3+K4);
m1(j)=max(I);
m2(j)=t(max(I)==I);
end
%Bucle if: si la diferencia entre el valor de tiempo para el máximo y 5 es
%menor o igual que 0.1 muestra el valor de a, maximo de I y tiempo
%correspondiente para los que esto sucede
if abs(m2(j)-5)<=0.1
disp('Valor de a para el cual el máximo está mas próximo a 5:');
disp(a(j));
disp('Valor máximo de I correspondiente:');
disp((m1(j)));
disp('Valor de tiempo correspondiente:');
disp((m2(j)));
end
end
%Grafica con los máximos de I y los tiempos correspondientes
plot(m2,m1,'*');
legend('Máximos de infectados','Location','best');%OptimizarEn la siguiente gráfica se muestran los máximos valores de infectados y los tiempos correspondientes para los que se produce el máximo según los parámetros a equidistribuidos del vector inicial.El objetivo es encontrar el más próximo al tiempo 5.
Observación: para el valor de a mas pequeño le corresponde el nº menor de infectados y que a medida que va incrementando el valor del parámetro, también crece el valor máximo de infectados, lo que es lógico. En la gráfica hemos representado los valores máximos para cada parámetro y son independientes uno de otro, digo esto porque podría llevarnos al error de pensar que para el valor más pequeño de a, hay un número mayor de infectados. En nuestro caso, hay un valor de 'a' para el cual el tiempo del máximo es exactamente 5, siendo 1192 el máximo de infectados como queda reflejado en el gráfico y en el recuadro.
