Symulacja grawitacji N-ciał: jak tańczą planety

Prawo powszechnego ciążenia Newtona

Każda masa we wszechświecie przyciąga każdą inną masę siłą daną prawem Newtona:

F = G × m₁ × m₂ / r²

Gdzie G to stała grawitacyjna (6,674 × 10⁻¹¹ N m² kg⁻²), m₁ i m₂ to masy dwóch ciał, a r to odległość między ich środkami. Siła jest przyciągająca i działa wzdłuż linii łączącej oba ciała.

Dla dwóch ciał daje to system w pełni rozwiązywalny. Ciała orbitują wokół wspólnego środka masy po elipsach — prawa Keplera wynikają bezpośrednio. Dodaj trzecie ciało i wszystko się zmienia.

Dlaczego N ciał jest trudne: problem trzech ciał

W 1887 roku Henri Poincaré udowodnił, że grawitacyjny problem trzech ciał nie ma ogólnego rozwiązania w postaci zamkniętej. Nie da się zapisać równania dającego pozycję trzech grawitujących mas w czasie t na podstawie ich warunków początkowych. System jest z reguły chaotyczny — małe różnice w warunkach początkowych prowadzą do zupełnie różnych długoterminowych trajektorii.

To nie jest porażka matematyki. To fundamentalna właściwość równań. Trzy ciała oddziałujące grawitacyjnie mogą wykazywać każdy rodzaj zachowania: stabilne orbity, ruch quasi-okresowy, rezonanse, chaotyczne błądzenie i wyrzucenie jednego ciała w nieskończoność.

Jedynym praktycznym sposobem przewidywania ewolucji N ciał jest całkowanie numeryczne: wykonuj małe kroki czasowe, oblicz siły, zaktualizuj prędkości i pozycje, powtarzaj. To właśnie robi symulacja — i to, co robi każdy kod astrofizyczny, od modeli układu słonecznego po symulacje galaktyk.

Naiwny algorytm: złożoność O(n²)

Prostym podejściem do obliczania sił grawitacyjnych jest sprawdzenie każdej pary ciał:

dla każdego ciała i:
    siła[i] = (0, 0, 0)
    dla każdego ciała j ≠ i:
        r = pozycja[j] - pozycja[i]
        odl = długość(r)
        siła[i] += G * masa[i] * masa[j] / odl² * normalizuj(r)

To wymaga n × (n−1) obliczeń siły na krok czasowy — złożoność O(n²). Dla 100 ciał: 9900 obliczeń. Dla 1000 ciał: 999 000 obliczeń. Dla 10 000 ciał: prawie 100 milionów obliczeń. Dla prawdziwej galaktyki ze 100 miliardami gwiazd naiwny algorytm jest całkowicie niewykonalny.

Softening zapobiegające osobliwościom

W równaniu siły czai się problem numeryczny: gdy dwa ciała zbliżają się bardzo blisko, r dąży do zera, a siła dąży do nieskończoności. W rzeczywistości gwiazdy nie są masami punktowymi — bliskie spotkania wiążą się ze złożoną fizyką. W symulacji dodajemy długość softeningu ε:

F = G × m₁ × m₂ / (r² + ε²)

To ogranicza maksymalną siłę i zapobiega numerycznym eksplozjom podczas bliskich spotkań. Długość softeningu jest parametrem symulacji — zbyt duża i fizyka staje się nierealistyczna; zbyt mała i bliskie spotkania destabilizują całkowanie.

Barnes-Hut: redukcja do O(n log n)

W 1986 roku Josh Barnes i Piet Hut opublikowali algorytm, który redukuje złożoność N-ciał z O(n²) do O(n log n) — transformacyjną poprawę. Kluczowy wgląd: odległy klaster wielu ciał można przybliżyć jako pojedyncze ciało w środku masy klastra. Błąd grawitacyjny wynikający z tego przybliżenia jest mały, gdy klaster jest daleko.

Algorytm Barnesa-Huta buduje oktadrzewo (drzewo, które rekurencyjnie dzieli przestrzeń 3D na osiem oktantów w każdym węźle). Każdy węzeł wewnętrzny przechowuje całkowitą masę i środek masy wszystkich ciał w swoim poddrzewie.

Przy obliczaniu siły na ciało i drzewo jest przeszukiwane. W każdym węźle podejmowana jest decyzja na podstawie stosunku s/d — gdzie s to rozmiar przestrzenny węzła, a d to odległość od ciała i do środka masy węzła. Jeśli s/d < θ (próg kąta otwarcia, zwykle 0,5 do 1,0), węzeł jest traktowany jako pojedyncza masa punktowa. W przeciwnym razie dzieci węzła są otwierane i oceniane indywidualnie.

