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.