Integradores simplécticos y errores a largo plazo

Año 2 · Módulo 12 · Nivel T

119.1. La diferencia entre calcular una trayectoria y representar una dinámica

Una órbita numérica puede aproximar muy bien la solución verdadera durante los primeros períodos y acumular después un error energético sistemático. También puede conservar la energía con gran precisión y describir de forma equivocada el desplazamiento de fase. En simulaciones prolongadas de sistemas hamiltonianos, estos errores tienen consecuencias distintas. Por ello, la precisión posicional inmediata debe evaluarse junto con las propiedades geométricas que conserva el integrador.

El error local mide lo que pierde una aproximación en un paso. El error global compara trayectorias durante un intervalo fijado. La conservación de energía, momento y estructura simpléctica exige pruebas adicionales. Una trayectoria muy larga permite observar propiedades dinámicas, pero también amplifica los defectos de la discretización; la convergencia al reducir el paso no reemplaza la comprobación de invariantes ni viceversa.

Los integradores simplécticos están construidos para conservar la forma simpléctica del flujo en cada paso, bajo las condiciones pertinentes. Esa estructura permite explicar mediante análisis retrospectivo por qué sus errores de energía suelen permanecer acotados durante tiempos largos, aunque la precisión de la trayectoria puntual siga deteriorándose. Compararemos métodos concretos y examinaremos qué tipo de error importa según la pregunta física que motiva la simulación.

La integración geométrica selecciona algunas propiedades estructurales del problema y construye una evolución discreta compatible con ellas, aunque las demás magnitudes tengan error en cada paso. Veremos por qué esta decisión suele ser preferible para horizontes largos, en qué condiciones permite controlar el error de energía y cuándo una alternativa no simpléctica de orden alto sigue siendo la elección más razonable.

119.2. La estructura que intenta respetar el integrador

Para un hamiltoniano autónomo regular en coordenadas canónicas $z=(q,p)\in\mathbb R^{2n}$ , las ecuaciones de Hamilton son

\dot q=\nabla_pH(q,p),\qquad \dot p=-\nabla_qH(q,p).

Introduciendo la matriz

J=\begin{pmatrix}0&I\\-I&0\end{pmatrix},\qquad \dot z=J\nabla H(z),

expresamos en una sola ecuación todas las parejas posición-momento. Para el flujo exacto $\Phi^t$ , definido en una región donde las soluciones existen, su derivada satisface

(D\Phi^t)^T J(D\Phi^t)=J.

Es la condición simpléctica en coordenadas canónicas. Implica la preservación de la forma diferencial $\omega=\sum_i dq_i\wedge dp_i$ y, en consecuencia, del volumen canónico. La implicación inversa no vale en general si hay más de un grado de libertad, porque conservar el determinante no garantiza que se preserven todas las áreas orientadas en los planos canónicos.

Un método de un paso define una aplicación $\Psi_h:z_n\mapsto z_{n+1}$ . Es simpléctico si, para cada punto donde es diferenciable,

\boxed{(D\Psi_h)^T J(D\Psi_h)=J.}

Es una propiedad del mapa completo, compatible con que las aproximaciones de $q$ y de $p$ tengan error por separado. La composición de dos mapas simplécticos también es simpléctica. Por eso los métodos de composición pueden resolver subflujos simples y concatenarlos sin perder la estructura simpléctica, siempre que cada subflujo la conserve.

La preservación de la forma simpléctica y la conservación exacta de la energía son propiedades diferentes del esquema numérico. En general, un integrador simpléctico no mantiene fijo el valor de $H$ ; la estructura que preserva produce, bajo hipótesis adicionales, un comportamiento energético característico. Tampoco asegura la conservación exacta del momento angular, que en determinados métodos se obtiene cuando la discretización respeta la simetría rotacional apropiada.

Las propiedades de los mapas dependen además de sus variables. Si formulamos el problema en coordenadas no canónicas, la matriz $J$ constante puede dejar de describir la estructura correcta. Que el modelo tenga una energía no hace simpléctico al algoritmo de integración. El cambio de coordenadas debe tratarse mediante la forma simpléctica transformada o mediante una formulación canónica consistente.

119.3. Qué hace mal Euler explícito en el oscilador más sencillo

