Pół szklanki mocnego kodu

Koniec świata

W czasach niepewności, wielkich zmian, kryzysów ludzie więcej myślą o sprawach ostatecznych. Czy to nie zbliża się koniec cywilizacji, a może nawet całego świata? W dawnych czasach dobrym wzorem ładu i uporządkowania był Kosmos w dostatecznie dużej skali. No bo przecież nie Ziemia, ze swoją przyziemną nieprzewidywalnością - ale już jej wspólna podróż z Księżycem wokół Słońca, od "zawsze" taka sama, mogłaby stanowić jakiś punkt odniesienia. Albo jeszcze lepiej: popatrzmy na cały Układ Słoneczny! Czy jego leniwie przemierzające przestrzeń planety są oazą spokoju, przewidywalności i stabilności - jeśli nie na wieczność, to może przynajmniej na miliony, lub lepiej miliardy, lat?

Jak trudno rozstrzygnąć to pytanie na gruncie matematyki, przekonał się sam wielki Henri Poincaré. Ale od czego są komputery... zwłaszcza że praktycznie sprawa jest bardzo prosta: należy zbadać trajektorie ruchu planet w interesującym nas okresie. Skoro planety poruszają się w próżni, na ich ruch ma wpływ jedynie siła grawitacji ze strony pozostałych ciał: przede wszystkim Słońca, którego masa jest 750 razy większa od łącznej masy wszystkich planet.

W najprostszej sytuacji - gdy jest tylko Słońce, o masie $m1,$ i jedna planeta, o masie $|m2,$ - przyciągają się one z siłą grawitacji o wartości

m1⋅m2 r2, F = G

działającą wzdłuż łączącego je promienia długości $r.$ W ogólnym przypadku, gdy mamy do czynienia z zagadnieniem $|N$ ciał, jest analogicznie: wypadkowa siła działająca na $i$ -te ciało będzie sumą sił przyciągania go przez pozostałe.

Ponieważ, jak pamiętamy z lekcji fizyki, przyspieszenie jest proporcjonalne do siły, $|⃗F = m$ a przecież przyspieszenie $|⃗a = ∆⃗v- ∆t$ i podobnie $|⃗v = ∆-⃗x , ∆t$ to mnożąc te równości przez $∆ t$ dostajemy układ równań

⎧ ⎪⎪∆x⃗i = ∆t⋅⃗vi, . ⎨⎪⎪∆⃗vi =∆ t⋅P F⃗i j/mi,dla i = 1,...,N ⎩ j xi

Ponieważ $∆ ⃗x = ⃗xnowe− ⃗xstare$ i podobnie $∆ ⃗v =⃗ vnowe− ⃗vstare,$ a siła przyciągania $|⃗Fi j$ pomiędzy ciałem $i$ -tym a $| j$ -tym zależy właśnie od $|⃗xi$ i $x⃗j ,$ to powyższe jest zamkniętym układem $|2N$ równań (wektorowych w przestrzeni trójwymiarowej), które wystarczy rozwiązać, na przykład takim oto algorytmem:

⎧ nowe stare stare ⎪⎪⎨⃗xi = ⃗xi +∆ t⋅⃗vi , ⎪⎪⃗vniowe= ⃗vstiare+ ∆t ⋅Pj xi⃗Fsita jre/mi. ⎩

Rzeczywiście, znając stare wartości prędkości i położenia wszystkich ciał, możemy na ich podstawie wyznaczyć z powyższego wartości nowe, tzn. po czasie $∆ t.$ Okazuje się, że wyprowadzona przez nas metoda to nic innego, jak znany i popularny schemat Eulera. O tym, że wspaniale radzi sobie np. w modelowaniu chorób zakaźnych, można było przeczytać w Delcie 5/2018.

Naszym celem będzie teraz zbadanie, czy Układ Słoneczny nie rozpadnie się przez, powiedzmy, najbliższe 10 tysięcy lat - w ludzkiej skali to dostatecznie długi czas. Masy oraz dzisiejsze początkowe położenia i prędkości wszystkich obiektów pobierzemy ze strony internetowej NASA, dotyczącej projektu Horizons. Przy tak długim czasie symulacji użyjemy niesamowicie krótkiego kroku czasowego, $∆ t = 2 godz.,$ i bardzo szybkiego komputera. Fani rubryki o mocnym kodzie wybaczą mi, że jedynie w pseudokodzie naszkicuję sedno powyższego algorytmu (prawdziwą implementację, z której pochodzą obrazki, zrobiłem w MATLAB-ie):

