Runge–Kutta y paso adaptativo

1. Motivación: una sonda en órbita muy excéntrica

Una sonda en una órbita de excentricidad $e = 0.9$ pasa el 90 % del tiempo lejos del planeta moviéndose despacio y, durante una fracción breve, cruza el pericentro a velocidad $\sqrt{(1+e)/(1-e)} = \sqrt{19} \approx 4.4$ veces mayor que en el apocentro. Si integramos con un paso fijo pequeño para resolver el pericentro, desperdiciamos casi todo el trabajo en el apocentro; si lo elegimos grande, el pericentro se resuelve mal y la órbita se abre o se cierra. Este es el problema que resuelven los métodos de paso adaptativo: estimar el error de cada paso y dejar que el propio integrador decida cuán largo puede ser el siguiente.

En esta lección usamos el problema de Kepler plano con parámetro gravitatorio $\mu = G(m_1+m_2)$ como banco de pruebas, porque tiene solución exacta (lección sobre la ecuación de Kepler) con la que medir errores. Escribimos $y = (x, y, v_x, v_y)$ y

Unidades: en el código ($\mathtt{src/lib/math/rungekutta.ts}$) todo es adimensional con $\mu = 1$, $a = 1$; el periodo es $T = 2\pi$. Para pasar a SI basta fijar $\mu$ y $a$ y usar $[t] = \sqrt{a^3/\mu}$.

2. Métodos de un paso: error local y global

Sea $\dot y = f(t,y)$ con $f$ lipschitziana en $y$ con constante $L$, es decir $\lVert f(t,u) - f(t,v)\rVert \le L\lVert u - v\rVert$. Un método de un paso avanza con una función de incremento $\Phi$:

Definición. El error local en el paso $n$ es la diferencia entre lo que da el método partiendo del valor exacto y el valor exacto un paso después: $\ell_{n+1} = y(t_n) + h\Phi(t_n, y(t_n), h) - y(t_{n+1})$. El método es de orden $p$ si $\lVert\ell_{n+1}\rVert \le C h^{p+1}$. El error global es $e_n = y_n - y(t_n)$.

Teorema (convergencia). Si el método tiene orden $p$ y $\Phi$ es lipschitziana en $y$ con constante $\Lambda$, entonces en un intervalo fijo $[t_0, t_0 + T_f]$

Demostración. Restando la definición de error local de la del método, $e_{n+1} = e_n + h\left[\Phi(t_n, y_n, h) - \Phi(t_n, y(t_n), h)\right] + \ell_{n+1}$, luego $\lVert e_{n+1}\rVert \le (1 + h\Lambda)\lVert e_n\rVert + C h^{p+1}$. Iterando esta recurrencia lineal con $e_0 = 0$,

usando $(1 + h\Lambda)^n \le e^{nh\Lambda}$. Con $nh \le T_f$ se obtiene la cota. $\blacksquare$

Un error local $O(h^{p+1})$ se convierte en global $O(h^p)$ porque se cometen $n = T_f/h$ errores. La constante $e^{\Lambda T_f}$ es pesimista pero real: en órbitas, los errores de fase crecen linealmente con el tiempo, y en sistemas caóticos exponencialmente (lección anterior).

2.1 Euler

$\Phi = f(t_n, y_n)$. Comparando con la serie de Taylor $y(t_n + h) = y + h\dot y + \tfrac{h^2}{2}\ddot y + \dots$, el error local es $\tfrac{h^2}{2}\ddot y + O(h^3)$: orden 1. Para una órbita circular, $\ddot y$ apunta hacia el centro y Euler sale siempre «por la tangente»: la órbita crece en espiral (ejercicio 1). En nuestra tabla (§5) hace falta $h \approx 10^{-3}$ para un error del 20 % tras una sola vuelta.

2.2 Condiciones de orden: el caso de dos etapas

Un método de Runge–Kutta explícito de $s$ etapas se define por su tabla de Butcher $(c, A, b)$:

Para dos etapas con $c_1 = 0$, $c_2 = a_{21} = c$ (condición $\sum_j a_{ij} = c_i$ que asumimos siempre), y escribiendo $f$, $f_t$, $f_y$ evaluadas en $(t_n, y_n)$,

La solución exacta es $y + h f + \tfrac{h^2}{2}(f_t + f_y f) + O(h^3)$, de modo que el método es de orden 2 si y solo si