Tomemos como experimento analítico el oscilador armónico de masa unitaria,

H(q,p)=\frac12(p^2+\omega^2q^2),\qquad \dot q=p,\quad \dot p=-\omega^2q.

El flujo exacto describe elipses cerradas en el plano $(q,p)$ , sin pérdida de amplitud ni aumento de energía. Con paso $h$ , Euler explícito da

q_{n+1}=q_n+hp_n,\qquad p_{n+1}=p_n-h\omega^2q_n.

Su matriz de transición es

M_E=\begin{pmatrix}1&h\\-h\omega^2&1\end{pmatrix}, \qquad\det M_E=1+h^2\omega^2>1.

Para un sistema canónico de un grado de libertad, esta expansión de área demuestra que el método no es simpléctico. Un cálculo directo de la energía da

H_{n+1}=\frac12\left[(p_n-h\omega^2q_n)^2+ \omega^2(q_n+hp_n)^2\right] =(1+h^2\omega^2)H_n.

Así, la energía numérica se multiplica exactamente por un factor mayor que uno en cada iteración. Tras $N$ pasos, $H_N=(1+h^2\omega^2)^NH_0$ . Cuando el tiempo total $Nh$ es grande, la simulación dibuja una espiral hacia fuera, pese a que las soluciones exactas de la ecuación están confinadas a elipses. Esa inestabilidad la introduce el algoritmo; el sistema físico no la tiene.

El defecto persiste aunque $h$ sea pequeño si mantenemos el paso fijo y prolongamos la integración. Al disminuir $h$ el crecimiento artificial es más lento, pero sigue acumulándose. Una gráfica de pocos períodos puede ocultar esa deriva.

Euler implícito, aplicado al mismo modelo, produce el comportamiento cualitativo opuesto. Su matriz tiene determinante $1/(1+h^2\omega^2)$ y contrae áreas, de modo que aparece una disipación artificial. Ambos métodos convergen al flujo exacto cuando $h\to0$ para tiempos finitos, pero a paso fijo y tiempos largos perturban el mecanismo dinámico de maneras incompatibles con su geometría.

119.4. Dividir el hamiltoniano para obtener mapas exactos

Supongamos ahora un hamiltoniano separable,

H(q,p)=T(p)+V(q),\qquad T(p)=\frac12p^TM^{-1}p,

con $M$ simétrica definida positiva. El flujo generado por $T$ mantiene constante $p$ y actualiza las posiciones según

D_h:(q,p)\longmapsto(q+hM^{-1}p,p).

Lo llamaremos deriva (drift). El flujo generado por $V$ mantiene fija la posición y cambia los momentos según

K_h:(q,p)\longmapsto(q,p-h\nabla V(q)).

Es el impulso (kick). Cada transformación es un flujo hamiltoniano exacto y, en consecuencia, es simpléctica. Aunque $T$ y $V$ no evolucionen de manera independiente en el sistema real, las composiciones $D_h\circ K_h$ y $K_h\circ D_h$ son integradores de primer orden que conservan la estructura simpléctica.

El primer orden de composición actualiza el momento antes que la posición,

p_{n+1}=p_n-h\nabla V(q_n),\qquad q_{n+1}=q_n+hM^{-1}p_{n+1}.

Es una de las variantes de Euler simpléctico. El orden contrario define otra variante y genera un error de fase diferente. Que ambos mapas sean simplécticos no los hace idénticos ni de segundo orden; la diferencia se debe a que los flujos no conmutan. Si intercambiamos las etapas sin analizar su composición, cambiamos el algoritmo y su hamiltoniano modificado.

En el oscilador, la matriz del esquema impulso-deriva es

M_S=\begin{pmatrix}1-h^2\omega^2&h\\-h\omega^2&1\end{pmatrix}, \qquad\det M_S=1.

En dos dimensiones, esta condición equivale a la simplecticidad. Sus autovalores están sobre el círculo unitario cuando $0<h\omega<2$ , condición que garantiza la estabilidad lineal del esquema. Para pasos mayores el área se sigue conservando, pero aparece inestabilidad numérica; la geometría simpléctica no reemplaza la elección de un paso adecuado.

