Barnes-Hut: Моделювання 100 000 гравітуючих тіл у реальному часі

Стіна О(n²)

Гравітація – це парний взаємодія. Щоб просунути систему N тіл на один крок, вам потрібна прискореність кожного тіла, і пряма сума дає її точно: для кожного тіла додається тягар від кожного іншого тіла. Це n(n−1)/2 парів – O(n²) роботи на крок. Тисяча тіл означає півмільйона взаємодій за крок, і це працює комфортно зі швидкістю 60 кадрів в секунду у браузері. Сто тисяч тіл означає п'ять мільярдів взаємодій за крок, і ви тепер чекаєте хвилини на кадр.

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

Дерево

Рекурсивно діліть об’єм моделювання на частини. У 2D кожен квадрат ділиться на чотири квадранти – це квінтодра́т; у 3D кожен куб ділиться на вісім октант – це окві́нтодра́т. Ви продовжуєте ділити клітинку доти, поки вона не містить максимум одну частинку. Кожен внутрішній вузол зберігає дві агреговані величини, обчислені під час підйому: загальну масу частинок нижче нього та їх центр маси.

insert(вузол, частинка): якщо вузол є порожнім листком: вузол.частинка = частинка; повертати якщо вузол є листком з однією частинкою: розділити(вузол); повторно вставити стару частинку вузол.маса += маса частинки // агреговані величини оновлюються під час спуску вузол.центр_маси = зважена середня величина центрів мас вставити(дитина, що містить позицію частинки, частинка) Побудова дерева коштує O(n log n) для відносно рівномірного розподілу – кожен з n вставлень спускається по дереву глибиною ~log n. Дерево будується з нуля на кожному кроці; це простіше, ніж оновлювати його, і, насправді, побудова займає невелику частку від загальної вартості порівняно з перетворенням сил.

insert(node, body):
  if node is empty leaf:            node.body = body; return
  if node is leaf with one body:    subdivide(node); reinsert the old body
  node.mass += body.mass            // aggregates updated on the way down
  node.com  = weighted average of centres of mass
  insert(child_containing(body.pos), body)

Критерій θ

Щоб визначити силу, що діє на тіло, потрібно пройтися деревом від кореня. На кожному вузлі порівнюється довжина сторони вузла s з відстанню d від тіла до центру маси вузла. Якщо вузол достатньо малий і знаходиться достатньо далеко – тобто, якщо

s / d < θ // θ = кут відкриття, зазвичай 0,5 — вважати весь вузол єдиною точковою масою в його центрі маси, додавати цей внесок і не спускатися далі. Інакше відкрити вузол та рекурсивно викликати його дочірні вузли. Лист завжди оцінюється безпосередньо, а тіло ніколи не притягує себе.

force(тіло, вузол): якщо вузол є листком: якщо вузол.тіло !== тіло: повертати pointForce(тіло, вузол.тіло); повернути 0; const s = вузол.розмір, d = відстань(позиція(тіла), позиція(центр_маси_вузла)); якщо s / d < theta повернути pointForce(тіло, вузол); // один терм, ціле піддерево нехай f = 0; для кожного c у node.дітки якщо c: f += force(тіло, c); повернути f Оскільки кожне тіло спускається лише на логарифмічну кількість рівнів до тих пір, поки критерій не буде задоволено, обхід коштує O(log n) на тіло та O(n log n) загалом на крок. Це вся алгоритм.

s / d < θ           // θ = opening angle, typically 0.5

Торговельна угода, яку ви обираєте з θ

θ – це один регулятор між точністю та швидкістю. θ = 0 ніколи не задовольняє критерію, тому кожен вузол відкривається, кожна пара оцінюється і ви повертаєтеся до точного O(n²) прямого підсумовування. Збільшення θ відкриває менше вузлів: менше взаємодійних термінів, швидші кроки та більша похибка апроксимації сили. Значення близько 0,5 є звичайним залишком в астрофізиці, 0,3 для роботи, що вимагає високої точності, і 1,0 для швидкої візуальної демонстрації, де невелике неправильне гравітаційне поле не має значення.

Важливо знати дві речі. По-перше, монопольне (центромасштабне) наближення ігнорує тетрапольний момент комірки, тому похибка зростає швидше, ніж можна було б очікувати, якщо б це було простою оцінкою; деякі коди додають тетрапольний термін, щоб дешево повернути точність. По-друге, помилка Barnes–Hut не є випадковим шумом – вона корелюється з геометрією дерева, тому тіло, що знаходиться біля межі комірки, може бачити коливання сили, поки дерево будується. Стандартною практикою є пом’якшення потенціалу (заміна 1/r² на 1/(r² + ε²)), щоб запобігти вибуху сили, коли дві тіла проходять близько одна до одної, що в іншому випадку вимагало б неможливо малого кроку часу.

Як ця симуляція N-тіла використовує метод Barnes–Hut

Ця симуляція на сайті створює нове квадрантове дерево кожного кадру, обстежує його один раз на тіло з θ = 0.5 та пом’якшенням (softening length) кількох пікселів, а потім просуває стан за допомогою Leapfrog – симплектичного інтегратора, оскільки гравітаційна система, інтегрована з RK4 або явним Еулером, втрачає або надмірно набирає енергії, і галактика помітно стискається або випаровується. Сила та інтегратор обираються разом: немає потреби витрачати чотири обчислення сили на крок на RK4, коли сама сила є оцінкою, контрольованою θ.

на кожному кадрі: 1. побудова квадрантового дерева з поточних положень – O(n log n) 2. накопичення маси + центру ваги вгору – O(n) 3. обстеження кожного тіла з θ = 0.5 – O(n log n) 4. Leapfrog kick–drift–kick – O(n) Для повноти інформації: Barnes–Hut не єдиний спосіб вирватися з O(n²). Метод Швидкої Мультиполя (Fast Multipole Method) розширює поле в мультиполі на обох сторонах взаємодії та досягає O(n) для заданої точності, ціною значно важливої реалізації; методи на основі сітки з використанням FFT вирішують потенціал на сітці і є звичайним вибором для космологічних об’ємів. Для інтерактивної мозаїки з тисячами тіл квадрантове дерево та один добре підібраний θ – це оптимальне рішення.

each frame:
  1. build quadtree from current positions   O(n log n)
  2. accumulate mass + centre of mass upward O(n)
  3. traverse per body with θ = 0.5          O(n log n)
  4. Leapfrog kick–drift–kick                O(n)

Часті запитання

Спробуйте наживо

Усе, що вище, працює прямо у вашому браузері — відкрийте N-Body Gravity і змінюйте параметри під час роботи. Нічого не встановлюється, нічого не завантажується на сервер, уся модель живе в одній вкладці.