Закон всесвітнього тяжіння Ньютона
Кожна маса у всесвіті притягує кожну іншу масу із силою, що визначається законом Ньютона:
F = G × m₁ × m₂ / r²Де G — гравітаційна стала (6.674 × 10−11 Н·м²·кг−2), m₁ і m₂ — маси двох тіл, а r — відстань між їхніми центрами. Сила є силою притягання і діє вздовж лінії, що з'єднує обидва тіла.
Для двох тіл це дає систему, яку можна розв'язати повністю. Тіла обертаються навколо спільного центра мас по еліпсах — закони Кеплера випливають безпосередньо. Додайте третє тіло — і все змінюється.
Чому задача N тіл складна: задача трьох тіл
У 1887 році Анрі Пуанкаре довів, що гравітаційна задача трьох тіл не має загального розв'язку у замкненій формі. Неможливо записати рівняння, яке дає положення трьох гравітуючих мас у момент часу t на основі їхніх початкових умов. Система, як правило, хаотична — малі відмінності в початкових умовах призводять до цілком різних довгострокових траєкторій.
Це не провал математики. Це фундаментальна властивість самих рівнянь. Три гравітаційно взаємодіючих тіла можуть демонструвати будь-який режим поведінки: стабільні орбіти, квазіперіодичний рух, резонанси, хаотичне блукання і викидання одного з тіл у нескінченність.
Єдиний практичний спосіб передбачити еволюцію системи N тіл — чисельне інтегрування: робити малі кроки за часом, обчислювати сили, оновлювати швидкості й положення, повторювати. Саме це і робить симуляція — і саме це робить кожен астрофізичний код, від моделей Сонячної системи до симуляцій галактик.
Наївний алгоритм: складність O(n²)
Прямолінійний підхід до обчислення гравітаційних сил полягає в перевірці кожної пари тіл:
for each body i:
force[i] = (0, 0, 0)
for each body j ≠ i:
r = position[j] - position[i]
dist = length(r)
force[i] += G * mass[i] * mass[j] / dist² * normalize(r)Це вимагає n × (n−1) обчислень сили на кожен крок за часом — складність O(n²). Для 100 тіл: 9 900 обчислень. Для 1 000 тіл: 999 000 обчислень. Для 10 000 тіл: майже 100 мільйонів обчислень. Для справжньої галактики зі 100 мільярдами зірок наївний алгоритм повністю нездійсненний.
Пом'якшення для запобігання сингулярностям
У рівнянні сили криється чисельна проблема: коли два тіла зближуються дуже сильно, r прямує до нуля, а сила прямує до нескінченності. Насправді зірки не є точковими масами — близькі зближення включають складну фізику. У симуляції ми додаємо довжину пом'якшення ε:
F = G × m₁ × m₂ / (r² + ε²)Це обмежує максимальну силу і запобігає чисельним вибухам під час близьких зближень. Довжина пом'якшення — параметр симуляції: занадто велика — і фізика стає нереалістичною; занадто мала — і близькі зближення дестабілізують інтегратор.
Барнс-Хат: зниження до O(n log n)
У 1986 році Джош Барнс і Піт Хат опублікували алгоритм, що зменшує складність задачі N тіл з O(n²) до O(n log n) — трансформаційне покращення. Ключова ідея: віддалене скупчення багатьох тіл можна наблизити одним тілом у центрі мас скупчення. Гравітаційна похибка від такого наближення мала, коли скупчення знаходиться далеко.
Алгоритм Барнса-Хата будує октодерево (дерево, яке рекурсивно ділить тривимірний простір на вісім октантів у кожному вузлі). Кожен внутрішній вузол зберігає сумарну масу і центр мас усіх тіл у своєму піддереві.
Під час обчислення сили на тіло i дерево обходиться. У кожному вузлі приймається рішення на основі відношення s/d — де s це просторовий розмір вузла, а d — відстань від тіла i до центра мас вузла. Якщо s/d < θ (порогове значення кута розкриття, зазвичай від 0.5 до 1.0), вузол трактується як єдина точкова маса. Інакше розкриваються і обробляються окремо дочірні вузли.
function compute_force(body i, node):
if node is leaf:
add direct force from node's body
else:
s = node size
d = distance from i to node's centre of mass
if s/d < θ:
add approximate force from node's mass
else:
for each child of node:
compute_force(body i, child)При θ ≈ 0.5 похибка в оцінці сили нижче 1%, тоді як обчислення скорочуються з годин до секунд для тієї ж кількості частинок. Сучасні космологічні симуляції використовують варіанти Барнса-Хата для симуляції мільярдів частинок темної матерії.
Стрибкоподібна інтеграція: стабільність важливіша за точність
Знання сил — лише половина задачі. Потрібно ще й інтегрувати рівняння руху — оновлювати положення і швидкості — так, щоб система залишалась стабільною протягом тисяч кроків за часом.
Просте інтегрування Ейлера (x += v*dt; v += a*dt)
здається природним, але має фатальний недолік для орбітальної
механіки: воно не зберігає енергію. Повна енергія системи з часом
дрейфує вгору, а орбіти повільно розкручуються назовні. У
довгостроковій симуляції планети покидають свої зірки.
Стрибкоподібний інтегратор (leapfrog) вирішує це, зсуваючи оновлення положення і швидкості на половину кроку за часом:
// kick-drift-kick (KDK) form:
v_half = v + a * (dt/2) // half-step velocity
x = x + v_half * dt // full-step position
a = compute_forces(x) // forces at new position
v = v_half + a * (dt/2) // second half-step velocityСтрибкоподібна схема є симплектичною — вона точно зберігає модифікований гамільтоніан, близький до справжнього. На практиці це означає, що енергія коливається навколо правильного значення замість дрейфу, що дозволяє симуляціям тривати мільярди кроків за часом без накопичення похибки.
Інтегратори вищого порядку (Рунге-Кутта 4, Йошида 6-го порядку) точніші на кожному кроці, але вимагають більше обчислень сили. Для гравітаційної задачі N тіл стрибкоподібна схема з малим кроком за часом зазвичай є найкращим компромісом.
Реальні застосування в астрофізиці
Симуляція N тіл — основний інструмент сучасної астрофізики:
- Формування планетних систем — симуляція ранньої історії Сонячної системи, коли сотні планетезималей зіштовхувались і зливались упродовж мільйонів років, формуючи планети.
- Зіткнення галактик — Чумацький Шлях зіткнеться з галактикою Андромеди приблизно через 4,5 мільярда років. Симуляції N тіл з 1970-х років передбачали часову шкалу та структуру цього злиття.
- Кульові скупчення — щільні сферичні скупчення від 100 000 до 1 мільйона зірок, які симулюють для розуміння зоряної динаміки, формування чорних дір та швидкостей утечі.
- Крупномасштабна структура — космологічні симуляції на кшталт Millennium та IllustrisTNG відстежують мільярди частинок темної матерії від Великого вибуху до сьогодення, відтворюючи космічну павутину з ниток, порожнин і скупчень галактик.
- Динаміка астероїдів — траєкторії навколоземних астероїдів обчислюються методом інтегрування N тіл із урахуванням збурень від усіх восьми планет, Місяця та Плутона.
Досліджуйте гравітаційну динаміку прямо у браузері на /nbody/. Спробуйте налаштувати стабільні орбіти трьох тіл у формі вісімки, або подивіться, що станеться, якщо додати масивного гостя до планетної системи.
Стабільність нашої Сонячної системи
Чи стабільна наша власна Сонячна система? Це виявляється напрочуд складним питанням. Чисельні інтегрування, виконані на мільярди симульованих років, показують, що внутрішні планети (від Меркурія до Марса) мають малу, але ненульову ймовірність увійти в хаотичні орбіти. Орбіта Меркурія може стати нестабільною за часові масштаби у кілька мільярдів років, потенційно зіштовхнувшись з Венерою або будучи викинутою.
Зовнішні планети стабільніші, але не ідеально. Резонанс чотирьох тіл між Юпітером, Сатурном, Ураном і Нептуном означає, що малі збурення можуть накопичуватися. Сонячна система, яку ми знаємо, може бути лише однією з багатьох можливих конфігурацій, яка випадково проіснувала достатньо довго, щоб еволюціонувало складне життя.