function oblicz_siłę(ciało i, węzeł):
    jeśli węzeł jest liściem:
        dodaj bezpośrednią siłę od ciała węzła
    inaczej:
        s = rozmiar węzła
        d = odległość od i do środka masy węzła
        jeśli s/d < θ:
            dodaj przybliżoną siłę od masy węzła
        inaczej:
            dla każdego dziecka węzła:
                oblicz_siłę(ciało i, dziecko)

Przy θ ≈ 0,5 błąd w oszacowaniach siły jest poniżej 1%, podczas gdy obliczenia spadają z godzin do sekund przy tej samej liczbie cząstek. Nowoczesne symulacje kosmologiczne używają wariantów Barnesa-Huta do symulowania miliardów cząstek ciemnej materii.

Całkowanie leapfrog: stabilność ponad dokładność

Znajomość sił to tylko połowa problemu. Trzeba też scałkować równania ruchu — zaktualizować pozycje i prędkości — w sposób, który pozostaje stabilny przez wiele tysięcy kroków czasowych.

Proste całkowanie Eulera (x += v*dt; v += a*dt) wydaje się naturalne, ale ma fatalną wadę dla mechaniki orbitalnej: nie zachowuje energii. Całkowita energia systemu dryfuje w górę z czasem, a orbity powoli rozszerzają się spiralnie. Przy długich symulacjach planety uciekają swoim gwiazdom.

Całkowatr leapfrog rozwiązuje ten problem, przesuwając aktualizacje pozycji i prędkości o pół kroku czasowego:

// forma kick-drift-kick (KDK):
v_pół = v + a * (dt/2)      // prędkość w pół kroku
x = x + v_pół * dt           // pozycja w pełnym kroku
a = oblicz_siły(x)           // siły w nowej pozycji
v = v_pół + a * (dt/2)       // druga połowa kroku prędkości

Leapfrog jest symplektyczny — dokładnie zachowuje zmodyfikowany hamiltonian bliski prawdziwemu. W praktyce oznacza to, że energia oscyluje wokół poprawnej wartości, zamiast dryfować, pozwalając symulacjom działać przez miliardy kroków czasowych bez akumulacji błędu.

Całkowatory wyższego rzędu (Runge-Kutta 4, Yoshida 6. rzędu) są dokładniejsze na krok, ale wymagają więcej obliczeń siły. Dla grawitacyjnego N-ciał leapfrog z małym krokiem czasowym jest zwykle najlepszym kompromisem.

Prawdziwe zastosowania astrofizyczne

Symulacja N-ciał to podstawowe narzędzie współczesnej astrofizyki:

• Formowanie się systemów planetarnych — symulowanie wczesnej historii układu słonecznego, gdy setki planetozymali zderzały się i łączyły przez miliony lat, by uformować planety.
• Kolizje galaktyk — Droga Mleczna zderzy się z galaktyką Andromedy za około 4,5 miliarda lat. Symulacje N-ciał od lat 70. przewidziały harmonogram i strukturę tego połączenia.
• Gromady kuliste — gęste sferyczne zbiory od 100 000 do 1 miliona gwiazd, symulowane, by zrozumieć dynamikę gwiazdową, formowanie się czarnych dziur i tempo ucieczki.
• Struktura wielkoskalowa — symulacje kosmologiczne, takie jak Millennium i IllustrisTNG, śledzą miliardy cząstek ciemnej materii od Wielkiego Wybuchu do teraźniejszości, odtwarzając kosmiczną sieć włókien, pustek i gromad galaktyk.
• Dynamika asteroid — trajektorie asteroid bliskich Ziemi są obliczane za pomocą całkowania N-ciał, uwzględniającego perturbacje od wszystkich ośmiu planet, Księżyca i Plutona.

Stabilność naszego układu słonecznego

Czy nasz własny układ słoneczny jest stabilny? Okazuje się, że jest to zaskakująco trudne pytanie. Całkowania numeryczne uruchamiane na miliardy symulowanych lat pokazują, że wewnętrzne planety (od Merkurego do Marsa) mają niewielkie, ale niezerowe prawdopodobieństwo wejścia w chaotyczne orbity. Orbita Merkurego mogłaby stać się niestabilna w skalach czasowych kilku miliardów lat, potencjalnie zderzając się z Wenus lub zostając wyrzuconą.

Planety zewnętrzne są bardziej stabilne, ale nie idealnie. Czterociałowy rezonans między Jowiszem, Saturnem, Uranem i Neptunem oznacza, że małe perturbacje mogą się kumulować. Układ słoneczny, jaki znamy, może być jedną z wielu możliwych konfiguracji, które akurat przetrwały wystarczająco długo, by ewoluowało złożone życie.