Barnes-Hut: Symulowanie Galaktyki Bez Obliczania Każdej Pary Ciał

Problem N-ciałek to mur parzystości

Symulowanie grawitacji dla n ciał w sposób bezpośredni oznacza obliczanie siły między każdą parą: n(n-1)/2 wyliczeń sił na kroku czasowym. Podwojenie liczby ciał powoduje cztery razy większe obciążenie. Ten algorytm sumowania elementów jest dokładny i, dla kilku tysięcy ciał na nowoczesnym GPU, jest nawet praktyczny – ale dla dziesiątek tysięcy a nawet milionów cząstek w symulacji galaktyki lub dużego skupiska grawitacyjnego, złożoność O(n²) staje się niewykonalna. Josh Barnes i Piet Hut opublikowali algorytm, który ją pokonuje w 1986 roku.

Wgląd: odległy kaskadzie przypomina punkt

Newton już wiek temu pokazał, że symetrycznie kulista masa przyciąga zewnętrzne ciała dokładnie tak, jakby cała jej masa była skoncentrowana w jej środku. Barnes-Hut opiera się na przybliżonej wersji tej samej idei: z wystarczająco dużej odległości cały kaskadzie ciał jest grawitacyjnie niezróżnicowany od jednego punktu masy, który znajduje się w środkowym punkcie kaskady. Zamiast sumować wpływ każdego indywidualnego ciała w odległej kaskadzie, sumuj je raz z wyprzedzeniem i traktuj kaskadę jako jedno ciało dla wszystkiego, co jest daleko.

Aby zdecydować, co stanowi "kaskadę" a co stanowi „wystarczająco blisko”, Barnes-Hut buduje drzewo kwadratowe (drzewo ósemkowe w 3D) we wszystkich ciałach na początku każdego kroku czasowego. Każdy węzeł przechowuje całkowitą masę i środek masy każdego ciała pod nim — obliczany od dołu do góry w jednym przejściu po zbudowaniu drzewa.

Test theta

Obliczanie siły działającej na dany obiekt polega na przechodzeniu przez drzewo od korzenia i, w każdym węźle, przeprowadzaniu jednego testu: porównywaniu rozmiaru węzła (s) z jego odległością (d) od badanego obiektu.

force(body, node): s = node.width // Rozmiar obszaru, którego obejmuje dany węzeł d = distance(body, node.centerOfMass) if s / d < theta OR węzeł jest liściem z jednym innym obiektem: return gravity(body, node.totalMass, node.centerOfMass) // Traktuj jako jeden punkt else: return sum(force(body, child) for child in node.children) // Sprawdź dokładniej Jeśli s/d jest małe – węzeł jest mały w porównaniu z jego odległością – cała poddrzewa jest skumulowana jako pojedynczy punkt-masa i rekurencja się zatrzymuje, niezależnie od tego, ile obiektów zawiera dany węzeł. Jeśli s/d jest duże, węzeł jest albo zbyt blisko lub za duży, aby traktować go jako jeden punkt, a algorytm przechodzi do jego czterech potomków w celu uzyskania dokładniejszej odpowiedzi. Próg theta (θ) to jedynny regulator, który kontroluje cały kompromis: θ = 0 zmusza algorytm do przechodzenia przez wszystkie węzły aż po pojedyncze obiekty, odzyskiwając dokładne O(n²) sumowanie; θ w zakresie około 0.5 do 1.0 jest typowym zakresem pracy, oferując niewielką, kontrolowaną ilość dokładności w zamian za ogromny wzrost prędkości.

force(body, node):
    s = node.width               // size of the region this node covers
    d = distance(body, node.centerOfMass)
    if s / d < theta OR node is a leaf with one other body:
        return gravity(body, node.totalMass, node.centerOfMass)   // treat as 1 point
    else:
        return sum(force(body, child) for child in node.children) // look closer

Dlaczego to jest O(n log n)

Budowa drzewa nad zbiorami n ciał kosztuje O(n log n), tak jak w przypadku każdej budowy zbalansowanego drzewa przestrzennego. Obliczanie siły działającej na pojedyncze ciało, które jest w kontakcie, średnio dotyka O(log n) węzłów drzewa: większość przestrzeni znajduje się daleko i jest podsumowana w jednym węźle, podczas gdy tylko kilka naprawdę blisko położonych ciał zmusza rekursję do głębszego przeszukiwania. Pomnożenie O(log n) dla każdego ciała przez n ciał daje całkowity koszt obliczeń siły w jednym kroku czasowym, czyli O(n log n) – to samo klasa asymptotyczna jak sortowanie porównawcze i znacznie lepsze niż O(n²) wynikające z bezpośredniego sumowania, gdy n osiągnie dziesiątki tysięcy.

Jak Barnes-Hut wpływa na wyniki

Przyspieszenie nie jest bezinteresowne. Przybliżenie środka masy wprowadza rzeczywisty błąd siły, gdy θ jest większe od zera, a ten błąd jest największy w przypadku bliskich spotkań między grupami o podobnej wielkości — dokładnie takie sytuacje, jak bliskość kolizji dwóch gromad gwiazd, gdzie często liczą się precyzja. Kody astrofizyczne produkcyjne rozwiązują to poprzez rozszerzenie wielopołowe (dodawanie momentu kwadrupolowego i wyżej, nie tylko masy całkowitej i środka masy) aby zmniejszyć błąd dla danego θ, oraz przez adaptacyjne θ, które automatycznie się kurczy w obszarach o dużej gęstości. Drzewo musi być również całkowicie zbudowane (lub starannie zaktualizowane) za każdym krokiem czasowym, co samo w sobie jest dominującym kosztem w wielu implementacjach — dobry powód, dla którego konstrukcje drzew przestrzennych otrzymują tak dużo uwagi w bazach kodu N-ciagów.

Frequently asked questions

Wypróbuj na żywo

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