Lección 8: Fuerzas e integración

"""Three bodies: the figure-eight orbit, from the starting values computed by Carles Simo
(in Chenciner and Montgomery, 2000). G = 1, three masses of 1.
Run: python three-body-orbits.py   (needs: pip install numpy matplotlib)
"""
import numpy as np
import matplotlib.pyplot as plt

plt.style.use("dark_background")

x1 = np.array([0.97000436, -0.24308753])
V = np.array([-0.93240737, -0.86473146])  # the velocity of body 3
pos = np.array([x1, -x1, [0.0, 0.0]])
vel = np.array([-V / 2, -V / 2, V])
period, dt = 6.32591398, 0.0005
turn = 0  # degrees to turn the first velocity by (the film breaks the eight with 15)

a = np.radians(turn)
vel[0] = [vel[0, 0] * np.cos(a) - vel[0, 1] * np.sin(a), vel[0, 0] * np.sin(a) + vel[0, 1] * np.cos(a)]


def pull(p):
    """Each body's acceleration: the sum of the other two's gravity, G*m / r**2 towards each"""
    acc = np.zeros_like(p)
    for i in range(3):
        for j in range(3):
            if i != j:
                d = p[j] - p[i]
                acc[i] += d / np.linalg.norm(d) ** 3
    return acc


track = []
acc = pull(pos)
for _ in range(int(3 * period / dt)):  # three periods, leapfrog steps (they keep the energy well)
    vel += acc * dt / 2
    pos += vel * dt
    acc = pull(pos)
    vel += acc * dt / 2
    track.append(pos.copy())
track = np.array(track)

colours = ["#5eead4", "#f0abfc", "#f5c76b"]
plt.figure(figsize=(7, 4))
for body in range(3):  # all three run along the same track, a third of a period apart
    plt.plot(track[:, body, 0], track[:, body, 1], color=colours[body], linewidth=1)
plt.scatter(track[-1, :, 0], track[-1, :, 1], c=colours, s=80, zorder=3)
plt.axis("equal")
plt.axis("off")
plt.show()

# Try: turn = 15 -> the three tracks no longer close; the eight falls apart

¿No tienes Python en este equipo? El cuaderno ejecuta el mismo script, los ejercicios y las respuestas en tu navegador, en Google Colab. El cuaderno está en inglés.

Línea a línea

1. def pull(p):

   La aceleración de cada cuerpo, debida a la gravedad de los otros dos.

2. acc[i] += d / np.linalg.norm(d) ** 3

   La gravedad G·m / r² a lo largo de la dirección d / r da, en conjunto, d / r³; aquí G y las masas valen 1.

3. vel += acc * dt / 2

   Salto de rana: medio paso de velocidad, un paso entero de posición y luego el otro medio paso con las fuerzas nuevas.

4. pos += vel * dt

   Las posiciones avanzan con las velocidades actualizadas a medias.

5. track = np.array(track)

   Todos los pasos como un solo array: track[paso, cuerpo, x o y].

6. plt.scatter(track[-1, :, 0], track[-1, :, 1], c=colours, s=80, zorder=3)

   Los tres cuerpos donde terminaron, en una misma pista compartida, separados por un tercio de periodo.

Pruébalo tú

1. ¿Dónde está el centro de masas al final?

   Pista

   track[-1].mean(axis=0).

   Respuesta

   En (0, 0), con un margen de unos 10⁻¹³: el momento total es cero, así que nunca se mueve.

2. Pon turn = 15, como en el vídeo. ¿Qué pasa, y adónde va ahora el centro de masas?

   Pista

   Compara el inicio con las posiciones tras un periodo.

   Respuesta

   Las trayectorias ya no se cierran: tras un periodo, los tres cuerpos están, en conjunto, a unas 2 unidades de donde empezaron. Girar una velocidad también le da al sistema algo de momento, así que el centro de masas se desplaza, hasta cerca de (−0.81, 0.67) al final.

3. Haz dt diez veces mayor, 0.005. ¿Sigue cerrándose el ocho?

   Pista

   Mide a qué distancia del inicio están los cuerpos tras un periodo.

   Respuesta

   Sí, pero con menos exactitud: tras un periodo, los tres están, en conjunto, a unos 0,006 del inicio, frente a unos 0,0001 con dt = 0.0005.

Cuestionario