El método sigue una dinámica hamiltoniana ligeramente distinta de la exacta, sin fuerza viscosa numérica. La energía original puede oscilar alrededor del valor correcto sin el crecimiento secular que tenía Euler explícito en este ejemplo. Una composición simétrica aumenta el orden local y mantiene la preservación de la forma simpléctica.

119.5. Simetría temporal y segundo orden en el método de Störmer–Verlet

Una composición muy usada es

\Psi_h=K_{h/2}\circ D_h\circ K_{h/2},

entendiendo que las operaciones se ejecutan de derecha a izquierda. Su forma operacional, conocida como velocity Verlet o impulso-deriva-impulso, es

\begin{aligned} p_{n+1/2}&=p_n-\frac h2\nabla V(q_n),\\ q_{n+1}&=q_n+hM^{-1}p_{n+1/2},\\ p_{n+1}&=p_{n+1/2}-\frac h2\nabla V(q_{n+1}). \end{aligned}

Cada subetapa es simpléctica y, por lo tanto, lo es la composición. El esquema es además simétrico en el tiempo, es decir, al reemplazar $h$ por $-h$ e invertir el orden de ejecución recuperamos el estado anterior, suponiendo aritmética exacta. Es una propiedad distinta de la conservación de energía y tiene consecuencias fuertes sobre la estructura del error. La composición simétrica cancela los términos pares pertinentes de la expansión local del operador; de ahí el orden global dos para problemas suficientemente regulares.

El mismo método puede escribirse como deriva-impulso-deriva, $D_{h/2}\circ K_h\circ D_{h/2}$ . Son esquemas emparentados y de igual orden, pero no necesariamente entregan idénticos estados al mismo tiempo si las variables se colocan en instantes diferentes. Para comparar implementaciones hay que declarar esa convención de colocación temporal.

La estructura de pasos de Verlet hace económico su uso en dinámica molecular, sistemas de partículas y modelos gravitacionales. Para una fuerza que depende solo de las posiciones, cada paso completo necesita una sola evaluación nueva de fuerza, porque reutiliza la evaluación previa. Eso suele ser mucho más barato que un esquema Runge–Kutta con varias evaluaciones por paso. A esa economía se suma una ventaja geométrica, pues en regímenes donde la teoría de error retrospectivo es aplicable las trayectorias de largo plazo suelen mostrar errores energéticos oscilatorios y acotados.

De ello no se sigue que Verlet sea más preciso en cualquier tarea. Si necesitamos determinar la posición de un cometa a una fecha próxima con un margen ínfimo, un algoritmo adaptativo de orden alto puede ser más eficiente. Si aparecen choques, discontinuidades o pasos muy diferentes entre encuentros cercanos, Verlet con paso fijo puede resultar peligroso. Con fuerzas disipativas, la física deja de ser hamiltoniana en las variables ordinarias y exigir simplecticidad sin reformular el problema carece de sentido.

Figura 119.1. C119 Figura 01 Euler Verlet energia trayectorias Ampliar figura ↗

Figura 119.1. Comparación para el oscilador armónico $H=(q^2+p^2)/2$ , con $q(0)=1$ , $p(0)=0$ y paso fijo $h=0{,}15$ . El panel A contrapone la expansión espuria de Euler explícito con las trayectorias acotadas de Euler simpléctico y Verlet; el recuadro amplía los movimientos cercanos a la órbita exacta. El panel B muestra la energía: en Euler explícito, $H_n=H_0(1+h^2)^n$ , mientras los dos integradores simplécticos mantienen un error energético acotado en la ventana presentada. La preservación simpléctica no significa conservación exacta de $H$ .

119.6. Una energía exactamente conservada que no coincide con la original

En el oscilador podemos ir más allá de afirmar que Verlet «funciona bien» e identificar una cantidad cuadrática que su mapa preserva exactamente. Para $T=p^2/2$ y $V=\omega^2q^2/2$ , las ecuaciones del método dan

\begin{pmatrix}q_{n+1}\\p_{n+1}\end{pmatrix} =\begin{pmatrix} 1-\frac{h^2\omega^2}{2}&h\\ -h\omega^2\left(1-\frac{h^2\omega^2}{4}\right)& 1-\frac{h^2\omega^2}{2} \end{pmatrix} \begin{pmatrix}q_n\\p_n\end{pmatrix}.

