Fórmula RK4 vs. fórmula de Euler

8/12

Lectura

El ejemplo de ecuación diferencial que usamos en la clase pasada es muy sencillo, pero lo empleamos porque necesitábamos un ejemplo donde conociéramos la solución exacta para calcular el error real del cálculo numérico. Esta clase está dividida en dos partes:

  1. En la primera parte, veremos un método que es mucho más eficiente que Euler y que además es uno de los más usados en la vida real.
  2. Luego en la segunda parte, resolveremos el sistema de Lorenz que vimos en clases pasadas el cual representa un modelo de la convección. También discutiremos porque la convección es un claro ejemplo de caos a la luz de una definición apropiada que daremos en esta sección.

Métodos predictor-corrector

El método de Euler es la forma más sencilla de resolver numéricamente una ecuación diferencial de la forma:

CleanShot 2020-07-23 at 02.02.25.png

Donde recordemos que se ejecuta con la fórmula iterativa:

CleanShot 2020-07-23 at 03.29.43.png

A partir de este método surgen varias mejoras, incluyendo los métodos predictor-corrector que ejecutan un paso inicial “predictor” y posteriormente ejecutan pasos adicionales “corrector” o correctivos. En este curso no profundizaremos en la justificación detallada detrás de estos métodos, solamente haremos una comparación con uno de los métodos más eficientes que es el Runge-Kutta de orden 4, cuya fórmula iterativa está dada en cuatro pasos así:

Screenshot_51.png

Haremos una comparación entre los métodos Runge-Kutta de orden 4 o RK4 y el método de Euler para ver cómo mejora la precisión numérica de un método a otro. Para ello implementaremos funciones que ejecuten paso a paso ambos algoritmos de una forma diferente a lo que hemos hecho antes (esto con el fin de ser más transparentes a la hora de comparar con las fórmulas matemáticas):

pasted image 0 (7).png

Observa que en el código:

  1. Dejé por fuera de todos los pasos “predictor-corrector” el 𝚫t, y lo puse al final en la fórmula iterativa que calcula yₙ₊₁. Matemáticamente son cosas equivalentes en virtud de la propiedad de factorización algebraica, así que el resultado es el mismo.
  2. Lo que en las fórmulas se llama yₙ, en el código es y₀. Así mismo yₙ se corresponde con y.
  3. Cada función ejecuta un solo paso de cada algoritmo, lo que quiere decir que debemos usarlas dentro de una estructura tipo for() o while() para ejecutar varios pasos sucesivamente.

Ahora, para comparar ambos métodos tomaremos como ejemplo el caso de la clase pasada donde f (t, y) = y cuya solución con el valor inicial y(0) = 1 es:

CleanShot 2020-07-23 at 03.37.32.png

Vamos a calcular la solución numérica por ambos métodos, para varios intervalos de tiempo y compararemos con la solución exacta, así:

pasted image 0 (8).png

El final del código son líneas que tienen que ver solamente con la visualización de las soluciones numéricas y la exacta, cuyo resultado se ve así:

Unknown-2.png

Y de esta visualización notamos inmediatamente que para un mismo 𝚫t el rendimiento de RK4 es superior al método de Euler (compara las curvas en color rojo que representan ambas soluciones para 𝚫t = 0.556. Así como hicimos en la clase pasada aquí también podemos calcular el error acumulado, ahora para ambos métodos dándonos como resultado lo siguiente:

image20.png
image22.png

El gráfico resultante nos muestra claramente que el error para RK4, es al menos dos órdenes de magnitud más pequeño que el de Euler, y esto refleja la evidente superioridad de RK4 incluso para intervalos 𝚫t relativamente grandes.

En nuestra próxima clase, veremos cómo usar RK4 para resolver sistemas complejos que presentan comportamiento caótico y el notebook de esta clases lo encuentras en este link.

Aportes 4

Preguntas 0

Ordenar por:

Los aportes, preguntas y respuestas son vitales para aprender en comunidad. Regístrate o inicia sesión para participar.

Admito que estoy confundido y tengo que trabajar en esta parte del contenido.
¿Hay algún libro o taller, video o alternativa para profundizar?
Estaría muy agradecido, me parece muy interesante

¿qué significa o qué implica los valores ki en el ajuste Runge-Kutta de orden 4?

Codigo.

import numpy as np 
import matplotlib.pyplot as plt 

def f(t,y): return y

def rk4(t0,y0,dt,f):
    k1 = f(t0, y0)
    k2 = f(t0 + dt / 2.0, y0 + dt * k1 / 2.0)
    k3 = f(t0 + dt / 2.0, y0 + dt * k2 / 2.0)
    k4 = f(t0 + dt, y0 + dt * k2)
    y = y0 + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
    return y

def euler(t0, y0, dt, f):
    f0 = f(t0, y0)
    y = y0 + dt * f0
    return y

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

if __name__=='__main__':
    tmax = 5
    plt.figure(figsize=(10,10))
    for n, c in zip([10,100], ['red','blue']):
        ys_rk4 = [1]
        ys_euler = [1]
        ts = np.linspace(0,tmax,n)
        dt = ts[1] - ts [0]
        for i in range(n-1):
            u_rk4 = rk4(ts[i], ys_rk4[i], dt,f)
            u_euler = euler(ts[i], ys_euler[i], dt, f)
            ys_rk4.append(u_rk4)
            ys_euler.append(u_euler)
        plt.plot(ts, ys_euler, '--', color=c, label='Euler dt = {}'.format(round(dt,3)))
        plt.plot(ts, ys_rk4, color=c, label='RK4 dt = {}'.format(round(dt,3)))
    plot_text_size = 20
    plt.plot(ts, exact_sol(ts), '-k', label='exact')
    plt.xticks(fontsize=plot_text_size)
    plt.yticks(fontsize=plot_text_size)
    plt.legend(fontsize = plot_text_size)
    plt.show()

    tmax = 100
    dt_arr = np.arange(0.0001, 0.5, 0.0001)
    error_rk4 = []
    error_euler = []
    for dt in dt_arr:
        ys_rk4 = [1]
        ys_euler = [1]
        ts = np.arange(0, tmax, dt)
        for  i in range(len(ts)):
            u_rk4 = rk4(ts[i], ys_rk4[i], dt, f)
            u_euler = euler(ts[i], ys_euler[i], dt, f)
            ys_rk4.append(u_rk4)
            ys_euler.append(u_euler)
        ys_exact = exact_sol(ts)
        error_rk4.append(np.abs(ys_rk4[-1]-ys_exact[-1]))
        error_euler.append(np.abs(ys_euler[-1]-ys_exact[-1]))

    plt. figure(figsize=(10,8))
    plot_text_size = 20
    plt.plot(dt_arr, error_rk4, label='RK4')
    plt.plot(dt_arr, error_euler, label='Euler')
    plt.plot(dt_arr[:-4900],error_rk4[:-4900], '-*r')
    plt.xticks(fontsize=plot_text_size)
    plt.yticks(fontsize=plot_text_size)
    plt.legend(fontsize=plot_text_size)
    plt.yscale('log')
    plt.xscale('log')
    plt.show()