10.5: Métodos Runge-Kutta

Los tres métodos que hemos encuestado hasta ahora (Forward Euler, Atrás Euler y Adams-Moulton) han consistido en muestrear la función derivada\(F(y,t)\) en uno de los pasos\(\{t_n\}\) de tiempo discretos y las soluciones en esos pasos de tiempo\(\{\vec{y}_n\}\). Es natural preguntarse si podemos mejorar la precisión muestreando la función derivada a valores “intermedios” de\(t\) y\(\vec{y}\). Esta es la idea básica detrás de una familia de métodos numéricos conocidos como métodos Runge-Kutta.

Aquí hay una versión simple del método, conocido como Runge-Kutta de segundo orden (RK2). Nuestro objetivo es sustituir el término derivado por un par de términos, de la forma

\[\vec{y}_{n+1} = \vec{y}_n + A h \vec{F}_A + B h \vec{F}_B,\]

donde

\[\begin{align}\vec{F}_A &= \vec{F}(\vec{y}_n, t_n)\\\vec{F}_B &= \vec{F}(\vec{y}_n + \beta \vec{F}_A, t_n + \alpha).\end{align}\]

Los coeficientes\(\{A, B, \alpha, \beta\}\) son parámetros ajustables cuyos valores elegiremos en breve, para minimizar el error de truncamiento local.

Durante cada paso de tiempo, comenzamos conociendo la solución\(\vec{y}_n\) a la vez\(t_{n}\), primero calculamos\(\vec{F}_{A}\) (que es el término derivado que entra en el método Forward Euler); luego lo usamos para calcular un término derivado “intermedio”\(\vec{F}_{B}\). Finalmente, utilizamos un promedio ponderado de\(\vec{F}_{A}\) y\(\vec{F}_{B}\) como término derivado para el cálculo\(\vec{y}_{n+1}\). A partir de esto, es evidente que se trata de un método explícito: para cada una de las subecuaciones, los “lados de la derecha” contienen cantidades conocidas.

Ahora tenemos que determinar los valores apropiados de los parámetros\(\{A, B, \alpha, \beta\}\). Primero, Taylor nos expandimos\(\vec{y}_{n+1}\)\(t_{n}\), usando la regla de la cadena:

\[\begin{align} \vec{y}_{n+1} &= \vec{y}_n + h \left.\frac{d\vec{y}}{dt}\right|_{t_n} + \frac{h^2}{2} \left.\frac{d^2\vec{y}}{dt^2}\right|_{t_n} + O(h^3)\\&= \vec{y}_n + h \vec{F}(\vec{y}_n, t_n) + \frac{h^2}{2} \left[ \frac{d}{dt}\vec{F}(\vec{y}(t), t)\right]_{t_n} + O(h^3) \\&= \vec{y}_n + h \vec{F}(\vec{y}_n, t_n) + \frac{h^2}{2} \left[ \sum_j \frac{\partial \vec{F}}{\partial y_j}\, \frac{dy_j }{dt} + \frac{\partial \vec{F}}{\partial t}\right]_{t_n} + O(h^3) \\&= \vec{y}_n + h \vec{F}_A + \frac{h^2}{2} \left\{ \sum_j \left[\frac{\partial \vec{F}}{\partial y_j} \right]_{t_n} \!\! F_{Aj} \; +\; \left[\frac{\partial \vec{F}}{\partial t}\right]_{t_n}\right\} + O(h^3) \end{align}\]

De la misma manera, nosotros Taylor expandimos el término derivado intermedio\(F_{B}\), cuya fórmula se dio anteriormente:

\[F_B = F_A + \beta \sum_j F_{Aj} \left[\frac{\partial F}{\partial y_j}\right]_{t_n} + \alpha \left[\frac{\partial F}{\partial t}\right]_{t_n}.\]

Si comparamos estas expansiones de Taylor con la fórmula RK2, entonces se puede ver que los términos se pueden hacer para que coincidan hasta (e incluyendo)\(O(h^{2})\), si se eligen los parámetros para obedecer las ecuaciones

\[A + B = 1, \quad \alpha = \beta = \frac{h}{2B}.\]

Un posible conjunto de soluciones es\(A = B = 1/2\) y\(\alpha = \beta = h\). Con estas condiciones cumplidas, el método RK2 tiene un error de truncamiento local de\(O(h^{3})\), un orden mejor que el Método Forward Euler (que también es un método explícito), y comparable al Método Adams-Moulton (que es un método implícito).

El error de truncamiento local se puede reducir aún más tomando muestras más intermedias de la función derivada. El método Runge-Kutta más utilizado es el método Runge Kutta de cuarto orden (RK4), que viene dado por

\[\begin{align} \vec{y}_{n+1} &= \vec{y}_n + \frac{h}{6}\left(\vec{F}_A + 2\vec{F}_B + 2\vec{F}_C + \vec{F}_D \right)\\ \vec{F}_A &= \vec{F}(\vec{y}_n,\, t_n), \\ \vec{F}_B &= \vec{F}(\vec{y}_n + \tfrac{h}{2} \vec{F}_A,\, t_n + \tfrac{h}{2}), \\ \vec{F}_C &= \vec{F}(\vec{y}_n + \tfrac{h}{2} \vec{F}_B,\, t_n + \tfrac{h}{2}), \\ \vec{F}_D &= \vec{F}(y_n + h\vec{F}_C,\, t_n + h). \end{align}\]

Esto tiene error de truncamiento local de\(O(h^{5})\). Es un método explícito, y por lo tanto tiene la desventaja de ser inestable si el problema es rígido y\(h\) es suficientemente grande.