Grawitacja N ciał: symulowanie kosmosu w przeglądarce

Problem N ciał to jeden z najstarszych problemów fizyki: mając N mas punktowych oddziałujących grawitacyjnie, przewidzieć ich przyszłe położenia. Dla N = 2 Newton rozwiązał go dokładnie już w XVII wieku — problem dwóch ciał ma rozwiązanie w postaci zamkniętej za pomocą elips, parabol i hiperbol. Dla N = 3 Poincaré udowodnił w 1889 roku, że nie istnieje ogólne rozwiązanie w postaci zamkniętej. Dla N = 10 000 — skali naszej symulacji galaktyki — potrzebny jest wydajny algorytm i staranne sztuczki numeryczne.

Podejście brutalnej siły: O(N²)

Najprostszym sposobem symulowania N grawitujących ciał jest obliczenie dla każdego ciała siły grawitacyjnej pochodzącej od każdego innego ciała. Prawo grawitacji Newtona daje siłę między ciałami i oraz j jako:

F_ij = G m_i m_j / |r_i − r_j|² w kierunku od i do j.

Dla każdego z N ciał sumujesz N − 1 sił parami. To O(N²) obliczeń siły na krok czasowy. Dla N = 1000 to milion obliczeń na klatkę. Dla N = 10 000 — sto milionów. Przy 60 klatkach na sekundę symulacja metodą brutalnej siły dla 10 000 ciał wymaga sześciu miliardów obliczeń siły na sekundę — zbyt wolno dla czasu rzeczywistego w przeglądarce.

Brutalna siła jest jednak przydatna dla małego N — poniżej około 500 ciał — gdzie jest maksymalnie dokładna i prosta w implementacji. Nasza symulacja wykorzystuje brutalną siłę dla małej liczby cząstek i przełącza się na Barnes-Hut powyżej pewnego progu.

Barnes-Hut: O(N log N) z drzewem czwórkowym

Algorytm Barnes-Hut, opublikowany przez Josha Barnesa i Pieta Huta w 1986 roku, redukuje złożoność do O(N log N), wykorzystując kluczowe przybliżenie: odległe grupy ciał można traktować jako pojedyncze ciało w ich środku masy.

Algorytm buduje drzewo czwórkowe (w 2D) lub ósemkowe (w 3D), które rekurencyjnie dzieli przestrzeń na komórki. Każda komórka przechowuje sumaryczną masę i środek masy wszystkich ciał wewnątrz niej. Aby obliczyć siłę działającą na ciało i, przechodzisz przez drzewo. Dla każdego węzła: jeśli węzeł jest wystarczająco daleko w stosunku do swojego rozmiaru — konkretnie, jeśli szerokość węzła s podzielona przez odległość do niego d spełnia s/d < θ (parametr strojenia, zwykle 0,5–1,0) — traktujesz cały węzeł jako pojedyncze ciało w jego środku masy. W przeciwnym razie rekurencyjnie schodzisz do dzieci węzła.

Algorytm	Złożoność	N = 1000	N = 10 000
Brutalna siła	O(N²)	1M operacji	100M operacji
Barnes-Hut	O(N log N)	~10 tys. operacji	~133 tys. operacji
Fast Multipole	O(N)	~1 tys. operacji	~10 tys. operacji

Parametr θ kontroluje kompromis między dokładnością a szybkością. θ = 0 daje dokładne wyniki brutalnej siły. θ = 1,0 daje znaczne błędy przybliżenia. θ = 0,5 zapewnia dobrą równowagę — błędy siły poniżej 1% dla większości ciał przy zachowaniu interaktywności symulacji.

Zmiękczanie grawitacyjne: zapobieganie osobliwościom

Prawo grawitacji Newtona ma osobliwość: gdy dwa ciała bardzo się do siebie zbliżają, siła między nimi rośnie bez ograniczeń (1/r² → nieskończoność, gdy r → 0). W prawdziwej galaktyce rzeczywiste bliskie spotkania gwiazd są niezwykle rzadkie, ponieważ gwiazdy są znikome w porównaniu z odległościami między nimi. W symulacji cząstki mogą zbliżać się dowolnie, powodując numeryczne eksplozje — jeden krok czasowy widzi ogromną siłę, generuje ogromną prędkość, a cząstka odlatuje w nieskończoność.