Llamemos $a=1-h^2\omega^2/2$ y $b=\omega^2(1-h^2\omega^2/4)$ . La matriz tiene la forma $\left(\begin{smallmatrix}a&h\\-hb&a\end{smallmatrix}\right)$ y satisface $a^2+h^2b=1$ . De ahí se verifica por multiplicación que

\boxed{E_h(q,p)=\tfrac12\bigl(p^2+bq^2\bigr)}

es exactamente invariante a lo largo de las iteraciones. Para $0<h\omega<2$ se tiene $b>0$ , de modo que las trayectorias numéricas viven en elipses cerradas de esta energía modificada. El hamiltoniano físico original, $H=\tfrac12(p^2+\omega^2q^2)$ , difiere de $E_h$ en $h^2\omega^4q^2/8$ . Por ello su valor oscila, pero la amplitud de la diferencia se controla mientras la órbita permanece acotada.

El criterio de estabilidad también se deduce de la traza, ya que los autovalores cumplen $\lambda+\lambda^{-1}=2a$ . Cuando $|a|<1$ son complejos conjugados de módulo uno; cuando $a<-1$ uno de ellos crece en módulo. En la frontera $h\omega=2$ aparecen autovalores $-1$ y, genéricamente, una matriz no diagonalizable que permite crecimiento lineal de amplitud. Haber conservado el determinante uno no previene este fallo.

Queda la fase. Para un oscilador numéricamente estable, el ángulo avanzado por paso satisface

\cos\theta=1-\frac{h^2\omega^2}{2},\qquad \omega_{\rm num}=\frac{\theta}{h}= \frac2h\arcsin\left(\frac{h\omega}{2}\right).

Al expandir la función inversa encontramos $\omega_{\rm num}=\omega+h^2\omega^3/24+O(h^4)$ . La frecuencia calculada está ligeramente desplazada, de manera que el error de fase acumulado aumenta aproximadamente como $t h^2$ en ese régimen. Una conservación energética excelente es compatible con una posición orbital cada vez más desfasada, así que un error energético pequeño no prueba precisión orbital.

119.7. El análisis retrospectivo y la ecuación que parece resolver el cálculo

El análisis de error retrospectivo (backward error analysis) invierte la pregunta habitual sobre cuánto se aparta un paso de la trayectoria exacta del problema original y busca una ecuación diferencial modificada cuyo flujo reproduzca, hasta un orden dado, el mapa discreto. Para un método simpléctico y suficientemente regular, la ecuación modificada puede escribirse formalmente como un sistema hamiltoniano con función

\widetilde H=H+h^rH_r+h^{r+1}H_{r+1}+\cdots,

si el método tiene orden $r$ y se cumplen las condiciones que justifican la expansión. Para un método simétrico de orden dos aparecen, en la expansión pertinente, potencias pares de $h$ . El hamiltoniano modificado explica por qué el algoritmo puede conservar casi exactamente una energía cercana a la física incluso cuando la trayectoria exacta ya no coincide punto por punto.

En el esquema impulso-deriva-impulso para $H=p^2/2+V(q)$ , con masas unitarias y potencial suave, los primeros términos son

\widetilde H(q,p)=H(q,p)+h^2\left[ \frac1{12}p^T\nabla^2V(q)p- \frac1{24}\lVert\nabla V(q)\rVert^2\right]+O(h^4).

Esta serie no tiene por qué converger para cualquier $h$ ni definir globalmente una energía exacta; habitualmente es una expansión asintótica, válida dentro de un dominio regular, que se trunca de forma óptima para obtener estimaciones útiles. En sistemas analíticos con soluciones confinadas en conjuntos compactos apropiados y con restricciones sobre el paso, los resultados clásicos muestran que el defecto de la energía original puede permanecer del orden del error del método durante tiempos que crecen exponencialmente al disminuir $h$ , antes de que otras limitaciones entren en juego.

También importan el tamaño de las constantes, la aparición de frecuencias altas, la regularidad y la proximidad de singularidades. La interpretación mediante un hamiltoniano modificado tiene un dominio perturbativo y temporal limitado; los teoremas de error retrospectivo no prometen conservación exacta de la energía original para tiempos arbitrarios. Los errores de redondeo pueden acumularse, un algoritmo mal implementado puede introducir sesgos y las resonancias numéricas o los cambios bruscos de régimen pueden invalidar el uso simple de una expansión local.

