Comparación de Métodos: Runge-Kutta 4 vs Euler en Ecuaciones Diferenciales

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:

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

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í:

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):

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:

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í:

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í:

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:

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.