Integradores simplécticos y caos en el problema de N cuerpos

1. El problema de N cuerpos

Sean $N$ masas puntuales $m_i$ con posiciones $\mathbf r_i\in\mathbb R^3$ y velocidades $\mathbf v_i = \dot{\mathbf r}_i$, $i = 1,\dots,N$. La ley de gravitación de Newton da

donde $G$ es la constante de gravitación. En toda la lección (y en src/lib/math/nbody.ts) usamos unidades adimensionales con $G = 1$: fijada una masa de referencia $M$ y una longitud $R$, la unidad de tiempo es $\sqrt{R^3/(GM)}$.

Con momentos $\mathbf p_i = m_i\mathbf v_i$, el sistema es hamiltoniano con

y $\dot{\mathbf q} = \partial H/\partial\mathbf p$, $\dot{\mathbf p} = -\partial H/\partial\mathbf q$. Es separable: $T$ depende solo de $\mathbf p$ y $V$ solo de $\mathbf q$. Las diez integrales clásicas son la energía $E = H$, el momento lineal $\mathbf P = \sum_i\mathbf p_i$, el movimiento uniforme del centro de masas y el momento angular

Para $N\ge 3$ no hay integrales adicionales suficientes para resolver el problema por cuadraturas (resultados de Bruns y Poincaré), así que la integración numérica es inevitable. La pregunta es: ¿qué estructura del flujo exacto debe respetar un buen integrador?

2. Dos integradores

Escribimos $y = (\mathbf q,\mathbf p)$, $\dot y = f(y)$, y un paso de tamaño $h$ de $y_n$ a $y_{n+1}\approx y(t_n+h)$.

2.1 Runge–Kutta clásico (RK4)

Error local $O(h^5)$ y global $O(h^4)$ en un intervalo fijo. Pero RK4 no es simpléctico: no hay nada que impida que los errores de energía de cada paso se acumulen con el mismo signo. El resultado típico en órbitas acotadas es una deriva secular, $\lvert\Delta E\rvert \sim C\,t\,h^{p}$ (con $p\le 4$ según el problema; véase el ejercicio 3, donde para el oscilador armónico sale $t\,h^5$).

2.2 Leapfrog (Störmer–Verlet, kick-drift-kick)

Como $H = T + V$, los flujos exactos de $T$ y de $V$ por separado son triviales:

• flujo de $V$ (kick): $\varphi^V_\tau(\mathbf q,\mathbf p) = (\mathbf q,\ \mathbf p - \tau\,\nabla V(\mathbf q))$;
• flujo de $T$ (drift): $\varphi^T_\tau(\mathbf q,\mathbf p) = (\mathbf q + \tau M^{-1}\mathbf p,\ \mathbf p)$,

con $M = \operatorname{diag}(m_i)$. El leapfrog es la composición simétrica (splitting de Strang)

es decir, en términos de aceleraciones $\mathbf a(\mathbf r)$:

Cuesta una evaluación de fuerzas por paso (la última se reutiliza en el siguiente), frente a cuatro de RK4. Tiene tres propiedades estructurales:

1. Simpléctico. Cada factor es el flujo exacto de un hamiltoniano, luego es simpléctico, y la composición de aplicaciones simplécticas lo es (ejercicio 1). Recordemos que $\Phi$ es simpléctica si su jacobiana cumple $\Phi'(y)^{\mathsf T} J\,\Phi'(y) = J$ con $J = \begin{pmatrix}0 & I\\ -I & 0\end{pmatrix}$.
2. Simétrico (reversible): $\Phi_{-h} = \Phi_h^{-1}$; por eso su orden es par, aquí $2$.
3. Conserva $\mathbf L$ exactamente (salvo redondeo), porque cada kick aplica fuerzas centrales por pares y cada drift mueve $\mathbf r_i$ paralelo a $\mathbf v_i$ (ejercicio 2).

3. ¿Por qué la energía no deriva? El hamiltoniano modificado

El análisis de error regresivo (backward error analysis) busca un hamiltoniano $\tilde H$ cuyo flujo exacto coincida con el método numérico. Si $X_F$ denota el campo hamiltoniano de $F$, el flujo de $V$ durante $\tau$ es $\exp(\tau X_V)$ y la fórmula simétrica de Baker–Campbell–Hausdorff da