En problemas caóticos, dos órbitas con estados iniciales vecinos divergen exponencialmente durante ciertos intervalos, incluso si el integrador reproduce bien cantidades estadísticas o geométricas. El análisis retrospectivo no suprime la sensibilidad inicial, pero permite separar la desviación por truncamiento, la amplificación dinámica de una desviación y la alteración de la estructura del sistema por la discretización.

Figura 119.2. C119 Figura 02 Energia modificada desfase Ampliar figura ↗

Figura 119.2. Oscilador armónico integrado con Störmer–Verlet, $h=0{,}2$ , $q_0=1$ y $p_0=0$ . En el panel A, la energía física oscila alrededor de su valor inicial mientras $E_h=[p^2+(1-h^2/4)q^2]/2$ permanece constante salvo redondeo. En el panel B, el error de fase aumenta de forma aproximadamente lineal y coincide con la relación exacta $\omega_{num}=2\arcsin(h/2)/h$ . A $t=1000$ el desfase es de 1,6742 radianes. Conservar el invariante modificado no garantiza precisión de fase a largo plazo.

119.8. Cuatro preguntas sobre error global, error de fase y conservación

La evaluación de una integración larga depende del observable que se pretende aproximar. La primera pregunta es geométrica y consiste en saber si el integrador preserva la estructura simpléctica, una simetría o un volumen relevante. La segunda es mecánica (¿cómo fluctúan la energía y los momentos integrales?). La tercera es orbital y mide la diferencia en posición y fase con respecto a una solución de referencia. La cuarta, estadística, examina si las frecuencias de visita, las distribuciones o las tasas de transporte cambian al refinar el paso. La respuesta a una de ellas deja abiertas las otras tres.

Para una órbita de Kepler, por ejemplo, conservar aproximadamente energía y momento angular ayuda a mantener casi constantes el semieje mayor y la excentricidad. Sin embargo, una pequeña desviación en la frecuencia media puede desplazar la anomalía orbital después de muchas vueltas. Si solo comparamos posiciones en un instante muy tardío, podríamos diagnosticar una gran divergencia, aunque el método haya seguido con fidelidad una familia de elipses cercanas. Que ese desfase sea o no el error decisivo depende de si queremos una efeméride o un retrato estadístico de largo plazo.

Una medida útil de energía es $\max_{0\le n\le N}|H(z_n)-H(z_0)|$ , acompañada de su evolución temporal y de un contraste entre distintos $h$ . El valor final aislado es insuficiente porque una oscilación puede regresar accidentalmente cerca de cero. Cuando el modelo conserva momento angular, se compara además $L_n$ con $L_0$ . Para detectar defectos de fase en movimientos casi periódicos podemos registrar los tiempos de paso por una sección. Para investigar caos o difusión interesa añadir exponentes finitos de Lyapunov, recorridos de acciones, histogramas y análisis de convergencia de cada observable.

Un criterio imprescindible es la convergencia de las conclusiones. Si una sección de Poincaré muestra un mar caótico con paso $h$ , pero al emplear $h/2$ ese patrón desaparece, aún no tenemos evidencia fiable de caos físico. Del mismo modo, una energía aparentemente estable puede ocultar una resonancia espuria que modifica la tasa de transporte. La integración estructural reduce algunos riesgos, pero sigue siendo necesario investigar la sensibilidad numérica.

119.9. Pasos adaptativos, encuentros cercanos y límites operacionales

Los integradores simplécticos de composición tradicionales tienen la ventaja de una estructura fija. Sin embargo, los problemas reales mezclan escalas temporales diferentes. Un planeta puede pasar largos intervalos lejos de su estrella y acelerar mucho cerca del periastro; dos cuerpos pueden aproximarse y después separarse; un sistema molecular puede combinar vibraciones rápidas con movimientos lentos. Parece natural disminuir el paso solo cuando el problema se vuelve exigente.

