Leçon 8 : Forces et intégration

"""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

Pas de Python sur cet ordinateur ? Le notebook exécute le même script, les exercices et les réponses dans ton navigateur, sur Google Colab. Le notebook est en anglais.

Ligne par ligne

1. def pull(p):

   L’accélération de chaque corps, due à la gravité des deux autres.

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

   La gravité G·m / r², le long de la direction d / r, donne en tout d / r³ ; ici, G et les masses valent 1.

3. vel += acc * dt / 2

   Saute-mouton : un demi-pas de vitesse, un pas entier de position, puis l’autre demi-pas avec les nouvelles forces.

4. pos += vel * dt

   Les positions avancent avec les vitesses mises à jour à moitié.

5. track = np.array(track)

   Tous les pas dans un seul tableau : track[pas, corps, x ou y].

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

   Les trois corps là où ils ont fini, sur une même trajectoire partagée, à un tiers de période les uns des autres.

À toi de jouer

1. Où se trouve le centre de masse à la fin ?

   Indice

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

   Réponse

   En (0, 0), à environ 10⁻¹³ près : la quantité de mouvement totale est nulle, donc il ne bouge jamais.

2. Mets turn = 15, comme dans le film. Que se passe-t-il, et où va maintenant le centre de masse ?

   Indice

   Compare le départ avec les positions après une période.

   Réponse

   Les trajectoires ne se referment plus : après une période, les trois corps sont, ensemble, à environ 2 unités de leur point de départ. Tourner une vitesse donne aussi au système une quantité de mouvement, donc le centre de masse dérive, jusqu’à environ (−0.81, 0.67) à la fin.

3. Rends dt dix fois plus grand, 0.005. Le huit se referme-t-il toujours ?

   Indice

   Mesure à quelle distance les corps sont du départ après une période.

   Réponse

   Oui, mais moins exactement : après une période, les trois sont, ensemble, à environ 0,006 du départ, contre environ 0,0001 avec dt = 0.0005.

Quiz