Es una familia uniparamétrica: $c = \tfrac12$ da el punto medio ($b = (0, 1)$) y $c = 1$ el método de Heun ($b = (\tfrac12, \tfrac12)$). Para orden 3 aparecen 4 condiciones, para orden 4 son 8 y para orden 5 son 17; el número de condiciones crece como los árboles con raíz de $p$ vértices (teoría de Butcher; Hairer–Nørsett–Wanner §II.2). Por eso el RK4 clásico con 4 etapas y orden 4 es un punto dulce: a partir de orden 5 hacen falta más etapas que el orden ($6$ etapas para orden 5, $7$ para orden 6).

2.3 RK4

Con $c = (0, \tfrac12, \tfrac12, 1)$, $a_{21} = a_{32} = \tfrac12$, $a_{43} = 1$ y $b = \tfrac16(1, 2, 2, 1)$ se cumplen las 8 condiciones de orden 4: error local $O(h^5)$, global $O(h^4)$. Aplicado a $\dot y = \lambda y$ da $y_{n+1} = R(h\lambda) y_n$ con $R(z) = 1 + z + z^2/2 + z^3/6 + z^4/24$, el polinomio de Taylor de $e^z$: la región de estabilidad absoluta $\lvert R(z)\rvert \le 1$ corta el eje real negativo en $z = -2.785$ y el imaginario en $\pm 2.83\,i$. Para un oscilador de frecuencia $\omega$ eso exige $h\omega \lt 2.83$; en el pericentro de una órbita la «frecuencia local» es $\sqrt{\mu/r_p^3}$, y es ella la que fija el paso.

3. Pares encajados y control del paso

3.1 La idea

Si disponemos de dos aproximaciones $y_{n+1}$ (orden $p$) e $\hat y_{n+1}$ (orden $\hat p = p - 1$) que comparten las mismas etapas $k_i$ (solo cambian los pesos $b_i$ frente a $\hat b_i$), la diferencia

estima el error local del método de orden inferior sin ninguna evaluación adicional. Los dos pares más usados:

Es el método ode45 de MATLAB y RK45 de SciPy. Las dos tablas están en rungekutta.ts; el test comprueba las identidades $\sum b_i = \sum\hat b_i = 1$, $\sum_j a_{ij} = c_i$, y mide el orden 5 de la solución propagada y el orden 5 del estimador de error (§5). Nuestra implementación no aprovecha FSAL (7 evaluaciones por paso), lo que infla un 17 % las cuentas de evaluaciones de Dormand–Prince en las tablas.

3.2 Ley de control del paso

Queremos que el error local por componente no supere una tolerancia mixta $sc_i = \mathrm{atol} + \mathrm{rtol}\cdot\max(\lvert y_{n,i}\rvert, \lvert y_{n+1,i}\rvert)$. Definimos la norma escalada

Como $\mathrm{err} \approx K h^{q}$ con $q = \hat p + 1$ (el orden del error local del método encajado de orden $\hat p$), el paso «óptimo» que daría $\mathrm{err} = 1$ es $h_{\text{opt}} = h\,\mathrm{err}^{-1/q}$. Con un factor de seguridad y límites para evitar oscilaciones:

Si $\mathrm{err} \le 1$ el paso se acepta y se avanza; si no, se rechaza y se repite con $h_{\text{nuevo}}$. Es un controlador proporcional (Hairer–Nørsett–Wanner §II.4, ec. 4.13). En integrarAdaptativo usamos $q = 5$ para ambos pares.

Conviene ser honesto sobre lo que se controla: el error local por paso, no el global. El error global acumulado puede ser órdenes de magnitud mayor que la tolerancia (tabla de §5: con $\mathrm{tol} = 10^{-6}$ el error de posición tras una vuelta es $\sim 10^{-3}$) porque los errores de fase se suman y porque la órbita excéntrica amplifica perturbaciones cerca del pericentro. La regla práctica es que el error global escala linealmente con la tolerancia para un método de orden $p$ propagado con extrapolación local ($\mathrm{tol}^{p/(p+1)}$ en teoría, pero en la práctica el exponente suele estar entre 0.8 y 1): la tabla da $\mathrm{tol}\, 10^{-6}\to10^{-12}$ y error $1.6\times10^{-3}\to 4\times10^{-10}$, es decir un factor $4\times 10^{6}$ en error por $10^{6}$ en tolerancia.

4. Simulador