Standardowym rozwiązaniem jest zmiękczanie grawitacyjne: zastąpienie mianownika r² wyrażeniem r² + ε², gdzie ε to długość zmiękczenia. To ogranicza maksymalną siłę do wartości skończonej, gdy dwie cząstki znajdują się w odległości ε od siebie. Siła wygląda teraz tak:

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

Fizycznie jest to równoważne rozmazaniu każdej masy punktowej w kulę o promieniu ε. Dla wartości r znacznie większych niż ε siła jest nieodróżnialna od masy punktowej. Dla bardzo małych r siła jest ograniczona do G m_i m_j / ε².

Dobór ε to sztuka. Zbyt duża wartość tłumi rzeczywiste efekty dynamiczne — cząstki, które powinny się wzajemnie „procować", zamiast tego gładko przez siebie przechodzą. Zbyt mała wartość — i bliskie spotkania zaczynają dominować nad krokiem czasowym, drastycznie spowalniając symulację. Typowym wyborem dla symulacji galaktyki jest ε równe około 1% średniego odstępu między cząstkami.

Całkowanie czasowe: metoda leapfrog

Mając siły, jak przesuwać pozycje i prędkości cząstek w czasie? Proste całkowanie Eulera — aktualizacja prędkości, a potem pozycji, przy użyciu bieżącej siły — jest szybkie, ale z czasem traci energię, przez co orbity powoli rozkręcają się na zewnątrz. Dla symulacji N ciał uruchomionej przez tysiące kroków czasowych ten dryf energii niszczy wyniki.

Integrator leapfrog, zwany też metodą Störmera-Verleta, rozwiązuje ten problem za pomocą przesuniętej aktualizacji. Prędkości są obliczane w krokach czasowych o połówkowej wartości, podczas gdy pozycje — w krokach całkowitych:

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

Ta drobna zmiana — prędkości i pozycje nie są już aktualizowane w tym samym momencie — sprawia, że integrator jest symplektyczny: dokładnie zachowuje nieznacznie zmodyfikowany hamiltonian (funkcję energii) układu. W praktyce oznacza to, że całkowita energia oscyluje wokół prawdziwej wartości, zamiast dryfować — orbity pozostają stabilne przez miliony kroków czasowych, nie stając się niefizyczne.

Dla jeszcze lepszej dokładności w wysokoprecyzyjnych symulacjach astronomicznych wykorzystuje się symplektyczne integratory wyższego rzędu, takie jak metoda czwartego rzędu Yoshidy czy klasyczna metoda Rungego-Kutty czwartego rzędu (niesymplektyczna, ale bardzo dokładna na każdym kroku). Nasza symulacja przeglądarkowa wykorzystuje leapfrog ze względu na jego prostotę, stabilność i doskonałe zachowanie energii.

Formowanie się galaktyk i struktura emergentna

Uruchom wystarczającą liczbę grawitujących cząstek z odpowiednimi warunkami początkowymi, a spontanicznie pojawi się niezwykła struktura. Nadaj chmurze cząstek pewien początkowy moment pędu i pozwól grawitacji ją zapaść: chmura spłaszcza się w dysk, pod wpływem niestabilności Jeansa powstają zgrupowania, a różnicowa rotacja generuje ramiona spiralne. To nie są zaprogramowane cechy — powstają z N prostych cząstek podążających za prawem Newtona.

Prawdziwe symulacje kosmologiczne, takie jak projekty Illustris i EAGLE, wykorzystują te same fundamentalne algorytmy — Barnes-Hut lub Fast Multipole dla grawitacji, SPH lub ruchomą siatkę dla dynamiki gazu — skalowane do miliardów cząstek na klastrach superkomputerowych. Jakościowe zachowania, które generują, zaskakująco dobrze pasują do obserwowanej morfologii galaktyk. Nasza symulacja przeglądarkowa jest miniaturową wersją tej samej fizyki obliczeniowej, której zawodowi kosmolodzy używają do zrozumienia wielkoskalowej struktury wszechświata.