Orbita liczona Verletem: skale, moment pędu i kontrola toru

Potrzebne wiadomości: Grawitacja Newtona, energia i moment pędu. Potrzebna jest wcześniejsza lekcja o metodzie Verleta; znajomość elipsy pomaga w ostatnim zadaniu.

Jedno centrum siły i jawne uproszczenia

Małe ciało porusza się w polu nieruchomej masy centralnej M, przy m znacznie mniejszym od M. Pomijamy atmosferę, napęd, inne ciała, niesferyczność źródła i poprawki relatywistyczne. μ = GM ma jednostkę m³/s². Dla r = (x; y) i r = |r| > 0:

\frac{d^2\mathbf r}{dt^2}=-\frac{\mu\mathbf r}{r^3},\qquad \varepsilon=\frac{v^2}{2}-\frac{\mu}{r},\qquad \ell_z=xv_y-yv_x.

ε to energia na jednostkę masy w J/kg, a ℓ_z moment pędu na jednostkę masy w m²/s. W modelu ciągłym obie są stałe. Potęga r³ w mianowniku wektora zapewnia, że wartość przyspieszenia wynosi μ/r². Wstawienie r² w mianowniku wektora dałoby inne prawo siły.

Jednostki usuwamy przez skalę, nie przez skasowanie etykiet

Wybieramy długość R₀ oraz czas T₀ = √(R₀³/μ). Zapisujemy r = R₀R, t = T₀τ, v = (R₀/T₀)V. Wielkie R i V oraz τ są bezwymiarowe. Dwa różniczkowania po czasie wprowadzają czynnik 1/T₀²:

\frac{d^2\mathbf R}{d\tau^2}=-\frac{\mathbf R}{|\mathbf R|^3},\qquad E^*=\frac{V^2}{2}-\frac{1}{|\mathbf R|},\qquad L^*=X V_y-Y V_x.

Energię właściwą odzyskujemy jako ε = (μ/R₀)E*, a moment pędu jako ℓ_z = √(μR₀)L*. Dla orbity kołowej o promieniu R₀ szybkość skali V₀ = R₀/T₀ = √(μ/R₀) jest szybkością kołową. Okres obiegu wynosi 2πT₀, nie T₀.

Krok centralnego przyspieszenia

W bezwymiarowych jednostkach stosujemy tę samą kolejność co dla sprężyny:

\mathbf V_{n+1/2}=\mathbf V_n-\frac{h}{2}\frac{\mathbf R_n}{|\mathbf R_n|^3},

\mathbf R_{n+1}=\mathbf R_n+h\mathbf V_{n+1/2},

\mathbf V_{n+1}=\mathbf V_{n+1/2}-\frac{h}{2}\frac{\mathbf R_{n+1}}{|\mathbf R_{n+1}|^3}.

Teraz h oznacza przyrost τ, a fizyczny krok jest hT₀. Nie mieszamy go z sekundami. Przyspieszenie obliczamy dopiero po sprawdzeniu, że promień jest dodatni i obliczenie nie weszło w obszar wyłączony z modelu, np. wnętrze ciała centralnego. Osobliwości w R = 0 nie naprawia dopisanie małej liczby do mianownika: zmieniłoby to prawo siły.

Koło jest niezależnym wzorcem

Dla początkowego stanu bezwymiarowego R(0) = (1; 0) oraz V(0) = (0; 1) dokładny tor to R(τ) = (cos τ; sin τ), E* = −1/2 i L* = 1. W kroku h = 0,10 przyspieszenie początkowe wynosi (−1; 0), V_1/2 = (−0,05; 1), a nowa pozycja R₁ = (0,995; 0,100).

Po ponownym wyliczeniu promienia i przyspieszenia otrzymujemy V₁ = (−0,099748134; 0,995000187). Energia E*₁ ≈ −0,499999969 jest bliska, lecz nierówna −1/2. Moment pędu X₁V_y1 − Y₁V_x1 wynosi 1 z dokładnością arytmetyki.

Dlaczego moment pędu przechodzi test niemal idealnie?

Półimpuls dodaje do V wektor równoległy do R, więc iloczyn wektorowy R × V nie zmienia się. Przesunięcie dodaje hV_1/2 do R, a (R + hV_1/2) × V_1/2 = R × V_1/2, bo wektor pomnożony wektorowo przez siebie daje zero. Drugi centralny półimpuls również nie zmienia momentu pędu.

W arytmetyce dokładnej ten argument działa dla każdego kroku, o ile jest on zdefiniowany. Dlatego mały błąd L* jest potrzebnym sprawdzeniem implementacji, lecz sam nie dowodzi dokładności orbity. Energia fizyczna nie ma tu analogicznego dokładnego zachowania. Dla dostatecznie małego stałego kroku w regularnym ruchu jej błąd zwykle oscyluje bez systematycznego dryfu na obserwowanym przedziale; to nie gwarancja dla bliskiego przelotu przez centrum lub dowolnie wielkiego kroku.

Trzy kontrole na tym samym przedziale

Dla opisanej orbity do τ = 20 maksimum względnego błędu energii w węzłach przy h = 0,10 wynosi około 0,0000246282, ale błąd końcowego położenia względem (cos 20; sin 20) wynosi 0,0618719 jednostki długości. Dla h = 0,05 błąd położenia spada do 0,0155147, a dla h = 0,025 do 0,00388153. Zbliżony czynnik cztery odpowiada drugiemu rzędowi.

W tym szczególnym starcie kołowym obserwowany błąd energii może maleć szybciej niż błąd toru. Nie przenosimy tego na ogólne orbity eliptyczne. Kontrola powinna obejmować konkretne wielkości potrzebne do zadania: położenie, czas powrotu, minimalny promień, energię i moment pędu. „Zachowane prawa” nie zastępują porównania siatek.

Dlaczego elipsa stawia większe wymagania?

Przy zbliżaniu do centrum rośnie przyspieszenie, a z zachowania momentu pędu rośnie także składowa poprzeczna prędkości. Krok dobry przy apocentrum może być za duży przy perycentrum. Wstępnym wskaźnikiem jest iloraz h|V|/|R|: informuje, jaką część lokalnego promienia ciało przebywa w jednym kroku, ale nie daje gwarantowanego błędu.

W tej lekcji zachowujemy stały krok i dobieramy go przez zagęszczanie z kontrolą najbliższego przejścia. Prosta zmiana h w zależności od stanu może zniszczyć symetrię czasową i własności geometryczne standardowego schematu. Adaptacyjne całkowanie jest osobną konstrukcją; nie uzyskuje tych zalet automatycznie przez warunek „gdy blisko, dziel h przez dwa”.

Zadania do tego tematu

• Jednostka czasu orbity nie jest jej okresem
• Moment pędu równy jeden nie gwarantuje trafienia w punkt
• Perycentrum ujawnia za duży krok

Dalsza lektura

Materiał referencyjny