Lektion 8: Kräfte und Integration

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

Kein Python auf diesem Computer? Das Notebook führt dasselbe Skript, die Übungen und die Lösungen in deinem Browser aus, in Google Colab. Das Notebook ist auf Englisch.

Zeile für Zeile

1. def pull(p):

   Die Beschleunigung jedes Körpers, durch die Gravitation der beiden anderen.

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

   Die Gravitation G·m / r² entlang der Richtung d / r ergibt zusammen d / r³; hier sind G und die Massen 1.

3. vel += acc * dt / 2

   Leapfrog: ein halber Schritt der Geschwindigkeit, ein ganzer Schritt der Position, dann der andere halbe Schritt mit den neuen Kräften.

4. pos += vel * dt

   Die Positionen bewegen sich mit den zur Hälfte aktualisierten Geschwindigkeiten.

5. track = np.array(track)

   Alle Schritte als ein Array: track[Schritt, Körper, x oder y].

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

   Die drei Körper dort, wo sie am Ende sind, auf einer gemeinsamen Bahn, eine Drittelperiode auseinander.

Probier es selbst

1. Wo ist der Schwerpunkt am Ende?

   Hinweis

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

   Lösung

   Bei (0, 0), auf etwa 10⁻¹³ genau: Der Gesamtimpuls ist null, also bewegt er sich nie.

2. Setze turn = 15, wie im Film. Was passiert, und wohin wandert der Schwerpunkt jetzt?

   Hinweis

   Vergleiche den Start mit den Positionen nach einer Periode.

   Lösung

   Die Bahnen schließen sich nicht mehr: Nach einer Periode sind die drei Körper zusammen etwa 2 Einheiten von ihrem Ausgangspunkt entfernt. Das Drehen einer Geschwindigkeit gibt dem System außerdem etwas Impuls, also driftet der Schwerpunkt, am Ende bis etwa (−0.81, 0.67).

3. Mach dt zehnmal größer, 0.005. Schließt sich die Acht noch?

   Hinweis

   Miss, wie weit die Körper nach einer Periode vom Start entfernt sind.

   Lösung

   Ja, aber weniger genau: Nach einer Periode sind die drei zusammen etwa 0,006 vom Start entfernt, gegenüber etwa 0,0001 mit dt = 0.0005.

Quiz