La adaptación ingenua del paso según el estado suele romper la estructura simpléctica de los mapas canónicos ordinarios. Aunque cada transformación $\Psi_h$ sea simpléctica para un valor fijo de $h$ , el mapa $z\mapsto\Psi_{h(z)}(z)$ incorpora derivadas adicionales de $h(z)$ y no tiene por qué satisfacer $D\Psi^TJD\Psi=J$ . La adaptación puede hacerse compatible con la integración geométrica (hay formulaciones de tiempo extendido, transformaciones de Sundman y métodos especializados), pero exige más que reemplazar una constante por una función y suponer que la garantía anterior sigue valiendo.

Los encuentros cercanos introducen otra dificultad. En gravitación newtoniana, $V=-\mu/r$ y sus derivadas aumentan sin cota al aproximarnos a una colisión. Las hipótesis de suavidad y acotación en un dominio compacto alejado de singularidades dejan de cumplirse. Un paso fijo adecuado para una órbita casi circular puede atravesar la región periástrica de una órbita muy excéntrica sin resolverla. Un pequeño error energético durante el resto de la trayectoria no justifica ese salto.

En astronomía computacional se usan integradores de escisión que resuelven exactamente una parte kepleriana e incorporan interacciones más débiles como impulsos. Los esquemas de Wisdom–Holman se basan precisamente en esa separación. Cuando aparecen encuentros cercanos o fuerzas no conservativas, pueden ser preferibles algoritmos adaptativos de orden alto, métodos regularizados o estrategias híbridas. Los métodos simplécticos no son la mejor opción en todos los escenarios; la elección depende del hamiltoniano, del objetivo observacional, de las escalas y del presupuesto de error.

También importa la aritmética de coma flotante. Las sumas sucesivas, los cambios repetidos de coordenadas y las operaciones que restan números casi iguales pueden introducir sesgos. Un algoritmo geométricamente correcto en aritmética exacta no necesariamente conserva sus virtudes sin límite en un ordenador. La implementación es parte del experimento científico.

119.10. Un experimento reproducible con dos mapas y una energía modificada

El primer ensayo usa el oscilador, cuyos resultados pueden cotejarse con identidades exactas. Elegimos $\omega=1$ , $q_0=1$ , $p_0=0$ , paso $h=0,15$ y $40\,000$ iteraciones. El experimento calcula Euler explícito y Verlet, pero detiene la parte de Euler después de un intervalo moderado para evitar que el crecimiento exponencial desborde la aritmética. Luego mide, para Verlet, la energía física y la energía modificada definida en la sección 119.6.

import numpy as np

omega, h = 1.0, 0.15
q, p = 1.0, 0.0
H0 = 0.5 * (p*p + omega**2*q*q)
Eh0 = 0.5 * (p*p + omega**2*(1-(h*omega)**2/4)*q*q)
error_H = 0.0
error_Eh = 0.0

for _ in range(40_000):
    pmedio = p - 0.5*h*omega**2*q
    qnuevo = q + h*pmedio
    pnuevo = pmedio - 0.5*h*omega**2*qnuevo
    q, p = qnuevo, pnuevo
    H = 0.5*(p*p + omega**2*q*q)
    Eh = 0.5*(p*p + omega**2*(1-(h*omega)**2/4)*q*q)
    error_H = max(error_H, abs(H-H0))
    error_Eh = max(error_Eh, abs(Eh-Eh0))

print('Máximo defecto en H:', error_H)
print('Máximo defecto en Eh:', error_Eh)
print('Frecuencia numérica:', 2*np.arcsin(h*omega/2)/h)

El código debe confirmar que el defecto de la energía modificada es comparable al redondeo acumulado, mientras que el error de la energía física es pequeño, acotado y mayor. No hemos supuesto igualdad de los errores ni asignado cifras que dependan de la plataforma. El valor de la frecuencia calculada se contrasta directamente con la expresión analítica.

En un segundo ensayo, conceptualmente independiente, podemos integrar una órbita de Kepler mediante Verlet y un método adaptativo de referencia. Para $\mu=1$ , $a=1$ , excentricidad $e=0,3$ , partimos del periastro con $r_0=(1-e,0)$ y velocidad transversal $v_0=(0,\sqrt{(1+e)/(1-e)})$ . La energía específica inicial es $-1/2$ y el momento angular específico es $\sqrt{1-e^2}$ . Ambos son invariantes que pueden comprobarse en cualquier punto de la órbita exacta. Una trayectoria con semieje mayor distinto del inicial, precesión espuria o error angular creciente necesitará explicaciones diferentes según el observable.

