Modelos numéricos al rescate

6/12

Lectura

Los modelos cuya solución está basada en un algoritmo numérico los denominamos modelos numéricos. El algoritmo numérico más sencillo que existe es el método de Euler el cual discutiremos a continuación.

Método de Euler

Cuando tenemos una ecuación diferencial de primer orden para una función desconocida y(t) lo más general posible, escribimos:

CleanShot 2020-07-23 at 01.05.35.png

Donde f(y, t) representa una expresión que contiene todos los demás términos de la ecuación diferencial que no acompañen a la derivada. Ahora, como vimos con el problema de la propagación de una epidemia, un sistema dinámico debe tener un punto de inicio o condición inicial. En general siempre debemos definir cual es esta condición inicial para nuestro sistema así:

CleanShot 2020-07-23 at 01.09.47.png

Con y₀ siendo el valor inicial de la variable que estamos buscando predecir. En este sentido, el método de Euler propone construir la solución de la ecuación diferencial paso a paso suponiendo que la derivada se aproxima por medio de diferencias finitas pero muy pequeñas:

modelos-numericos-6.png

Estas diferencias finitas representan intervalos de tiempo tan pequeños como uno los considere y paso a paso se puede ir calculando el siguiente valor de la variable y con base en el valor del tiempo anterior:

CleanShot 2020-07-23 at 01.10.27.png

Esto representa un proceso iterativo donde vamos construyendo la solución paso a paso. El precio de esto es que el aproximar una derivada por intervalos finitos implica que en cada paso se genera un error que se acumula en los siguientes pasos. Matemáticamente hablando, la única manera de reducir este error es considerando un intervalo de tiempo 𝚫t ➝ 0, pero computacionalmente no podemos considerar números infinitamente pequeños, entonces procuramos tomarlo tan pequeño como nuestros recursos de cómputo nos lo permitan. Esta es la gran diferencia entre el cálculo de infinitesimales y los métodos numéricos.

Dado que los métodos numéricos son muy sencillos conceptualmente hablando, pero pueden representar procesos que deban repetirse muchas veces, el uso de un computador para su ejecución es lo ideal y por lo tanto hoy en día todos los métodos numéricos se implementan en un lenguaje de programación para poder realizar múltiples pasos en tiempos razonablemente cortos y dejándole todo el trabajo a una máquina. Como dato curioso estos algoritmos eran calculados a mano antes de que aparecieran las grandes computadoras, pero ya te imaginarás lo engorroso que eso era en aquella época.

Solución de Euler en Python

Tomaremos como ejemplo una ecuación diferencial sencilla junto con una condición inicial que usaremos como base para aplicar el método de Euler:

CleanShot 2020-07-23 at 01.12.15.png

Es decir esto representa un problema donde la solución será una función que es igual a su derivada y que su valor para t = 0 sea 1. El algoritmo de Euler en este caso quedaría de la siguiente manera:

CleanShot 2020-07-23 at 01.13.08.png

Para inicial el algoritmo debemos definir el intervalo de tiempo, digamos que tomamos 𝚫t = 0.01 y calculamos los pasos así:

CleanShot 2020-07-23 at 01.14.16.png

Y así sucesivamente hasta el número de pasos que sea necesario. Por ejemplo si necesitamos construir la solución hasta un tiempo t = 10 , requerimos calcular al menos 1000 pasos con un 𝚫t = 0.01. Esto es mejor hacerlo en un computador, que es justamente lo que haremos a continuación en las siguientes líneas de código:

pasted image 0 (4).png

De aquí, la solución numérica correspondiente se vé así:

image24.png

Y para los que ya saben que la solución exacta (la que obtenemos de los métodos del cálculo) de esta ecuación ecuación diferencial es la función exponencial:

CleanShot 2020-07-23 at 01.17.40.png

Se darán cuenta que al comparar ambos tipos de soluciones, se observa una pequeña diferencia que va creciendo conforme avanzamos en el tiempo y que es mayor si escogemos intervalos de tiempo más largos:

image26.png

Y si quieres seguir jugando con este código te darás cuenta de que esta diferencia desaparece a medida que el intervalo de tiempo es cada vez más pequeño, solo que cada vez tendrás que realizar más pasos.

Y es así como vemos una alternativa para obtener soluciones a modelos de sistemas dinámicos a partir del paradigma numérico, con la advertencia que los métodos numéricos tienen el problema que el error es directamente proporcional al intervalo de tiempo que se usa y esto, a su vez, es inversamente proporcional al número de operaciones que deben realizarse. Es decir, que para obtener mejores soluciones es preciso realizar más pasos y por lo tanto consumir más recursos computacionales.

El notebook de esta clase lo encuentras en este link de Google Colab. En nuestra próxima clase discutiremos en detalle las diferencias entre las soluciones exactas que obtenemos con álgebra y cálculo contra las soluciones numéricas.

Aportes 5

Preguntas 0

Ordenar por:

¿Quieres ver más aportes, preguntas y respuestas de la comunidad? Crea una cuenta o inicia sesión.

Muy buen aporte respecto a los métodos numéricos, sería excelente ampliar o crear un curso de análisis numérico y optimización, para la ruta de data science. Pues tienen grandes aplicaciones y al parecer solo los que somos matemáticos estamos al tanto del tema

Por si no queda claro como llegamos a la formula.

import numpy as np 
import matplotlib.pyplot as plt 

def exact_sol(ts): return np.exp(ts)

def num_sol(ts, dt, tf=10, y0=1):
    ys = [y0]
    ts = [0]
    num_steps = int(tf/dt)
    for _ in range(num_steps):
        ts.append(ts[-1]+dt)
        ys.append((1+dt)*ys[-1])
    return ts, ys


 
if __name__=='__main__':

    ys = [1] # creamos un areglo inicial de valores en y
    ts = [0] # creamos un areglo inicial de valores en t
    dt = 0.01

    num_steps = 1000
    for i in range(num_steps):
        ts.append(ts[-1]+dt) # calculamos el proximo t y lo agregamos al arreglo
        ys.append((1+dt)*ys[-1]) # calculamos el proximo y y lo agregamos al arreglo

    #print (ts)
    #print(ys)

   
    plt.figure(figsize=(15,8))
    plt.plot(ts,ys)
    plt.xticks(fontsize=20)
    plt.yticks(fontsize=20)
    plt.show()

    plt.figure(figsize=(12,6))
    for dt in [0.1, 0.05, 0.01]:
        ts, ys = num_sol(ts, dt)
        plt.plot(ts, ys, '--', label='dt = {}'.format(dt))
    
    plt.plot(ts,exact_sol(ts), label='exact')
    plt.xticks(fontsize=20)
    plt.yticks(fontsize=20)
    plt.legend(fontsize=20)
    plt.show()

Genial!