Задача N тел: прямое суммирование

Несколько тел под взаимным притяжением: прямое суммирование и законы сохранения.

Задача N тел — движение N масс, каждая из которых притягивает каждую; при N≥3 точного решения нет, и единственный путь — численное интегрирование.

Почему три тела — это уже хаос

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

Прямое суммирование O(N²)

Самый честный метод: для каждого тела складываем притяжение от всех остальных. Если тел N, то на каждом шаге нужно N·(N-1) вычислений силы — сложность O(N²). Для трёх тел это пустяк, для миллиарда — катастрофа (об этом позже). Прогоним три тела, отслеживая энергию и импульс:

import math
G = 1.0
bodies = [
    {"m":2.0, "x":0.0,  "y":0.0,  "vx":0.0,  "vy":0.0},
    {"m":1.0, "x":1.0,  "y":0.0,  "vx":0.0,  "vy":1.3},
    {"m":1.0, "x":-1.0, "y":0.0,  "vx":0.0,  "vy":-1.3},
]
def accelerations(bs):
    n=len(bs); ax=[0.0]*n; ay=[0.0]*n
    for i in range(n):
        for j in range(n):
            if i==j: continue
            dx=bs[j]["x"]-bs[i]["x"]; dy=bs[j]["y"]-bs[i]["y"]
            r=math.hypot(dx,dy)+1e-9
            f=G*bs[j]["m"]/r**3
            ax[i]+=f*dx; ay[i]+=f*dy
    return ax, ay
def energy(bs):
    E=sum(0.5*b["m"]*(b["vx"]**2+b["vy"]**2) for b in bs)
    n=len(bs)
    for i in range(n):
        for j in range(i+1,n):
            dx=bs[j]["x"]-bs[i]["x"]; dy=bs[j]["y"]-bs[i]["y"]
            r=math.hypot(dx,dy)+1e-9
            E-=G*bs[i]["m"]*bs[j]["m"]/r
    return E
def momentum(bs):
    return sum(b["m"]*b["vx"] for b in bs), sum(b["m"]*b["vy"] for b in bs)
dt=0.0005
ax,ay=accelerations(bodies)
E0=energy(bodies); p0=momentum(bodies)
for step in range(20000):
    for i,b in enumerate(bodies):
        b["x"]+=b["vx"]*dt+0.5*ax[i]*dt*dt
        b["y"]+=b["vy"]*dt+0.5*ay[i]*dt*dt
    ax2,ay2=accelerations(bodies)
    for i,b in enumerate(bodies):
        b["vx"]+=0.5*(ax[i]+ax2[i])*dt; b["vy"]+=0.5*(ay[i]+ay2[i])*dt
    ax,ay=ax2,ay2
E1=energy(bodies); p1=momentum(bodies)
print(f"Шагов: 20000, прошло t={20000*dt:.1f}")
print(f"Энергия:  старт {E0:+.4f}  конец {E1:+.4f}")
print(f"Импульс X: {p0[0]:+.4f} -> {p1[0]:+.4f}")
print(f"Импульс Y: {p0[1]:+.4f} -> {p1[1]:+.4f}")

Вывод:

Шагов: 20000, прошло t=10.0
Энергия:  старт -2.8100  конец -2.8100
Импульс X: +0.0000 -> +0.0000
Импульс Y: +0.0000 -> +0.0000

За 20 000 шагов и довольно сложную динамику энергия не сдвинулась (−2.8100 → −2.8100), а импульс остался строго нулевым. Импульс сохраняется точно (третий закон Ньютона: силы между телами равны и противоположны, в сумме ноль), энергия — почти точно благодаря симплектическому Верле. Два независимых закона сохранения подтверждают, что симуляция верна.

Устойчивость конфигураций

Не всякая система N тел устойчива. Большинство случайных начальных условий приводят к тому, что одно тело набирает скорость и выбрасывается прочь, а остальные образуют тесную пару. Но есть и удивительно устойчивые решения: точки Лагранжа, орбиты-«восьмёрки» трёх равных масс. Симуляция — единственный способ их найти и проверить на устойчивость.

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

Сохранение импульса в коде — прямое следствие симметрии: вклад тела j в ускорение i равен по величине и противоположен вкладу i в ускорение j. При суммировании по всем парам эти вклады в полный импульс взаимно уничтожаются. Поэтому импульс сохраняется машинно точно, до последнего бита, независимо от шага dt. Энергия же сохраняется лишь приближённо (зависит от метода и шага) — вот почему её дрейф используют как индикатор качества, а импульс — как индикатор корректности кода.

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

Считать силу дважды или забыть пару. Двойной цикл по i и j должен покрыть все упорядоченные пары, кроме i=j.
Не смягчать особенность при сближении. При r→0 сила взрывается; малое слагаемое в знаменателе (softening) спасает.
Обновлять позиции и скорости вперемешку с расчётом сил. Сначала считают все ускорения, потом обновляют — иначе нарушается синхронность.

Итоги

Задача N тел при N≥3 не имеет точного решения — только численно.
Прямое суммирование имеет сложность O(N²).
Импульс сохраняется точно (третий закон Ньютона), энергия — приближённо.
Большинство конфигураций неустойчивы; устойчивые ищут симуляцией.