Resolución de EDOs para Modelar Epidemias con Python
En la clase pasada vimos que un sistema de EDOs para describir la propagación de una epidemia está dado, en su forma más simple por:
Donde β es una constante cuyo significado revelaremos a continuación. Para poder resolver estas ecuaciones y obtener las funciones que nos permitirán predecir cómo evolucionará la epidemia usaremos métodos muy populares en matemáticas. Puedes profundizar en este curso de ecuaciones diferenciales aquí. En nuestro caso particular usaremos trucos poco comunes pero necesarios para esta situación específica.
Revisemos la solución en un paso a paso detallado:
Suponemos que la población total N es constante por lo que siempre se cumple que S(t) + I(t)=N. De esto despejamos así S(t) = N - I(t) y lo usamos para obtener una sola ecuación con una sola incógnita.
De lo cual la segunda ecuación diferencial se convierte en:
La ecuación que acabamos de obtener pertenece a una clase especial de ecuaciones denominada ecuaciones diferenciales de Bernoulli, y para estas podemos definir una nueva variable así:
Donde el cambio para la derivada se traduce en y'(t) = -I'(t)/I(t)². Ahora, reemplazando todo lo anterior en la ecuación [3] nos permite obtener la siguiente ecuación lineal:
La cual se resuelve fácilmente con el método de separación de variables para obtener una solución con la forma:
Donde, A es una constante que se determina por la condición inicial del problema. En nuestro caso la condición inicial problema es la cantidad de personas infectadas al inicio de nuestro sistema. Lo usual es definir que al principio siempre existe al menos 1 persona infectada, y esto lo podemos definir así I(t) = 1
Al reemplazar la condición inicial en la ecuación [5] obtenemos ahora:
Finalmente, de la ecuación [6] podemos obtener también la función para el número de individuos susceptibles con la relación:
Y así, hemos obtenido funciones que nos permiten saber cómo evolucionan las variables del sistema dinámico para cualquier valor del tiempo. Ahora, igual que hicimos con el ejercicio del péndulo, podemos escribir código en Python que nos permita visualizar las predicciones que se derivan de estas funciones.
Primero importamos las librerías necesarias y definimos los parámetros y funciones de nuestro modelo:
Así, hemos definido la población total con 100 personas y un solo individuo inicialmente infectado, mientras que el valor de β lo hemos definido de forma arbitraria, por ahora. Luego creamos listas de elementos que constituyen las tablas de datos sobre las cuales graficamos nuestras funciones:
En el ejemplo vemos que estamos tomando el tiempo desde 0 hasta 49 días y así mismo podemos ver la otra lista con el número de infectados total, día a día. Aquí lo extraño es ver que un día tenga 1.80 infectados y estos números no deben preocuparnos pues son simplemente aproximaciones que se derivan en que estamos usando una variable continua para aproximar un fenómeno que es discreto. Así por ejemplo, que en el día 2 haya 1.80 infectados y luego en el día 3 haya 2.42 infectados nos indica que aproximadamente le toma a la epidemia entre dos y tres días alcanzar un total de dos personas infectadas. Ahora, una visualización completa del comportamiento en el tiempo la obtenemos así:
Donde el eje horizontal es el tiempo en días y el eje vertical es el número de individuos. Varios hechos observamos de este gráfico:
Después de aproximadamente 30 días, toda la población resulta afectada.
El número de infectados siempre aumenta porque no estamos asumiendo dentro del modelo la posibilidad de recuperación de los individuos.
Como ves, las predicciones que se derivan de este modelo son algo obvias pero totalmente coherentes a la luz de las simplificaciones que incluimos en nuestro modelo.
Ahora te propongo un ejercicio adicional, calculando las soluciones de nuestro sistema para diferentes valores del parámetro β como se ve a continuación:
En este caso solo usamos las curvas de población de infectados. Lo que vemos es que a medida que el parámetro β crece , la curva alcanza el máximo en menos tiempo, esto indica que el parámetro β es una medida de la rapidez de propagación de la epidemia y ese es el modo en que debemos interpretar este número.
El hecho de que podamos predecir el valor de una variable para cualquier instante de tiempo futuro o pasado es lo que entendemos por un sistema determinístico, ya que su evolución está completamente determinada por una regla claramente definida. Es así, como el lenguaje de las ecuaciones diferenciales, aplicadas a estas situaciones, nos permite codificar el concepto de determinismo en números concretos. En este link puedes ver el notebook completo con el código.
Me llamó la antención el artificio de la ecuación de Bernoulli. Les comparto este video para ententer mejor este tipo de ecuaciones deferenciales:
https://www.youtube.com/watch?v=Z-78QrHg2fM
Gracias por el aporte!
Estoy disfrutando mucho este taller, repasando y aprendiendo cosas geniales.
¡Tengo una duda!
En este problema, ¿cómo se elige el valor final de la rapidez de propagación (beta)?
¿Simplemente lo comparamos con los datos reales y decimos que ese valor se ajusta?
Saludos Pacho :)
Tengo cuaderno y lapiz para escribir, dibujar y graficar el ejercicio como lo voy entendiendo para darle forma a lo que estoy leyendo y me está ayudando a aprender y comprender lo que estoy haciendo.
Hola, alguien me podría por fa explicar como se resolvió la ecuación diferencial de Bernoulli [3]? Es no me es muy clara la estructura.
Gracias :)
Buen día Andrés,
En el caso de la ecuación [3] lo que se busca es encontrar una sustitución que permita reducir el orden de la ecuación a una ecuación lineal de primer orden. Esto se consigue gracias a un teorema el cual puedes encontrar mas detallado aquí.
La sustitución que propone el teorema es la misma que se indica en el desarrollo del curso, es decir y(t) = 1 / I(t); una vez transformada la ecuación con la sustitución tenemos una nueva ecuación [4] de primer orden.
A continuación te dejo mi desarrollo de la ecuación [4] por medio del método de variables separables:
Reescribe las derivadas con la notación de Leibniz
(dy/dt) = β (1 - Ny) [4]
Organiza la ecuación en términos de dy y dt
(dy / β (1 - Ny)) = dt
Resuelve la integral de cada uno de los lados
(1/-Nβ) ln|1-Ny| = t + C
Aplica exponencial a cada uno de los lados
(1 - Ny) ^ (1/-Nβ) = e ^ (t + C)
Despeja y
1 - Ny = ((e ^ t) (e ^ C)) ^ (-Nβ)
1 - Ny = A (e ^ (-tNβ)) Nota: Vamos a llamar A a (e ^ C) ^ (-Nβ) ya que es una constante
Ny = (1 - A(e ^ (-Nβt)))
y = (1 - A(e ^ (-Nβt))) / N
Deshacemos la sustitución original
I(t) = N / (1 - A(e ^ (-Nβt)))
Espero haber respondido tu pregunta, hasta pronto :smile:
¡Nunca dejes de aprender!
Muchas gracias @dafmontenegro ! por la explicación y el tiempo, me aclaraste la pregúnta :)