donde los términos de orden $h^3$ son conmutadores dobles de $X_V$ y $X_T$. Como el conmutador de campos hamiltonianos es (salvo signo) el campo del corchete de Poisson $\{F,G\} = \sum_i\left(\partial_{\mathbf q_i}F\cdot\partial_{\mathbf p_i}G - \partial_{\mathbf p_i}F\cdot\partial_{\mathbf q_i}G\right)$, el método es (formalmente) el flujo exacto de

Los coeficientes se comprueban en el oscilador armónico, donde todo es explícito (ejercicio 4). Al ser $\Phi_h$ simpléctico, la serie $\tilde H$ es hamiltoniana —esto es lo que falla para RK4— y la solución numérica conserva $\tilde H$ salvo términos exponencialmente pequeños. El resultado preciso (Hairer, Lubich & Wanner, cap. IX, teorema 8.1, para $H$ analítico) es

con $C$, $h_0$ independientes de $n$ y $h$. Es decir: error de energía acotado en $O(h^2)$ durante tiempos exponencialmente largos, sin deriva. El error oscila (crece en los pericentros y vuelve a bajar), pero no se acumula.

Dos advertencias honestas:

• La constante $C$ y la escala $h_0$ dependen de lo analítico que sea $H$ a lo largo de la órbita. En un encuentro cercano ($r_{ij}\to 0$) la escala de tiempo local $\sim\sqrt{r_{ij}^3/(G m)}$ se hace minúscula y un paso fijo deja de resolverla: el leapfrog de paso fijo pierde sus garantías. En la práctica se usan regularizaciones (Kustaanheimo–Stiefel), pasos adaptativos simétricos o integradores especializados.
• Un paso adaptativo ingenuo destruye la simplecticidad: la aplicación cambia de paso a paso con un $h(y)$.

3.1 Orden superior sin perder la estructura: el integrador de Yoshida

El leapfrog es de orden 2. Yoshida (1990) observó que, por ser simétrico ($\Phi_{-h} = \Phi_h^{-1}$), su desarrollo en $h$ solo contiene potencias impares, $\Phi_h = \exp\!\left(h X_H + h^3 X_3 + h^5 X_5 + \dots\right)$, y que una composición simétrica de tres pasos de tamaños $w_1 h$, $w_0 h$, $w_1 h$,

anula el término $h^3$: la primera condición garantiza la consistencia (avanzar $h$ en total) y la segunda cancela $X_3$ (los conmutadores de orden $h^3$ se suman linealmente en los $w_i^3$ al primer orden). La solución es

El paso central es negativo y mayor que $h$ en valor absoluto: el método va hacia delante, hacia atrás y de nuevo hacia delante. Es inevitable (Suzuki 1991): ninguna composición de pasos positivos del leapfrog supera el orden 2. El resultado es simpléctico y simétrico (composición de aplicaciones simplécticas y simétricas), de orden 4, con 3 evaluaciones de fuerza por paso. El mismo truco iterado da órdenes 6 y 8 (Yoshida da los coeficientes), con 7 y 15 leapfrogs por paso.

Está en src/lib/math/yoshida.ts. Medidas con nuestro código (test yoshida.test.ts):

Al dividir $h$ por 2 el error del leapfrog se divide por 4 y el de Yoshida por 16: órdenes 2 y 4. Y en la órbita excéntrica larga ambos quedan acotados (el máximo de la segunda mitad no supera al de la primera), mientras que RK4 con el mismo paso termina con $6.7\times10^{-3}$ y sigue creciendo. A igual coste (3 evaluaciones de Yoshida frente a 3 pasos de leapfrog de tamaño $h/3$, cuyo error sería $2.2\times10^{-2}/9 = 2.4\times10^{-3}$), Yoshida gana por un factor 4 en este problema, y la ventaja crece al reducir $h$.

3.2 El teorema de conservación de energía en general

El resultado de la sección anterior no es exclusivo del orden 2. Para cualquier integrador simpléctico de orden $p$ aplicado a un hamiltoniano analítico, con paso suficientemente pequeño (Hairer, Lubich & Wanner, cap. IX, teorema 8.1),