En el laboratorio numérico se evaluará qué invariantes conserva cada método y cuál es el error de fase que acumula. Los dos resultados deben documentarse por separado, porque una curva visualmente convincente puede esconder deriva energética, mientras que una trayectoria con fase desplazada puede respetar muy bien la geometría del sistema. La distinción volverá a aparecer en los ejemplos orbitales con regularidad e inestabilidad propias de las ecuaciones.

119.11. Ejercicios resueltos sobre estructura, estabilidad y significado del error

Ejercicio 1. Comprobar que una deriva es simpléctica. Para un grado de libertad y masa unitaria, $D_h(q,p)=(q+hp,p)$ . Su jacobiano es $A=\left(\begin{smallmatrix}1&h\\0&1\end{smallmatrix}\right)$ . Al calcular $A^TJA$ , con $J=\left(\begin{smallmatrix}0&1\\-1&0\end{smallmatrix}\right)$ , la multiplicación produce exactamente $J$ . La transformación es el flujo exacto del hamiltoniano cinético $T(p)=p^2/2$ , aunque no conserva la energía total $p^2/2+V(q)$ de un sistema general, porque desplaza $q$ sin actualizar $p$ conforme a la fuerza. El ejemplo distingue entre preservar una forma geométrica y preservar una función particular.

Ejercicio 2. Comparar los errores de Euler y Verlet. Para $\omega=1$ y $h=0,1$ , Euler explícito multiplica la energía por $1,01$ en cada paso. Después de cien iteraciones, el factor es $1,01^{100}\approx2,705$ . El defecto está demostrado analíticamente, sin necesidad de observar la simulación. En Verlet, en cambio, $E_h=\frac12[p^2+0,9975q^2]$ se mantiene exactamente en aritmética ideal. Para condiciones iniciales $q_0=1$ y $p_0=0$ , tenemos $E_h=0,49875$ . Como $q^2\le2E_h/0,9975=1$ , el error absoluto de energía física es $H-E_h=0,00125q^2\le0,00125$ . Esta diferencia pequeña y oscilatoria tiene un significado muy distinto del crecimiento exponencial anterior.

Ejercicio 3. Localizar una inestabilidad numérica sin inestabilidad física. En el oscilador con $\omega=2$ , tomemos $h=1,1$ . Entonces $h\omega=2,2>2$ . La matriz de Verlet todavía tiene determinante uno, pero su traza es $2[1-(2,2)^2/2]=-2,84$ . Sus autovalores son reales y recíprocos, con uno de módulo mayor que uno. La órbita numérica genérica crece; la exacta continúa siendo una elipse. La prueba matemática excluye atribuir el crecimiento a caos o a una bifurcación del modelo físico.

Ejercicio 4. ¿Conservar energía implica una buena efeméride? Considérese el oscilador con $h=0,2$ y $\omega=1$ . Su frecuencia numérica es $\omega_{\rm num}=10\arcsin(0,1)\approx1,0016742$ . En un tiempo $t=1000$ , el desfase acumulado respecto al oscilador exacto se aproxima a $1,6742$ radianes, aun cuando la energía modificada sea constante. Una gráfica de energía favorable no certifica que las posiciones sean utilizables como predicciones puntuales a ese horizonte. La conclusión sería diferente si el objetivo fuera reproducir una distribución de amplitudes durante un período prolongado.

Ejercicio 5. Evaluar una afirmación sobre pasos variables. Alguien sostiene que «cada paso es simpléctico, luego adaptar $h$ a la rapidez local mantiene la simplecticidad». La deducción es inválida, porque el jacobiano del mapa compuesto con $h(z)$ contiene términos proporcionales a $\nabla h$ y la identidad demostrada para $h$ constante deja de aplicarse. Para confirmar que una adaptación particular preserva la geometría necesitamos una nueva demostración, posiblemente en un espacio de fases extendido.

Estos ejercicios separan tres tareas que una investigación completa necesita, que son verificar una propiedad mediante una identidad exacta, medir otra mediante errores numéricos y decidir si el resultado es pertinente para la pregunta científica. El criterio de elección depende del observable, del tiempo de interés y de los mecanismos que deseamos estudiar.

