Comparación de Métodos: Runge-Kutta 4 vs Euler en Ecuaciones Diferenciales
Contenido del curso
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:
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:
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.
Lorena Mondragón
Eduardo Monzón
Gracias por el resumen.
Roger Christian Cansaya Olazabal
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()
Arturo Barrios
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
Carlos Perilla
Para extender la parte teoríca: https://web.archive.org/web/20151123033428/http://ocw.um.es/ingenierias/ampliacion-de-calculo-numerico/material-de-clase-1/ivp-ode.pdf https://sites.google.com/site/metodosnumericosmecanica/home/unidad-i
Un poco más practico: https://www.myphysicslab.com/explain/runge-kutta-en.html
Juan Nuñez
El código se vuelve extenso porque el profesor hace una comparativa entre 2 métodos y la solución exacta y cada método lo corre con 2 diferentes dt con el fin de hacer el versus entre Euler y RK. Dejo el código solo con la parte de RK4
import numpy as np import matplotlib.pyplot as plt def f(t, y): return y #la funcion de retorno cambia dependiendo de la ecuacion diferencial #esta ecuación es dy/dt = y por lo que f(dt, y) es la ecuacion diferencial a resolver def rk4(t0, y0, dt): #la funcion f no es necesario pasarsela como parámetro ya que es función global kuta1 = f( t0, y0 ) kuta2 = f( t0 + dt/2 , y0 + dt*kuta1/2 ) kuta3 = f( t0 + dt/2 , y0 + dt*kuta2/2 ) kuta4 = f( t0 + dt , y0 + dt*kuta3 ) y = y0 + (dt/6)*(kuta1 + 2*kuta2 + 2*kuta3 + kuta4) return y if __name__=="__main__": tmax=50 #es el intervalo en eje X máximo al que queremos llegar ya que #al ser método numérico no podemos obtener valores para todo t n=100 #este valor se puede cambiar dependiendo del numero de pasos que querramos dar #de este valor depende el intervalo discreto dt ys=[1] # aqui se pone las condiciónes de frontera o iniciales de la ED ts=np.linspace(0, tmax, n) #esta función regresa un arreglo de n valores entre 0 y tmax espaciados equitativamente dt=ts[1]-ts[0] #dado que ts esta equitativamente espaciado se puede hacer la resta entre cualesquiera valores contiguos #cada valor generado por rk4 será almacenado en el arreglo de resultados for i in range(n-1): ys.append(rk4(ts[i], ys[i], dt)) plt.figure(figsize=(10,8)) plt.plot(ts, ys, 'o') # plot del arreglo resultado rk4 plt.plot(ts, np.exp(ts), '--') #plot del arreglo de la solución exacta plt.show()
Arturo Barrios
¿qué significa o qué implica los valores ki en el ajuste Runge-Kutta de orden 4?
Carlos Perilla
Son pendientes, y al final se realiza un promedio ponderado de las pendientes. Aquí hay más info: https://es.wikipedia.org/wiki/M%C3%A9todo_de_Runge-Kutta#M%C3%A9todos_de_Runge-Kutta_de_cuarto_orden