La demostración sigue el mismo camino: el método es el flujo exacto, salvo términos exponencialmente pequeños, de un hamiltoniano modificado $\tilde H = H + h^p H_{p+1} + h^{p+1}H_{p+2} + \dots$ (truncado óptimamente), y $\tilde H$ se conserva exactamente por ese flujo; la diferencia $\tilde H - H = O(h^p)$ es acotada en una región acotada del espacio de fases. Para Yoshida, $p = 4$: eso es lo que muestra la tabla. La condición «paso suficientemente pequeño» es seria: Yoshida tiene pasos internos de tamaño $1.7h$, así que su límite de estabilidad en el oscilador es $h\omega \lt 2/1.7024 = 1.17$ en vez de $2$ (ejercicio 7).

Qué observar. Con el escenario «2 cuerpos, e = 0.8» y «Ambos superpuestos», el error del leapfrog (línea continua) oscila con picos en cada pericentro pero su envolvente no crece; el de RK4 (discontinua) arranca más pequeño en muchos tramos y crece de forma monótona hasta superarlo. Sube la velocidad para ver decenas de periodos. En los tests de src/lib/math/__tests__/nbody.test.ts (e = 0.7, $h = 0.03$, 100 periodos) el error de RK4 crece linealmente con $t$, mientras que el máximo del leapfrog en la segunda mitad no supera al de la primera. En nuestras corridas la deriva de RK4 escala aproximadamente como $h^5$ (paso $0.02\to0.03$ multiplica el error por $\approx 7.5 \approx 1.5^5$), en línea con el ejercicio 3.

4. Sensibilidad a las condiciones iniciales

Sea $y(t)$ una solución y $y(t) + \delta y(t)$ otra que parte a distancia $\lVert\delta y(0)\rVert = \delta_0$ en el espacio de fases. El máximo exponente de Lyapunov es

Si $\lambda > 0$ la separación crece como $\delta_0 e^{\lambda t}$ y la órbita es caótica; $t_L = 1/\lambda$ es el tiempo de Lyapunov. En un sistema integrable (Kepler) las órbitas vecinas se separan solo polinómicamente (por la diferencia de periodos) y $\lambda = 0$.

Numéricamente la separación satura cuando deja de ser pequeña, así que el método de Benettin et al. (1980) renormaliza: cada intervalo $\tau$ mide $\delta_k$, acumula $\ln(\delta_k/\delta_0)$ y reescala la trayectoria vecina a distancia $\delta_0$ a lo largo de la dirección actual. Tras $K$ intervalos,

Está implementado en lyapunovBenettin; el simulador muestra además la estimación sin renormalizar $\ln(\delta/\delta_0)/t$ (válida mientras $\delta \ll 1$). La distancia usada es la euclídea en $\mathbb R^{6N}$, $\delta^2 = \sum_i \lVert\mathbf r_i-\mathbf r_i'\rVert^2 + \lVert\mathbf v_i-\mathbf v_i'\rVert^2$.

Cómo medirlo en la práctica. Tres decisiones afectan al número que sale: (1) el intervalo de renormalización $\tau$ debe ser menor que el tiempo en que $\delta$ sale del régimen lineal ($\delta_0 e^{\lambda\tau} \ll 1$) pero no tan corto que domine el ruido de redondeo; (2) $\delta_0$ debe superar con holgura el redondeo ($10^{-8}$ frente a $10^{-16}$) para que la dirección inicial de la perturbación sea significativa; (3) $\lambda$ se define como límite $t \to \infty$, pero la estimación a tiempo finito converge como $\sim 1/t$ (ley de los grandes números sobre los $\ln(\delta_k/\delta_0)$) y, en sistemas que alternan tramos regulares y encuentros cercanos, fluctúa fuertemente de un intervalo a otro. Por eso debe darse siempre con su tiempo de integración y, si es posible, promediada sobre condiciones iniciales vecinas. Una comprobación útil: en un sistema integrable (dos cuerpos) el estimador debe decaer como $\ln t/t$ hacia cero; si no lo hace, el error numérico está siendo tomado por caos. Y otra: $\lambda$ no depende del integrador si $h$ es pequeño; si cambia al cambiar $h$ (como ocurre en el problema pitagórico con $h = 10^{-3}$ frente a $10^{-4}$), los encuentros cercanos no están resueltos.