119.12. Criterios para comparar métodos

La validación preliminar consta de cuatro comprobaciones. Primero, fijar una solución de referencia o un caso exactamente soluble donde puedan cuantificarse errores de orden y fase. Segundo, elegir invariantes que realmente pertenezcan al modelo (energía si es autónomo, momento angular si tiene simetría rotacional, otras integrales solo si están demostradas). Tercero, repetir el cálculo con pasos menores y registrar, además de la trayectoria, los observables cuya interpretación científica importa. Cuarto, examinar si el comportamiento persiste durante horizontes sucesivamente mayores o si aparece al cruzar una resonancia numérica.

El criterio de éxito cambia con el problema. Un seguimiento de satélite durante dos horas puede exigir máxima precisión posicional inmediata. La exploración de transporte lento en un hamiltoniano casi integrable puede requerir sobre todo preservación geométrica, análisis retrospectivo y estabilidad estadística sobre tiempos enormes. Un encuentro cercano puede requerir regularización o paso adaptativo. Incluso en un mismo proyecto puede ser legítimo combinar métodos, siempre que cada transición esté explicitada y auditada. El costo computacional también cuenta, porque duplicar la precisión no siempre da información nueva si la incertidumbre del modelo físico domina desde el comienzo.

La disciplina metodológica exige distinguir los defectos del algoritmo de las propiedades del sistema. Una trayectoria numérica irregular puede obedecer a caos real, pero también a una mala resolución de fuerzas intensas o a inestabilidad del esquema. Una trayectoria visualmente estable puede esconder grandes errores de fase. La preservación simpléctica da una garantía precisa, con límites precisos. No sustituye el análisis de convergencia, ni la comparación de modelos ni el estudio de la incertidumbre de las condiciones iniciales.

Una comprobación estructural adicional puede hacerse con el momento angular en el problema central. Si la fuerza es radial, el impulso de Verlet cambia el momento lineal por un vector paralelo a la posición, $p\mapsto p+\alpha q$ . El momento angular $L=q\times p$ no cambia en esa subetapa, porque $q\times(\alpha q)=0$ . Durante la deriva, $q\mapsto q+\beta p$ y, nuevamente, $(q+\beta p)\times p=q\times p$ . En consecuencia, la composición preserva exactamente el momento angular en aritmética ideal, mientras aproxima la energía mediante una oscilación pequeña. El resultado no depende del orden del método. Se debe a que cada subflujo respeta una simetría del problema. Una fuerza tangencial, un potencial anisótropo o una modificación arbitraria de los impulsos alterarían la demostración.

Una comparación convincente incluye reducción de paso, contraste con soluciones conocidas cuando existan, repetición con métodos alternativos y registro de invariantes y observables físicos relevantes. Los resultados que persisten bajo esas comprobaciones son más sólidos que los que dependen de una sola gráfica de trayectoria.

Referencias

Hairer, E., Lubich, C., y Wanner, G. (2006). Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations (2.ª ed.). Springer. https://doi.org/10.1007/3-540-30666-8
Hairer, E. (2005). Important aspects of geometric numerical integration. Journal of Scientific Computing, 25, 67–81. https://doi.org/10.1007/BF02728983
Leimkuhler, B., y Reich, S. (2004). Simulating Hamiltonian Dynamics. Cambridge University Press.
Sanz-Serna, J. M., y Calvo, M. P. (1994). Numerical Hamiltonian Problems. Chapman & Hall.
Wisdom, J., y Holman, M. (1991). Symplectic maps for the N-body problem. The Astronomical Journal, 102, 1528–1538. https://doi.org/10.1086/115978
Rein, H., y Tamayo, D. (2015). WHFast: A fast and unbiased implementation of a symplectic Wisdom–Holman integrator for long-term gravitational simulations. Monthly Notices of the Royal Astronomical Society, 452(1), 376–388. https://doi.org/10.1093/mnras/stv1257
Rein, H., y Spiegel, D. S. (2015). IAS15: A fast, adaptive, high-order integrator for gravitational dynamics, accurate to machine precision over a billion orbits. Monthly Notices of the Royal Astronomical Society, 446(2), 1424–1437. https://doi.org/10.1093/mnras/stu2164