N-body гравітація: симулюємо космос у браузері

Задача N тіл — одна з найстаріших у фізиці: маючи N точкових мас, що взаємодіють через гравітацію, передбачити їхні майбутні положення. Для N = 2 Ньютон розв'язав її точно ще у сімнадцятому столітті — задача двох тіл має замкнутий аналітичний розв'язок у вигляді еліпсів, парабол і гіпербол. Для N = 3 Пуанкаре довів у 1889 році, що загального замкнутого розв'язку не існує. Для N = 10 000 — масштабу нашої симуляції галактики — потрібен ефективний алгоритм і ретельні числові прийоми.

Метод грубої сили: O(N²)

Найпростіший спосіб симулювати N гравітуючих тіл — обчислити для кожного тіла гравітаційну силу від кожного іншого тіла. Закон всесвітнього тяжіння Ньютона дає силу між тілами i та j як:

F_ij = G m_i m_j / |r_i − r_j|² у напрямку від i до j.

Для кожного з N тіл потрібно підсумувати N − 1 парних сил. Це O(N²) обчислень сили на кожен часовий крок. Для N = 1 000 це один мільйон обчислень на кадр. Для N = 10 000 — сто мільйонів. При 60 кадрах за секунду симуляція методом грубої сили для 10 000 тіл вимагає шести мільярдів обчислень сили за секунду — надто повільно для реального часу в браузері.

Проте метод грубої сили корисний при невеликому N — менше приблизно 500 тіл — де він максимально точний і простий у реалізації. Наша симуляція використовує грубу силу для невеликої кількості частинок і перемикається на Barnes-Hut, коли їх кількість перевищує певний поріг.

Barnes-Hut: O(N log N) за допомогою квадродерева

Алгоритм Barnes-Hut, опублікований Джошем Барнсом і Пітом Хатом у 1986 році, знижує складність до O(N log N), використовуючи ключове наближення: віддалені групи тіл можна розглядати як єдине тіло в їхньому центрі мас.

Алгоритм будує квадродерево (у 2D) або октодерево (у 3D), яке рекурсивно розбиває простір на комірки. Кожна комірка зберігає сумарну масу і центр мас усіх тіл усередині неї. Щоб обчислити силу, що діє на тіло i, потрібно пройти по дереву. Для кожного вузла: якщо вузол достатньо далеко відносно свого розміру — а саме якщо ширина вузла s, поділена на відстань до нього d, задовольняє s/d < θ (параметр налаштування, зазвичай 0,5–1,0) — весь вузол розглядається як єдине тіло в його центрі мас. Інакше відбувається рекурсивний спуск до дочірніх вузлів.

Алгоритм	Складність	N = 1 000	N = 10 000
Груба сила	O(N²)	1М операцій	100М операцій
Barnes-Hut	O(N log N)	~10 тис. операцій	~133 тис. операцій
Fast Multipole	O(N)	~1 тис. операцій	~10 тис. операцій

Параметр θ керує компромісом між точністю та швидкістю. θ = 0 дає точні результати грубої сили. θ = 1,0 дає значні похибки наближення. θ = 0,5 забезпечує хороший баланс — похибка сили менше 1% для більшості тіл при збереженні інтерактивності симуляції.

Гравітаційне згладжування: запобігання сингулярностям

Закон гравітації Ньютона має сингулярність: коли два тіла зближуються, сила між ними зростає необмежено (1/r² → нескінченність при r → 0). У реальній галактиці справжні близькі зустрічі між зірками надзвичайно рідкісні, оскільки зірки мізерно малі порівняно з відстанями між ними. У симуляції ж частинки можуть зближуватися як завгодно, спричиняючи числові вибухи — один часовий крок дає величезну силу, породжує величезну швидкість, і частинка вилітає в нескінченність.

Стандартне вирішення — гравітаційне згладжування: замінити знаменник r² на r² + ε², де ε — довжина згладжування. Це обмежує максимальну силу скінченним значенням, коли дві частинки перебувають на відстані ε одна від одної. Тепер сила виглядає так:

F = G m_i m_j / (r² + ε²)

Фізично це еквівалентно розмазуванню кожної точкової маси в сферу радіусом ε. Для значень r, набагато більших за ε, сила невідрізнима від сили точкових мас. Для дуже малих r сила обмежена значенням G m_i m_j / ε².

Вибір ε — це мистецтво. Занадто велике значення пригнічує реальні динамічні ефекти — частинки, які мали б «пращею» відкидати одна одну, натомість плавно проходять крізь одна одну. Занадто мале значення — і близькі зустрічі домінують над часовим кроком, різко сповільнюючи симуляцію. Типовий вибір для симуляції галактики — ε, що дорівнює приблизно 1% середньої відстані між частинками.

Інтегрування за часом: метод leapfrog

Маючи сили, як просувати положення та швидкості частинок у часі? Простий метод Ейлера — оновити швидкість, потім положення, використовуючи поточну силу — швидкий, але з часом втрачає енергію, через що орбіти повільно розкручуються назовні. Для N-body симуляції, що працює тисячі часових кроків, цей дрейф енергії псує результати.

Інтегратор leapfrog, також званий методом Штермера-Верле, вирішує це за допомогою зсунутого оновлення. Швидкості обчислюються в напівцілих часових кроках, а положення — у цілих:

• v_n+1/2 = v_n-1/2 + a_n Δt
• x_n+1 = x_n + v_n+1/2 Δt

Ця крихітна зміна — швидкості й положення більше не оновлюються одночасно — робить інтегратор симплектичним: він точно зберігає трохи модифікований гамільтоніан (функцію енергії) системи. На практиці це означає, що повна енергія коливається навколо істинного значення, а не дрейфує — орбіти залишаються стабільними протягом мільйонів часових кроків, не стаючи нефізичними.

Для ще кращої точності у високоточних астрономічних симуляціях використовують симплектичні інтегратори вищого порядку, як-от метод четвертого порядку Йошиди чи класичний метод Рунге-Кутти четвертого порядку (несимплектичний, але дуже точний на кожному кроці). Наша браузерна симуляція використовує leapfrog за його простоту, стабільність і відмінне збереження енергії.

Формування галактик і емерджентна структура

Запустіть достатню кількість гравітуючих частинок із відповідними початковими умовами — і спонтанно виникає вражаюча структура. Надайте хмарі частинок певний початковий кутовий момент і дайте гравітації її стиснути: хмара сплющується в диск, під дією нестійкості Джинса формуються згустки, а диференціальне обертання породжує спіральні рукави. Це не запрограмовані властивості — вони виникають із N простих частинок, що підкоряються закону Ньютона.

Справжні космологічні симуляції, як-от проєкти Illustris та EAGLE, використовують ті самі фундаментальні алгоритми — Barnes-Hut або Fast Multipole для гравітації, SPH або рухому сітку для динаміки газу — масштабовані до мільярдів частинок на суперкомп'ютерних кластерах. Якісна поведінка, яку вони відтворюють, напрочуд добре узгоджується зі спостережуваною морфологією галактик. Наша браузерна симуляція — це мініатюрна версія тієї самої обчислювальної фізики, яку професійні космологи використовують для розуміння великомасштабної структури Всесвіту.