1. El problema de siempre: ¿cómo avanzar en el tiempo?
Imagina que tienes un planeta girando en torno a una estrella. Tú conoces sus posiciones y velocidades en el presente, y la fuerza que actúa (gravedad). La pregunta eterna: ¿dónde estará un pasito después?
Aquí entra la integración numérica. En la compu, el tiempo no es continuo: vivimos de ticks discretos \(\Delta t\). Pero cuidado: si eliges mal el método, tu planeta puede terminar fugándose a otra galaxia solo porque el algoritmo le metió más energía de la cuenta. 😅
2. El primo bruto: Euler
El método de Euler dice:
\[
\vec x_{t+\Delta t} \approx \vec x_t + \vec v_t \,\Delta t,
\qquad
\vec v_{t+\Delta t} \approx \vec v_t + \vec a_t\,\Delta t.
\]
Es directo y fácil, pero no conserva la energía. Resultado: órbitas que deberían ser elipses terminan como espirales abiertas o colapsos al centro. Un desastre.
Analogía musical: es como tocar batería sin metrónomo. Empiezas bien, pero al tercer compás ya vas a destiempo. 🥁
3. La idea de Verlet (1967)
Loup Verlet, trabajando en dinámica molecular, pensó: “¿Y si uso tanto la aceleración presente como la futura para equilibrar la velocidad?” Boom:
\[
\vec x_{t+\Delta t} = \vec x_t + \vec v_t \,\Delta t
+ \tfrac{1}{2}\vec a_t\,\Delta t^2,
\]
\[
\vec v_{t+\Delta t} = \vec v_t
+ \tfrac{1}{2}\left(\vec a_t + \vec a_{t+\Delta t}\right)\Delta t.
\]
Esto es el Velocity-Verlet, el que usamos en el bloque N-cuerpos. Lo mágico: es un integrador simpléctico. Palabra rimbombante que significa que respeta la geometría del espacio de fases: la energía no se “derrite” tan fácil.
4. Comparativa visual 👀
Vamos a ver un planeta ficticio (masa \(m=1\), órbita circular alrededor de \(M=10\)). Tres métodos con el mismo \(\Delta t\):
- Euler: órbita se deforma en pocas vueltas.
- Runge-Kutta 4: bastante bueno, pero caro en cómputo.
- Velocity-Verlet: estable, casi conserva energía durante cientos de órbitas.
(Aquí entra tu gráfico comparativo de deriva de energía. 🎨)
5. Derivación exprés (sin tanto dolor) 📜
De Taylor expandido en serie:
\[
\vec x(t+\Delta t) = \vec x(t) + \vec v(t)\Delta t + \tfrac12\vec a(t)\Delta t^2
+ \mathcal O(\Delta t^3).
\]
Lo mismo hacia atrás:
\[
\vec x(t-\Delta t) = \vec x(t) - \vec v(t)\Delta t + \tfrac12\vec a(t)\Delta t^2
+ \mathcal O(\Delta t^3).
\]
Sumando ambos y eliminando la velocidad, obtienes:
\[
\vec x(t+\Delta t) = 2\vec x(t) - \vec x(t-\Delta t) + \vec a(t)\Delta t^2.
\]
Ese es el Verlet clásico (posición–posición). La versión velocity viene de añadir la velocidad explícita usando \(\vec a_{t+\Delta t}\). Diferentes sabores, mismo espíritu.
6. Laboratorio rápido en Python 🐍
import numpy as np
import matplotlib.pyplot as plt
G, M = 1.0, 10.0
DT = 0.05 # súbele para ver el desastre de Euler
STEPS = 4000
EPS = 0.0 # sin softening para que el error pegue más fuerte
# estado inicial
pos0 = np.array([1.0, 0.0], dtype=float)
v_circ = np.sqrt(G*M/np.linalg.norm(pos0))
v_factor = 1.10 # >1: espiral hacia afuera; <1: hacia adentro
vel0 = np.array([0.0, v_circ*v_factor], dtype=float)
def acc_from(p):
r2 = np.dot(p, p) + EPS**2
return -G*M*p / (r2**1.5)
def energy(p, v):
r = np.linalg.norm(p)
KE = 0.5*np.dot(v, v)
PE = -G*M/r
return KE + PE
# --- Euler ---
p_e = pos0.copy()
v_e = vel0.copy()
traj_e = []
E_e = []
# --- Velocity-Verlet ---
p_v = pos0.copy()
v_v = vel0.copy()
a_v = acc_from(p_v)
traj_v = []
E_v = []
for _ in range(STEPS):
# Euler
a_e = acc_from(p_e)
p_e = p_e + v_e*DT
v_e = v_e + a_e*DT
traj_e.append(p_e.copy())
E_e.append(energy(p_e, v_e))
# Velocity-Verlet
p_v = p_v + v_v*DT + 0.5*a_v*(DT**2)
new_a_v = acc_from(p_v)
v_v = v_v + 0.5*(a_v + new_a_v)*DT
a_v = new_a_v
traj_v.append(p_v.copy())
E_v.append(energy(p_v, v_v))
traj_e = np.array(traj_e); traj_v = np.array(traj_v)
t = np.arange(STEPS)*DT
E0_e, E0_v = E_e[0], E_v[0]
rel_e = (np.array(E_e)-E0_e)/abs(E0_e)
rel_v = (np.array(E_v)-E0_v)/abs(E0_v)
plt.figure(figsize=(11,5))
# trayectorias
plt.subplot(1,2,1)
plt.plot(traj_e[:,0], traj_e[:,1], label='Euler')
plt.plot(traj_v[:,0], traj_v[:,1], label='Velocity-Verlet')
plt.gca().set_aspect('equal', 'box')
plt.xlim(-8, 2.5); plt.ylim(-8, 2.5)
plt.xlabel('x'); plt.ylabel('y'); plt.title('Trayectorias')
plt.legend()
# deriva de energía
plt.subplot(1,2,2)
plt.plot(t, rel_e, label='Euler')
plt.plot(t, rel_v, label='Verlet')
plt.xlabel('t'); plt.ylabel('$(E-E_0)/|E_0|$'); plt.title('Deriva de energía')
plt.legend()
plt.tight_layout()
plt.show()
Este código traza dos trayectorias y la deriva de energía; con Euler verás la espiral, con Verlet la órbita estable.
Figura 1: Comparativa de Euler vs Velocity-Verlet: la órbita de Euler se espirala por no conservar la geometría del espacio de fases; Verlet mantiene energía casi constante.
Epílogo filosófico ✨
Lo bonito es que Verlet no es solo un truco: es un recordatorio de que respetar las simetrías y estructuras de un sistema físico es más importante que sumar con muchas cifras decimales.