Consecuencia práctica. Si queremos el estado con tolerancia $\Delta$ partiendo de una incertidumbre $\delta_0$, el horizonte de predicción es

que solo crece logarítmicamente con la precisión inicial. Ni siquiera el redondeo de doble precisión ($\delta_0\sim 10^{-16}$) escapa: también es una perturbación que el caos amplifica. Por eso en un sistema caótico no tiene sentido pedir la trayectoria exacta a tiempos largos, sino propiedades estadísticas o estructurales (energía, estado final), y por eso importa tanto que el integrador respete las cantidades conservadas.

5. Dos soluciones paradigmáticas

5.1 La figura en ocho

Moore (1993) encontró numéricamente una órbita periódica en la que tres masas iguales se persiguen sobre una misma curva en forma de ocho; Chenciner y Montgomery (2000) demostraron su existencia minimizando la acción $\int_0^T (T - V)\,dt$ sobre lazos con las simetrías del ocho. Con $m_i = 1$, $G = 1$ y las condiciones de Simó (2002),

el periodo es $T \approx 6.3259$, el momento angular es nulo y la energía vale $E \approx -1.28714$ (comprobado en los tests). Simó mostró que es linealmente estable: una perturbación de $10^{-6}$ no destruye el ocho, y en el simulador $\delta(t)$ crece despacio, no exponencialmente. Es una buena prueba de integradores: tras un periodo el estado debe volver al inicial (en los tests, a menos de $10^{-4}$ con $h = T/6000$).

5.2 El problema pitagórico de Burrau

Burrau (1913) planteó tres masas $3$, $4$ y $5$ en reposo en los vértices de un triángulo rectángulo de lados $3$, $4$, $5$, cada masa frente al lado de igual longitud. Con las coordenadas de Szebehely & Peters (1967),

el centro de masas está en el origen, $\mathbf L = \mathbf 0$ y $E = V = -\left(\tfrac{12}{5} + \tfrac{15}{4} + \tfrac{20}{3}\right) = -\tfrac{769}{60}$. Burrau solo pudo integrar a mano el comienzo; Szebehely y Peters (1967), ya con ordenador y regularización, mostraron que tras una serie de encuentros cercanos el sistema se rompe: hacia $t \approx 59.4$ (unidades $G = 1$) se produce un último encuentro triple y la masa $3$ es expulsada definitivamente, mientras las masas $4$ y $5$ quedan como binaria; Peters prolongó la integración hasta $t = 67$ para confirmar la escapada (movimiento «elíptico-hiperbólico»).

Aquí aparecen juntas las dos lecciones: los encuentros cercanos hacen que el leapfrog de paso fijo muestre picos de $\lvert\Delta E/E_0\rvert$ (prueba a variar $h$) y la dinámica es caótica. Con nuestra implementación (lyapunovBenettin, $h = 10^{-4}$, $t\in[0,20]$) la estimación a tiempo finito es $\lambda \approx 1$ en unidades $G = 1$, frente a $\approx 0.1$ para la órbita circular de dos cuerpos en $t\in[0,50]$ (en este caso el valor verdadero es $0$ y lo que se mide es el decaimiento lento $\sim\ln t/t$ de la estimación). Son estimaciones de nuestro código, no valores de referencia.

Problemas

Ejercicio 1 (simplecticidad del leapfrog). Demuestra que $\varphi^V_\tau$ y $\varphi^T_\tau$ son simplécticas y concluye que $\Phi_h$ lo es.

Ejercicio 2 (conservación exacta de $\mathbf L$). Demuestra que el leapfrog conserva $\mathbf L$ para cualquier $h$.

Ejercicio 3 (RK4 disipa energía en el oscilador). Para $\dot q = p$, $\dot p = -q$ escribe $z = q + i p$, de modo que $\dot z = -i z$. Demuestra que RK4 da $z_{n+1} = R(-ih)\,z_n$ con $R(w) = 1 + w + w^2/2 + w^3/6 + w^4/24$ y que

