Lição 8: Forças e integração

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

Não tem Python neste computador? O notebook executa o mesmo script, os exercícios e as respostas no navegador, no Google Colab. O notebook está em inglês.

Linha a linha

1. def pull(p):

   A aceleração de cada corpo, causada pela gravidade dos outros dois.

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

   A gravidade G·m / r² ao longo da direção d / r dá, ao todo, d / r³; aqui G e as massas valem 1.

3. vel += acc * dt / 2

   Leapfrog: meio passo de velocidade, um passo inteiro de posição e depois o outro meio passo com as forças novas.

4. pos += vel * dt

   As posições avançam com as velocidades já atualizadas em meio passo.

5. track = np.array(track)

   Todos os passos num só array: track[passo, corpo, x ou y].

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

   Os três corpos onde terminaram, numa única trajetória comum, separados por um terço de período.

Experimente

1. Onde está o centro de massa no fim?

   Pista

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

   Resposta

   Em (0, 0), com um desvio de cerca de 10⁻¹³: o momento linear total é nulo, por isso nunca se move.

2. Defina turn = 15, como no filme. O que acontece, e para onde vai agora o centro de massa?

   Pista

   Compare o início com as posições depois de um período.

   Resposta

   As trajetórias deixam de se fechar: depois de um período, os três corpos estão, somando os três, a cerca de 2 unidades do ponto onde começaram. Girar uma velocidade também dá algum momento linear ao sistema, por isso o centro de massa desloca-se, até cerca de (−0,81; 0,67) no fim.

3. Torne dt dez vezes maior, 0.005. O oito ainda se fecha?

   Pista

   Meça a que distância do início estão os corpos depois de um período.

   Resposta

   Sim, mas com menos exatidão: depois de um período, os três estão, somando os três, a cerca de 0,006 do início, contra cerca de 0,0001 com dt = 0.0005.

Quiz