Barnes-Hut: Symulacja 100 000 Ciężarowych Ciał w Czasie Rzeczywistym

Ściana O(n²)

Grawitacja to siła działająca między każdą parą ciał. Aby postępować w systemie N-ciał o jeden krok, potrzebujesz przyspieszenia każdego ciała, a bezpośrednie sumowanie daje to dokładnie: dla każdego ciała dodajesz siłę przyciągania od wszystkich pozostałych ciał. To oznacza n(n−1)/2 par – O(n²) pracy na kroku. Tysiąc ciał oznacza pół miliona interakcji na krok i uruchamia się komfortowo z 60 fps w przeglądarce. Sto tysięcy ciał oznacza pięć miliardów interakcji na krok, a teraz czekasz minuty na klatkę.

Algorytm Barnesa–Huta (Josh Barnes i Piet Hut, opublikowany w Nature w 1986 roku) pokonuje tę ścianę dzięki jednej obserwacji fizycznej: odległy grom gwiazd prawie dokładnie działa jak pojedyncza gwiazda w jego środku masy. Nie musisz wiedzieć, że Galaktyka Andromedy zawiera bilion gwiazd, aby obliczyć jej wpływ na Słońce – potrzebujesz tylko jej całkowitej masy i środka masy. Barnes–Hut przekształca tę obserwację w strukturę danych i uzyskuje O(n log n) pracy na kroku.

Drzewo

Rekurencyjnie podział na obszar symulacji. W 2D każdy kwadrat dzieli się na cztery ćwiartki — drzewo quadtree; w 3D każdy sześcian dzieli się na osiem oktantów — drzewo octree. Kontynuuj podział komórki, aż do momentu, gdy zawiera ona maksymalnie jedno ciało. Każdy węzeł wewnętrzny przechowuje dwie agregaty obliczone podczas powrotnej podróży w górę: całkowitą masę ciał poniżej niego oraz ich oś środkowa masy.

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 // agregaty aktualizowane podczas schodzenia się node.com = weighted average of centres of mass insert(child_containing(body.pos), body) Budowanie drzewa kosztuje O(n log n) dla rozsądnie jednorodnego rozkładu — każdy z n wstawień opada po drzewie o głębokości ~log n. Drzewo jest budowane od nowa za każdym krokiem; to jest prostsze niż jego aktualizowanie i, w praktyce, koszt budowy jest niewielki w porównaniu z podróżą siły.

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)

Kryterium θ

Aby uzyskać siłę działającą na jedno ciało, przejdź drzewo od korzenia. W każdym węźle porównaj długość boku węzła s z odległością d od pozycji ciała do środka masy węzła. Jeśli węzeł jest wystarczająco mały i wystarczająco daleko — to znaczy, jeśli

s / d < θ // θ = kąt otwarcia, zwykle 0.5 — traktuj cały węzeł jako pojedynczą masę punktową w jego środku masy, dodaj ten jeden wkład i nie zstępuj dalej. W przeciwnym razie otwórz węzeł i rekurencyjnie przejdź do jego potomków. Liść zawsze jest oceniany bezpośrednio, a ciało nigdy nie przyciąga się samo do siebie.

force(body, node): if node is a leaf: if node.body !== body: return pointForce(body, node.body); return 0; const s = node.size, d = distance(body.pos, node.com); if (s / d < theta) return pointForce(body, node); // jeden wyraz, cały poddrzewo let f = 0; for (const c of node.children) if (c) f += force(body, c); return f Ponieważ każde ciało zstępuje tylko o logarytmiczną liczbę poziomów zanim spełnione zostaje kryterium, koszt przeszukiwania wynosi O(log n) na ciało i O(n log n) w ogóle. To jest cały algorytm.

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

Wybór kompromisu z parametrem θ

θ to pojedynczy regulator pomiędzy dokładnością a szybkością. Wartość θ = 0 nigdy nie spełnia kryterium, dlatego każdy węzeł jest otwierany, każda para jest oceniana i wracamy do dokładnego obliczenia O(n²) sumy bezpośredniej. Zwiększanie θ otwiera mniej węzłów: mniejsza liczba interakcji, szybsze kroki i większy błąd przybliżony w sile. Wartości około 0,5 są standardowo używane w astrofizyce, z 0,3 dla pracy wymagającej wysokiej dokładności i 1,0 dla szybkiego demonstracyjnego wizualizacji, gdzie niewielka nieścisłość pola pływowego jest niezauważalna.

Należy jednak pamiętać o dwóch ważnych kwestiach. Po pierwsze, przybliżenie monopola (centrum masy) ignoruje moment kwadrupolowy komórki, co powoduje wzrost błędu wraz ze zwiększaniem θ szybciej niż sugeruje wstępny oszacowanie; niektóre implementacje dodają termin kwadrupolowy, aby tanio odzyskać dokładność. Po drugie, błąd Barnes–Hut nie jest przypadkowym szumem – jest on skorelowany z geometrią drzewa, co oznacza, że ciało znajdujące się w pobliżu granicy komórki może obserwować drgania siły podczas budowy drzewa. Standardową praktyką jest również zmiękczanie potencjału (zamiana 1/r² na 1/(r² + ε²)), które zapobiega eksplozji siły, gdy dwie ciele znajdują się blisko siebie, co w przeciwnym razie wymagałoby niemożliwie małego kroku czasowego.

Jak symulacja N-ciał tutaj z tym korzysta

Symulacja na tej stronie buduje świeżą kwadryczność co kadrę, przechodzi przez nią raz na ciało z θ = 0.5 i długością zmiękczenia kilku pikseli, a następnie porusza stan z Leapfroga – integratora symplectic, ponieważ system grawitacyjny zintegrowany z RK4 lub Eulera jawnym traci lub zyskuje energię i galaktyka widocznie się rozpada lub paruje. Aproksymacja siły i integrator są wybierane razem: nie ma sensu wydawać czterech ocen sił na krok w RK4, gdy każda siła jest sama w sobie przybliżeniem kontrolowanym przez θ.

co kadrę: 1. buduj kwadryczność z bieżących pozycji O(n log n) 2. akumulacja masy + środek ciężkości do góry O(n) 3. przechodź co ciało z θ = 0.5 O(n log n) 4. Leapfrog kick–drift–kick O(n) Dla kompletności: Barnes–Hut nie jest jedynym ucieczką przed O(n²). Metoda Fast Multipole Method rozszerza pole w wielopole na obu stronach interakcji i osiąga O(n) dla określonej dokładności, kosztem znacznie cięższej implementacji; metody cząstek-siatki rozwiązują dla potencjału na siatce z FFT i są zwykle wybierane do objętości kosmologicznych. Dla interaktywnego płótna z dziesiątami tysięcy ciał, kwadryczność i jeden dobrze dobrany θ to idealne połączenie.

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)

Często zadawane pytania

Wypróbuj na żywo

Wszystko powyżej działa bezpośrednio w Twojej przeglądarce — otwórz N-Body Gravity i zmieniaj parametry podczas działania. Nic nie jest instalowane ani przesyłane na serwer, cały model działa w jednej karcie.