Euler vs RK4: precision importa
Por que el metodo de Euler fracasa en sistemas caoticos y como RK4 lo resuelve
Cuando Euler no basta
En la leccion anterior vimos que Euler introduce un error de \(\mathcal{O}(h)\) por cada unidad de tiempo simulada. Para una ecuacion simple como \(dx/dt = x\), esto produce una desviacion gradual y predecible. Pero en un sistema caotico, la situacion es dramaticamente peor.
Recordemos la propiedad fundamental del caos: sensibilidad a condiciones iniciales. Dos trayectorias que empiezan infinitesimalmente cerca divergen exponencialmente. Esto significa que cualquier error numerico -- por pequeno que sea -- se amplifica como si fuera una perturbacion en las condiciones iniciales.
Con Euler aplicado a Lorenz, cada paso introduce un error que el sistema amplifica exponencialmente. Despues de unos cientos de pasos, la trayectoria numerica ya no tiene ninguna relacion con la trayectoria real. No es que se desvie un poco: se va a un ala completamente diferente del atractor.
Error acumulado en caos
En un sistema con exponente de Lyapunov \(\lambda > 0\), una perturbacion \(\delta_0\) crece como:
$$\delta(t) \approx \delta_0 \cdot e^{\lambda t}$$Para el atractor de Lorenz, \(\lambda \approx 0.9\). Si Euler introduce un error de \(\epsilon\) por paso, despues de un tiempo \(t\), el error acumulado no crece linealmente sino exponencialmente:
$$\text{Error}(t) \sim \epsilon \cdot e^{0.9t}$$Con \(dt = 0.01\) y Euler, \(\epsilon \approx 10^{-2}\). Despues de tan solo \(t = 20\) unidades de tiempo (2000 pasos), el error es del orden de \(10^{-2} \cdot e^{18} \approx 6.6 \times 10^5\). El error es un millon de veces mayor que la escala del atractor. La solucion numerica es basura.
Runge-Kutta de orden 4
El metodo de Runge-Kutta de cuarto orden (RK4) es el caballo de batalla de la integracion numerica. Fue desarrollado por los matematicos alemanes Carl Runge y Martin Kutta alrededor de 1900. Su idea central es evaluar la derivada en cuatro puntos dentro de cada paso, en lugar de uno solo.
Donde Euler pregunta "cual es la pendiente aqui?" y avanza, RK4 pregunta: "cual es la pendiente al inicio, a mitad de camino (dos veces, de formas ligeramente diferentes), y al final?" Luego combina las cuatro respuestas con una media ponderada inteligente:
$$x_{n+1} = x_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)$$donde:
$$k_1 = f(t_n,\; x_n)$$ $$k_2 = f\!\left(t_n + \tfrac{h}{2},\; x_n + \tfrac{h}{2}k_1\right)$$ $$k_3 = f\!\left(t_n + \tfrac{h}{2},\; x_n + \tfrac{h}{2}k_2\right)$$ $$k_4 = f\!\left(t_n + h,\; x_n + hk_3\right)$$El resultado es un metodo de orden 4: el error local es \(\mathcal{O}(h^5)\) y el error global es \(\mathcal{O}(h^4)\). Exploraremos la derivacion completa en la leccion 13.
Comparacion directa
La diferencia entre \(\mathcal{O}(h)\) y \(\mathcal{O}(h^4)\) parece modesta escrita asi. Pero cuando sustituimos numeros concretos, la ventaja de RK4 es aplastante:
$$\text{Para } h = 0.01: \quad \text{Euler} \sim 10^{-2}, \quad \text{RK4} \sim 10^{-8}$$RK4 es un millon de veces mas preciso que Euler con el mismo paso temporal. Volviendo al calculo del error en Lorenz: con RK4 y \(dt = 0.01\), el error por paso es \(\sim 10^{-8}\). Despues de \(t = 20\):
$$\text{Error}_{RK4}(20) \sim 10^{-8} \cdot e^{18} \approx 0.66$$Todavia hay divergencia (es caos, al fin y al cabo), pero la trayectoria de RK4 es fiel a la real durante mucho mas tiempo. Mientras Euler pierde la pista despues de unos pocos segundos de tiempo simulado, RK4 puede seguir la trayectoria correcta durante decenas de unidades de tiempo.
En sistemas caoticos, Euler diverge de la solucion real en pocas iteraciones. RK4, con su error \(\mathcal{O}(h^4)\), mantiene la fidelidad millones de veces mas tiempo para el mismo paso temporal.
Coste computacional
Hay un precio por la precision de RK4: necesita 4 evaluaciones de la funcion por paso, mientras Euler solo necesita 1. Pero esta comparacion es enganosa. Lo que realmente importa es la precision por unidad de coste computacional.
Si quieres la misma precision que RK4, Euler necesita un paso temporal mucho mas pequeno. Para igualar el error \(\mathcal{O}(h^4)\) de RK4 con \(h = 0.01\), Euler necesitaria \(h \approx 10^{-8}\). Eso es un millon de veces mas pasos, cada uno con una evaluacion. Resultado: Euler cuesta un millon de evaluaciones donde RK4 cuesta 4.
| Metodo | Orden | Eval/paso | Error (dt=0.01) | Estable en Lorenz? |
|---|---|---|---|---|
| Euler | 1 | 1 | ~10-2 | No (diverge rapido) |
| Euler mejorado | 2 | 2 | ~10-4 | Marginal |
| RK4 | 4 | 4 | ~10-8 | Si |
| RK45 (adaptativo) | 4-5 | 6 | ~10-10 | Si (optimo) |
En la practica, RK4 con \(dt = 0.005\) es la eleccion estandar para la mayoria de atractores caoticos 3D. Ofrece un equilibrio excelente entre precision, estabilidad y velocidad. ChaosLab usa exactamente esta configuracion.
Ejercicios
- Explica con tus palabras por que \(\mathcal{O}(h^4)\) es mucho mejor que \(\mathcal{O}(h)\) cuando \(h = 0.01\). Cuantos ordenes de magnitud de diferencia hay?
- Si integras Lorenz con Euler usando \(dt = 0.01\) durante 1000 pasos (\(t = 10\)), y luego con RK4 el mismo intervalo, estima el orden de magnitud del error acumulado en cada caso.
- Experimentando con el LAB, encuentra el valor aproximado de \(dt\) a partir del cual Euler produce una trayectoria claramente incorrecta para Lorenz (la trayectoria se sale del atractor o se va al infinito).
LAB: Euler vs RK4 en Lorenz
Dos integraciones del sistema de Lorenz desde la misma condicion inicial. Izquierda: Euler. Derecha: RK4. Ajusta dt y observa como Euler diverge mientras RK4 permanece estable.