Qué observar. El simulador integra con paso fijo (RK4 o leapfrog). Con «2 cuerpos, e = 0.8», mide cuántos pasos por periodo necesita RK4 para que la órbita se cierre visiblemente: ese número es el que un integrador adaptativo reparte de forma desigual, concentrando los pasos cortos en el pericentro (en nuestra corrida con $e = 0.9$ y $\mathrm{tol} = 10^{-10}$, el paso varía entre $4.3\times10^{-5}$ y $9.5\times10^{-2}$: un factor 2200).

5. Resultados numéricos

Todos los números de esta sección provienen de rungekutta.ts (el test rungekutta.test.ts verifica los órdenes; las tablas se generaron con el mismo código). Órbita $a = 1$, $\mu = 1$, $T = 2\pi$; error de posición $\lVert \mathbf r_N - \mathbf r(T)\rVert$ tras una vuelta frente a la solución de Kepler, y error relativo de energía $\lvert\Delta E/E\rvert$.

5.1 Orden medido, $e = 0.5$, paso fijo

Al dividir $h$ por 2, Euler divide el error por $\approx 1.5$–$2$ (orden 1, todavía fuera del régimen asintótico mientras el error es $O(1)$), RK4 por $\approx 17$–$21$ (orden 4, $2^4 = 16$) y Dormand–Prince por $\approx 29$–$37$ (orden 5, $2^5 = 32$). En la última fila Dormand–Prince ya ha tocado el suelo del redondeo en doble precisión ($\sim 10^{-13}$ tras $6400\times 7$ evaluaciones): reducir más el paso no mejora nada, porque el error de redondeo crece como $\epsilon_{\text{máq}}/h$.

La energía de RK4 deriva: con $h = T/200$ el error $\lvert\Delta E/E\rvert$ es $9.3\times10^{-7}$ tras 1 vuelta, $9.3\times10^{-6}$ tras 10 y $9.3\times10^{-5}$ tras 100: exactamente lineal en $t$, como se demostró para el oscilador en la lección anterior.

5.2 Paso adaptativo, $e = 0.9$, una vuelta

Comparación con RK4 de paso fijo en la misma órbita: para bajar el error de posición a $1.6\times10^{-7}$ hacen falta $16\,000$ pasos ($64\,000$ evaluaciones); Dormand–Prince lo consigue con tolerancia $10^{-10}$ y $1953$ evaluaciones, 33 veces menos. La ganancia es mayor cuanto más excéntrica es la órbita: en una circular el paso adaptativo apenas aporta.

Dos observaciones sobre la tabla. (1) A tolerancia laxa hay muchos rechazos (≈ 30 %): el controlador proporcional sobreestima el paso tras el pericentro; los controladores PI (Gustafsson 1991) lo reducen. (2) Dormand–Prince conserva mejor la energía que RKF45 al mismo coste porque sus coeficientes minimizan el error de la fórmula de orden 5, que es la que se propaga; RKF45 fue optimizado para la de orden 4.

6. Análisis de error: lo que el paso adaptativo no da

1. No es simpléctico. Cambiar $h$ según el estado convierte la aplicación en $y \mapsto \Phi_{h(y)}(y)$, que en general no conserva la forma simpléctica aunque $\Phi_h$ lo sea (Hairer–Lubich–Wanner §VIII.2). La energía vuelve a derivar: en la tabla, $\lvert\Delta E/E\rvert$ crece con el tiempo como en RK4. Para integraciones de miles de periodos se prefiere un simpléctico de paso fijo (o con paso variable reversible, que usa una transformación de tiempo $dt = g(y)\,d\tau$ y recupera la conservación).
2. Controla el error local, no el global. Vimos que el global puede ser $10^{3}$ veces la tolerancia. Para estimar el global hay que integrar dos veces (con $\mathrm{tol}$ y $\mathrm{tol}/10$) y comparar, o integrar hacia delante y hacia atrás.
3. Redondeo. El error total es $\approx C_1 h^{p} + C_2\,\epsilon_{\text{máq}}/h$; existe un $h$ óptimo y pedir $\mathrm{tol} \lt 10^{-13}$ en doble precisión es inútil (fila 6400 de la tabla 5.1).
4. Encuentros cercanos. El paso adaptativo los resuelve (es su gran ventaja frente al paso fijo), pero el coste se dispara: $h \sim \sqrt{r_{\min}^3/\mu}$. Para colisiones casi exactas hace falta regularizar (Kustaanheimo–Stiefel) o usar coordenadas que eliminen la singularidad.