Co wynika z przeprowadzonych symulacji? Wieści są nad wyraz niepokojące. Z wykresu obok wynika, że jeszcze w obecnym stuleciu Merkury zostanie wytrącony ze swojej orbity!

Co gorsza, symulacja przewiduje, że podobny los czeka Wenus i Ziemię: zaledwie za 30 lat promień naszej orbity wzrośnie dwa razy - co spowoduje, że będzie do nas docierać 8 razy mniej energii ze Słońca. Nic nie pomoże, nawet globalne ocieplenie: nasza planeta zacznie zamarzać. Za mniej niż 500 lat Ziemia znajdzie się 10 razy dalej od Słońca, a to już oznacza całkowity i nieodwołalny koniec... Tylko odległe planety-olbrzymy: Jowisz, Saturn, Uran i Neptun przetrwają ten kataklizm.

Czy to prawda? Czy rzeczywiście jeszcze za naszego życia czekają nas tak dramatyczne wydarzenia? Dlaczego tak musi się stać? Przyjrzyjmy się otrzymanym rozwiązaniom raz jeszcze, ale z innej strony - i zastanówmy się, czy czasem nie zostaliśmy oszukani. Na wykresach na następnej stronie możemy zobaczyć, jak bardzo wraz z upływem czasu całkowita energia układu odchyla się od wartości początkowej.

Okazuje się, że schemat Eulera, który przecież tak dobrze sprawdza się w wielu innych sytuacjach, tutaj napotyka poważny problem: w sposób sztuczny pompuje do układu energię, której już po kilkudziesięciu latach przybywa kilka procent! Jasne więc, że wiele ciał tego nie wytrzymuje i urywa się z orbit - tyle, że jest to wyłącznie artefakt symulacji: w końcu prawdziwe planety zasada zachowania energii przecież obowiązuje. Problem nierespektowania fizycznych zasad zachowania dotyczy nie tylko schematu Eulera, ale też wielu innych popularnych metod numerycznych wyznaczania trajektorii układów dynamicznych. Na przykład znana i lubiana metoda RK4, dużo dokładniejsza od Eulera, też cierpi na podobną chorobę: w niej z kolei energia powoli, ale nieubłaganie... znika.

Względna zmiana energii całkowitej układu w czasie. Od lewej: schemat Eulera (1% już po 50 latach), symplektyczny schemat Eulera ( $10 3$ przez 50 tysięcy lat), schemat Verleta ( $10 4$ przez 50 tysięcy lat)

Z drugiej strony, zgodnie z teorią zbieżności schematu Eulera, wystarczy dostatecznie zmniejszyć $∆ t,$ by z dowolną zadaną dokładnością uzyskać rozwiązanie. Jednak w moim przypadku zmniejszenie $∆ t$ do minut lub - być może - sekund nie wchodzi w grę z prozaicznych powodów: nie mam aż tak szybkiego komputera!

Jest jednak promyk nadziei: można ulepszyć schemat Eulera jednym subtelnym ruchem ręki - i nie jest to gest sięgania po kartę kredytową w celu kupna nowego komputera. Wystarczy po prostu zmienić kolejność dwóch instrukcji w pętli, o tak:

Wycinek zmian odległości od Słońca względem odległości początkowej przez najbliższe **50 tysięcy** lat, obliczony schematem Verleta z $∆ t 10$ dni. Merkury - linia kropkowana, Ziemia - pomarańczowa, Jowisz - czarna. Tym razem wszystko przebiega cyklicznie, jak w zegarku

Jak widać na środkowym wykresie powyżej, tym razem energia nie jest co prawda idealnie zachowana, ale przynajmniej oscyluje wokół pierwotnej wartości. Orbity uzyskane w ten sposób przez tysiące lat trzymają fason i nie rozwijają się w dramatyczne spirale. Jeszcze lepsza - i jednocześnie dokładniejsza - jest metoda leap-frog, zwana też schematem Verleta.

Wykres obok pokazuje, że w jej przypadku odchylenie nie przekroczyło jednej dziesięciotysięcznej na przestrzeni 50 tysięcy lat!

I tym sposobem, po zmianie schematu na symplektyczny schemat Eulera (a tym bardziej na metodę leap-frog), Świat został uratowany: przez kolejne dziesiątki tysięcy lat - jak na rysunku obok - ruch planet jest uspokajająco cykliczny, jak w zegarku. Układ Słoneczny trzyma się stabilnie!

Może więc nie będzie tak źle?