Derivacion de Runge-Kutta 4
Las cuatro pendientes que dominan la simulacion
La idea: de Simpson a RK4
El metodo de Euler evalua la pendiente en un solo punto y avanza en linea recta: una evaluacion, orden 1. El metodo de Heun (y el del punto medio) mejoran evaluando la pendiente en dos puntos, alcanzando orden 2. La pregunta natural es: que pasa si evaluamos en mas puntos?
Runge-Kutta de cuarto orden (RK4) lleva esta idea a su conclusion natural: usa cuatro evaluaciones de la pendiente por paso, cuidadosamente posicionadas y ponderadas. El resultado es un error local de \(\mathcal{O}(h^5)\), lo que significa que cada paso individual es extraordinariamente preciso.
La clave esta en la ponderacion. Los pesos no son arbitrarios: provienen de la regla de Simpson para integracion numerica. Simpson aproxima la integral de una funcion usando tres puntos (inicio, medio, final) con pesos \(1:4:1\). RK4 adapta esta idea a las EDO, evaluando el punto medio dos veces (con correcciones sucesivas) para capturar la curvatura de la solucion con notable fidelidad.
Las cuatro evaluaciones
Dada una EDO \(\frac{dy}{dt} = f(t, y)\) con estado actual \((t_n, y_n)\) y paso temporal \(h\), RK4 calcula cuatro pendientes intermedias:
$$k_1 = f(t_n,\; y_n)$$\(k_1\) es la pendiente al inicio del intervalo. Identica a la evaluacion de Euler. Es nuestra primera aproximacion de la velocidad del sistema.
$$k_2 = f\!\left(t_n + \frac{h}{2},\; y_n + \frac{h}{2}\,k_1\right)$$\(k_2\) es la pendiente en el punto medio, estimada usando \(k_1\) para llegar alli. Usamos la pendiente inicial para predecir donde estaremos a mitad de camino, y luego evaluamos la derivada en ese punto predicho. Esta pendiente captura informacion sobre la curvatura que \(k_1\) por si sola no contiene.
$$k_3 = f\!\left(t_n + \frac{h}{2},\; y_n + \frac{h}{2}\,k_2\right)$$\(k_3\) es una segunda evaluacion en el punto medio, pero usando \(k_2\) para llegar en lugar de \(k_1\). Esta es la innovacion clave de RK4: la primera estimacion del punto medio no es perfecta, asi que la usamos para obtener una prediccion corregida. \(k_3\) evaluada en este punto mejorado captura efectos de orden superior, incluyendo informacion sobre la tercera derivada.
$$k_4 = f\!\left(t_n + h,\; y_n + h\,k_3\right)$$\(k_4\) es la pendiente al final del intervalo, estimada avanzando un paso completo con la pendiente corregida \(k_3\). Nos da informacion sobre el comportamiento de la derivada al otro extremo del paso.
Finalmente, combinamos las cuatro pendientes con la media ponderada:
$$y_{n+1} = y_n + \frac{h}{6}\left(k_1 + 2k_2 + 2k_3 + k_4\right)$$Los pesos \(1, 2, 2, 1\) (que suman 6, de ahi el denominador) no son arbitrarios. Provienen directamente de la regla de Simpson para integracion numerica, como veremos a continuacion.
Interpretacion geometrica
Cada pendiente \(k_i\) representa una "propuesta" de direccion para avanzar:
\(k_1\) -- pendiente al inicio: Miramos donde estamos ahora. Si avanzaramos usando solo esta pendiente, tendriamos Euler. Pero sabemos que la curva real gira, asi que necesitamos mas informacion.
\(k_2\) -- punto medio via \(k_1\): Damos un medio paso tentativo usando \(k_1\) y miramos la pendiente alli. Esta pendiente refleja mejor la curvatura en la primera mitad del intervalo.
\(k_3\) -- punto medio corregido via \(k_2\): Repetimos el medio paso, pero usando \(k_2\) como guia. Esto refina nuestra estimacion del punto medio. En cierto sentido, \(k_3\) es la pendiente que hubieramos obtenido si hubieramos usado Heun para llegar.
\(k_4\) -- pendiente al final via \(k_3\): Avanzamos un paso completo con \(k_3\) y miramos la pendiente al final. La media ponderada de las cuatro captura la curvatura de la solucion con extraordinaria fidelidad.
Los pesos 1:2:2:1 vienen de la regla de Simpson. Las dos evaluaciones en el punto medio (\(k_2 + k_3\)) suman peso 4, que es el peso central de Simpson. Los extremos (\(k_1\) y \(k_4\)) reciben peso 1 cada uno.
Conexion con la cuadratura de Simpson
Toda EDO tiene una forma integral equivalente:
$$y(t_n + h) = y(t_n) + \int_{t_n}^{t_n + h} f(t, y(t))\, dt$$El problema se reduce a aproximar esta integral. La regla de Simpson de \(1/3\) aproxima una integral con tres puntos equiespaciados:
$$\int_a^b g(t)\, dt \approx \frac{b - a}{6}\left[g(a) + 4\,g\!\left(\frac{a+b}{2}\right) + g(b)\right]$$Los pesos son \(1:4:1\), sumando 6. En RK4, evaluamos el punto medio dos veces (\(k_2\) y \(k_3\)), cada una con peso 2. La suma \(2 + 2 = 4\) coincide con el peso central de Simpson:
$$\frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4) = \frac{h}{6}\!\left(1 \cdot k_1 + 4 \cdot \frac{k_2 + k_3}{2} + 1 \cdot k_4\right)$$RK4 es, en esencia, la regla de Simpson aplicada a la forma integral de la EDO. La complicacion es que el integrando \(f(t, y(t))\) depende de la solucion desconocida. Las evaluaciones sucesivas (\(k_1 \to k_2 \to k_3 \to k_4\)) resuelven esta circularidad prediciendo los valores intermedios de \(y\) paso a paso.
Error y convergencia
Expandimos \(y(t_n + h)\) en serie de Taylor:
$$y(t_n + h) = y_n + hy' + \frac{h^2}{2}y'' + \frac{h^3}{6}y''' + \frac{h^4}{24}y^{(4)} + \frac{h^5}{120}y^{(5)} + \cdots$$Expandiendo cada \(k_i\) en Taylor y sustituyendo en la formula de RK4, se demuestra que RK4 coincide con la expansion exacta hasta el termino \(h^4\). La primera discrepancia aparece en el termino \(h^5\):
$$y_{n+1}^{\text{RK4}} = y(t_n + h) + \mathcal{O}(h^5)$$Este es el error de truncamiento local: \(\mathcal{O}(h^5)\) por paso. Para recorrer un intervalo fijo \(T\) necesitamos \(N = T/h\) pasos, asi que el error global es:
$$\text{Error global} = N \cdot \mathcal{O}(h^5) = \frac{T}{h} \cdot \mathcal{O}(h^5) = \mathcal{O}(h^4)$$La consecuencia practica es extraordinaria. Cada vez que reducimos \(h\) a la mitad:
Euler \(\mathcal{O}(h)\): error se reduce por factor 2.
Heun \(\mathcal{O}(h^2)\): error se reduce por factor 4.
RK4 \(\mathcal{O}(h^4)\): error se reduce por factor 16. Cada halving del paso temporal regala cuatro bits extra de precision.
| Metodo | Orden | Error local | Error global | Eval/paso | Factor al halvar h |
|---|---|---|---|---|---|
| Euler | 1 | \(\mathcal{O}(h^2)\) | \(\mathcal{O}(h)\) | 1 | 2x |
| Heun (RK2) | 2 | \(\mathcal{O}(h^3)\) | \(\mathcal{O}(h^2)\) | 2 | 4x |
| RK4 | 4 | \(\mathcal{O}(h^5)\) | \(\mathcal{O}(h^4)\) | 4 | 16x |
RK4 usa 4 evaluaciones por paso, pero su eficiencia por evaluacion es muy superior. Euler necesitaria 16 veces mas pasos (16 evaluaciones) para lograr la misma mejora que RK4 consigue con un solo halving. El coste adicional de 4 evaluaciones se paga con creces.
Por que no RK5, RK6...?
Si cuatro evaluaciones dan orden 4, por que no usar cinco para orden 5? A partir de orden 5, el numero de evaluaciones crece mas rapido que el orden. RK4 es el "punto dulce": 4 evaluaciones para orden 4 (una evaluacion por orden). Un metodo RK de orden 5 necesita al menos 6 evaluaciones, y uno de orden 6 necesita al menos 7.
Por esto RK4 ha dominado la integracion numerica durante mas de un siglo. Cuando se necesita mas precision, en lugar de subir el orden se usa un metodo adaptativo como RK45 (Dormand-Prince), que combina un metodo de orden 4 y uno de orden 5 para estimar y controlar el error automaticamente. Lo veremos en la leccion 14.
Ejercicios
- Aplica un paso de RK4 a mano sobre \(\frac{dy}{dx} = -y\) con \(y_0 = 1\), \(x_0 = 0\) y \(h = 0.2\). Calcula \(k_1, k_2, k_3, k_4\) y \(y_1\). Compara con la solucion exacta \(e^{-0.2} = 0.81873\).
- Repite el ejercicio con Euler usando el mismo \(h = 0.2\). Cual es el error relativo de cada metodo? Cuantas veces mas preciso es RK4?
- Explica por que los pesos de RK4 son \(1:2:2:1\) y no \(1:1:1:1\) (promedio simple). Que regla de cuadratura justifica esta eleccion?
LAB: Visualizacion de las 4 pendientes
Un paso de RK4 sobre \(dy/dx = -y\) (solucion exacta: \(e^{-x}\)). Las cuatro pendientes aparecen en secuencia: \(k_1\) (cyan), \(k_2\) (violeta), \(k_3\) (rosa), \(k_4\) (ambar). La media ponderada se muestra en verde. Ajusta \(h\) para ver como cambia la geometria del paso.