Problemas

1. (Euler en órbita circular). Para $\dot y = f(y)$ con $\mu = 1$, $r = 1$ y velocidad circular $v = 1$, demuestra que un paso de Euler aumenta $r$ en $\tfrac{h^2}{2} + O(h^4)$ y la energía específica en $\tfrac{h^2}{2} + O(h^4)$, de modo que la órbita se abre en espiral.

2. (Condiciones de orden 2). Demuestra que no existe ningún método de Runge–Kutta explícito de 2 etapas de orden 3.

3. (Estabilidad lineal de RK4). Demuestra que para $\dot y = i\omega y$ el RK4 da $\lvert R(i h\omega)\rvert^2 = 1 - \tfrac{(h\omega)^6}{72} + \tfrac{(h\omega)^8}{576}$ y deduce el límite de estabilidad $h\omega \lt 2\sqrt2$.

4. (Estimador de error). Para un par encajado de órdenes $p$ y $\hat p = p - 1$, demuestra que $y_{n+1} - \hat y_{n+1}$ es una estimación asintóticamente correcta del error local de $\hat y_{n+1}$ pero no del de $y_{n+1}$, y explica por qué, aun así, se propaga $y_{n+1}$.

5. (Ley de control). Supón que $\mathrm{err}(h) = K h^{5}$ exactamente. Demuestra que, tras un rechazo, el nuevo paso $h' = 0.9\,h\,\mathrm{err}^{-1/5}$ se acepta siempre, y calcula cuánto desperdicia el factor 0.9 en el número de pasos.

6. (Computacional: tabla de orden). Implementa Euler, RK4 y Dormand–Prince de paso fijo y reproduce la tabla 5.1 para $e = 0.5$. Luego estima el orden empírico $p = \log_2\left(\mathrm{err}(h)/\mathrm{err}(h/2)\right)$ fila a fila y explica las desviaciones respecto de 1, 4 y 5.

import numpy as np
mu = 1.0
def f(t, y):
    r3 = np.hypot(y[0], y[1])**3
    return np.array([y[2], y[3], -mu*y[0]/r3, -mu*y[1]/r3])
def rk4(f, t, y, h):
    k1 = f(t, y); k2 = f(t+h/2, y+h/2*k1); k3 = f(t+h/2, y+h/2*k2); k4 = f(t+h, y+h*k3)
    return y + h/6*(k1 + 2*k2 + 2*k3 + k4)
def orbita(paso, n, e=0.5):
    y = np.array([1-e, 0, 0, np.sqrt((1+e)/(1-e))]); h = 2*np.pi/n
    for k in range(n): y = paso(f, k*h, y, h)
    return np.hypot(y[0]-(1-e), y[1])      # tras T la sonda vuelve al pericentro
for n in [100, 200, 400, 800, 1600]:
    print(n, orbita(rk4, n), orbita(rk4, n)/orbita(rk4, 2*n))

7. (Computacional: paso adaptativo). Implementa el controlador de §3.2 sobre Dormand–Prince e integra la órbita $e = 0.9$ una vuelta con $\mathrm{tol} = 10^{-6}, 10^{-8}, 10^{-10}$. Registra $h$ en función de la anomalía verdadera y comprueba que $h \propto r^{3/2}$ aproximadamente. ¿Por qué ese exponente?

Lecturas

• Hairer, Nørsett & Wanner (1993), Solving ODEs I: §II.1 (Runge–Kutta y teoría de Butcher), §II.3 (convergencia), §II.4 (estimación del error, control de paso, controlador PI de Gustafsson), §II.5 (tabla de Dormand–Prince).
• Press et al. (2007), Numerical Recipes, 3.ª ed.: §17.1–17.2 (StepperDopr5, control de paso).
• Dormand & Prince (1980), J. Comput. Appl. Math. 6, 19–26 (tabla 2: coeficientes RK5(4)7M). Fehlberg (1969), NASA TR R-315 (tabla III).
• Hairer, Lubich & Wanner (2006), Geometric Numerical Integration: §VIII.2 (por qué el paso variable destruye la simplecticidad y cómo recuperarla con transformaciones de tiempo).
• Murray & Dermott (1999), cap. 2: solución exacta del problema de dos cuerpos usada como referencia.