Задача N тел: гравитация каждого с каждым

Добавь третье тело — и точного решения больше нет; задача N тел живёт только в симуляции.

Задача N тел — расчёт движения системы из $N$ гравитирующих тел, каждое из которых притягивает каждое; для $N \ge 3$ общего аналитического решения не существует.

Почему N тел — это сложно

Две гравитирующие точки решаются аналитически — это орбиты Кеплера. Но стоит добавить третье тело, и замкнутого решения уже нет: система становится хаотической, чувствительной к начальным условиям. Зато численно задача N тел решается прямолинейно. Идея проста: ускорение каждого тела — это сумма притяжений со стороны всех остальных. Если тел $N$, то на каждом шаге для каждого тела складываем $N-1$ вкладов:

$$\vec a_i = \sum_{j \neq i} G\frac{m_j}{r_{ij}^2}\,\hat r_{ij},$$

где $r_{ij}$ — расстояние между телами $i$ и $j$, а $\hat r_{ij}$ — единичный вектор от $i$ к $j$. Это двойной цикл «каждый с каждым», стоимостью $O(N^2)$ операций на шаг.

Симуляция системы из трёх тел

Соберём звезду и две планеты. Звезду не закрепляем — пусть тоже движется (она чуть качается от притяжения планет, это и есть честная N-body постановка). Главный тест корректности — суммарный импульс системы: на него не действуют внешние силы, поэтому он обязан сохраняться.

import math
G = 4*math.pi**2
bodies = [
    {"m": 1.0,  "x": 0.0,  "y": 0.0, "vx": 0.0, "vy": 0.0},                # звезда
    {"m": 1e-3, "x": 1.0,  "y": 0.0, "vx": 0.0, "vy": math.sqrt(G/1.0)},   # планета A
    {"m": 1e-3, "x": -2.0, "y": 0.0, "vx": 0.0, "vy": -math.sqrt(G/2.0)},  # планета B
]
def accelerations(bodies):
    n = len(bodies)
    acc = [[0.0, 0.0] for _ in range(n)]
    for i in range(n):
        for j in range(n):
            if i == j:
                continue
            dx = bodies[j]["x"] - bodies[i]["x"]
            dy = bodies[j]["y"] - bodies[i]["y"]
            r = math.hypot(dx, dy)
            f = G*bodies[j]["m"]/r**3      # сила/масса_i на единицу dx,dy
            acc[i][0] += f*dx; acc[i][1] += f*dy
    return acc
def total_p(bodies):
    return (sum(b["m"]*b["vx"] for b in bodies),
            sum(b["m"]*b["vy"] for b in bodies))

dt = 0.001
for s in range(0, 2001):
    if s % 400 == 0:
        rA = math.hypot(bodies[1]["x"], bodies[1]["y"])
        rB = math.hypot(bodies[2]["x"], bodies[2]["y"])
        px, py = total_p(bodies)
        print(f"t={s*dt:.2f}  r_A={rA:.3f}  r_B={rB:.3f}  P=({px:+.6f}, {py:+.6f})")
    acc = accelerations(bodies)
    for i, b in enumerate(bodies):
        b["vx"] += acc[i][0]*dt; b["vy"] += acc[i][1]*dt
    for b in bodies:
        b["x"] += b["vx"]*dt; b["y"] += b["vy"]*dt

Вывод:

t=0.00  r_A=1.000  r_B=2.000  P=(+0.000000, +0.001840)
t=0.40  r_A=0.996  r_B=1.998  P=(+0.000000, +0.001840)
t=0.80  r_A=0.997  r_B=1.998  P=(-0.000000, +0.001840)
t=1.20  r_A=0.998  r_B=2.000  P=(+0.000000, +0.001840)
t=1.60  r_A=0.998  r_B=2.004  P=(+0.000000, +0.001840)
t=2.00  r_A=0.999  r_B=2.007  P=(+0.000000, +0.001840)

Обе планеты держат свои орбиты (радиусы около $1$ и $2$ а.е.), а суммарный импульс системы стоит как вкопанный — $P_y = 0.001840$ от начала до конца. Это бесценная проверка: если бы импульс «поплыл», в коде была бы ошибка (например, несимметричный учёт сил). Сохранение суммарного импульса — встроенный детектор багов любого N-body движка.

Как работает под капотом

Двойной цикл $O(N^2)$ хорош для десятков тел, но для миллионов звёзд в галактике он непосилен. Тогда применяют умные приближения: метод Барнса — Хатта группирует далёкие тела в «облака» и считает их как одну точку, снижая стоимость до $O(N\log N)$; particle-mesh методы раскладывают гравитацию через сетку. Но идея неизменна: ускорение — сумма вкладов. Ещё одна тонкость — смягчение (softening): при близком сближении $r \to 0$ сила $\frac{1}{r^2}$ взрывается и симуляция разлетается; поэтому в знаменатель добавляют малую константу $\varepsilon$, заменяя $r^2$ на $r^2 + \varepsilon^2$. Это убирает численные «выстрелы» при тесных контактах ценой крошечного искажения близкодействия.

Частые ошибки

Закреплять звезду без нужды. Для честной N-body все тела движутся; фиксация центра нарушает сохранение импульса.
Делить на $r^2$ без смягчения. При тесном сближении сила взрывается и тела «выстреливают» в бесконечность; нужен softening $\varepsilon$.
Не проверять сохранение импульса/энергии. Это главные индикаторы корректности; без них баг в суммировании сил незаметен.

Итог

Ускорение тела — сумма притяжений всех остальных ($O(N^2)$ двойной цикл).
Для $N \ge 3$ аналитики нет — задача живёт только в симуляции.
Сохранение суммарного импульса — встроенный тест корректности движка.
Для масштаба и устойчивости нужны приближения ($O(N\log N)$) и смягчение силы $\varepsilon$.