Ejercicio 4 (invariante exacto del leapfrog en el oscilador). Demuestra que el leapfrog aplicado a $H = (p^2 + q^2)/2$ es la aplicación lineal

que conserva exactamente $I_h = p^2 + \left(1 - \tfrac{h^2}4\right)q^2$, que es estable si y solo si $0 \lt h \lt 2$, y que esto concuerda con $\tilde H$ de la sección 3.

Ejercicio 5 (horizonte de predictibilidad). Con $\lambda = 1$ (unidades $G=1$), ¿cuánto se alarga el horizonte de predicción con tolerancia $\Delta = 10^{-2}$ si pasamos de $\delta_0 = 10^{-6}$ a $\delta_0 = 10^{-16}$?

Ejercicio 6 (con el simulador). (a) En el ocho con «Ambos», fija $h = 0.1$: ¿qué integrador mantiene la forma de la órbita durante más periodos? Relaciónalo con la sección 3. (b) Activa «Perturbar» en el ocho y en el problema pitagórico. Compara cómo crece $\delta(t)$ en la gráfica logarítmica: una recta de pendiente $\lambda/\ln 10$ indica crecimiento exponencial. (c) En el pitagórico, reduce $h$ de $10^{-2}$ a $10^{-4}$. ¿Cambia el estado final? ¿Por qué no contradice esto el teorema de la sección 3? (Pista: encuentros cercanos y caos.)

Ejercicio 7 (computacional: Yoshida). Implementa el leapfrog y el integrador de Yoshida de orden 4 para el problema de dos cuerpos ($G = 1$, $m_1 = 1$, $m_2 = 0.01$, $a = 1$, $e = 0.7$). (a) Reproduce la tabla de §3.1 midiendo $\max\lvert\Delta E/E_0\rvert$ en 100 periodos con $h = 0.03$. (b) Mide el orden de ambos en un periodo con $h = 0.02, 0.01, 0.005$. (c) Aplica Yoshida al oscilador armónico $\dot q = p$, $\dot p = -q$ y determina numéricamente el $h$ a partir del cual la solución explota; compáralo con $2/1.7024$.

import numpy as np
w1 = 1/(2 - 2**(1/3)); w0 = -2**(1/3)/(2 - 2**(1/3))
def acel(r):                      # r: (N,3); masas m globales
    a = np.zeros_like(r)
    for i in range(len(r)):
        d = r - r[i]; dist = np.linalg.norm(d, axis=1); dist[i] = np.inf
        a[i] = np.sum((m/dist**3)[:, None]*d, axis=0)
    return a
def leapfrog(r, v, h):
    v = v + 0.5*h*acel(r); r = r + h*v; v = v + 0.5*h*acel(r); return r, v
def yoshida4(r, v, h):
    for w in (w1, w0, w1): r, v = leapfrog(r, v, w*h)
    return r, v
def energia(r, v):
    K = 0.5*np.sum(m*np.sum(v**2, axis=1)); U = 0
    for i in range(len(r)):
        for j in range(i+1, len(r)): U -= m[i]*m[j]/np.linalg.norm(r[i]-r[j])
    return K + U

Lecturas

• Hairer, Lubich & Wanner (2006), Geometric Numerical Integration: cap. I (ejemplos y leapfrog), §II.4–II.5 (métodos de composición y splitting: construcción de Yoshida), cap. VI (integradores simplécticos), cap. IX (análisis de error regresivo; teorema 8.1 de conservación de energía).
• Yoshida (1990), Phys. Lett. A 150, 262: los coeficientes de órdenes 4, 6 y 8 y la demostración por BCH.
• Chenciner & Montgomery (2000), Ann. Math. 152, 881, y Simó (2002): existencia y estabilidad del ocho. Moore (1993), PRL 70, 3675: su descubrimiento numérico.
• Szebehely & Peters (1967), AJ 72, 876: la solución completa del problema pitagórico, con las figuras de la expulsión.
• Benettin, Galgani, Giorgilli & Strelcyn (1980), Meccanica 15, 9: el método de renormalización para los exponentes de Lyapunov.
• Aarseth, S. J. (2003), Gravitational N-Body Simulations, cap. 1–2: el punto de vista del especialista en N cuerpos sobre integradores y regularización.