Autorem poniższego opracowania jest dr Piotr A. Dybczyński z Instytutu Obserwatorium Astronomiczne UAM w Poznaniu.
Całkowanie równań ruchu niezaburzonego (keplerowskiego)
Równania nieperturbowanego ruchu ciała o zaniedbywalnej masie (planetoidy lub komety) wokół Słońca można zapisać następująco:
(1) `d^2/dt^2 vec r = -muvec r/r^3`
gdzie `vec r` jest heliocentrycznym wektorem położenia tego ciała a stała `mu = k^2` jest iloczynem masy Słońca (tu równej 1) i kwadratu tzw. stałej Gaussa: `k = 0.01720209895`. Używamy tu klasycznych jednostek: AU, doba i M⊙ .
Zapis tej jest równoważny układowi trzech równań skalarnych:
`d^2/dt^2 x = -mu*x/r^3`
(2) `d^2/dt^2 y = -mu*y/r^3`
`d^2/dt^2 z = -mu*z/r^3`
gdzie przez `x`, `y` i `z` oznaczyliśmy składowe wektora `vec r`. W mechanice nieba często drugą pochodną po czasie oznaczamy przy pomocy dwóch kropek nad zmienną, co daje nam taki zapis powyższych róznań:
`ddot x = -mu*x/r^3`
(3) `ddot y = -mu*y/r^3`
`ddot z = -mu*z/r^3`
Pierwsza z opisywanych tu metod całkowania numerycznego (Metoda Eulera) stosowana jest bezpośrednio do równań różniczkowych pierwszego rzędu a równania powyższe są układem trzech równań różniczkowych rzędu drugiego. Dla tej metody przepiszemy je więc jako układ sześciu równań różniczkowych rzędu pierwszego:
`dot x = u`
`dot y = v`
`dot z = w`
(4) `dot u = -mu*x/r^3`
`dot v = -mu*y/r^3`
`dot w = -mu*z/r^3`
gdzie przez `u`, `v` i `w` oznaczyliśmy składowe wektora prędkości naszego ciała (`dot vec r`) a pojedyncza kropka nad zmienną oznacza pierwszą pochodną względem czasu `(d/dt)`.
Metoda Eulera
Jest to chyba najprostsza metoda całkowania numerycznego równań różniczkowych. Jeżeli mamy znaleźć rozwiązanie równania różniczkowego w postaci:
(5) `dot x(t) = f(t,x(t))`
w momencie czasu `t_k` mając dane początkowe w momencie `t_0`: `x(0)` oraz `dot x(0)` to dzielimy interwał `(t_0,t_k)` na `n` odpowiednio małych kroków `h` i postępujemy według schematu:
(6) ` x(t+h) = x(t) + h*f(t,x(t))`
zaczynając od momentu startowego `t_0` a kończąc w momencie `t_k`.
Jeżeli każde z równań układu (4) potraktujemy jako równanie postaci (5) i wypiszemy dla niego odpowiednią formułę obliczeniową w oparciu o schemat (6) to dostaniemy komplet formuł do obliczeń:
` x(t+h) = x(t) + h*u(t)`
` y(t+h) = y(t) + h*v(t)`
` z(t+h) = z(t) + h*w(t)`
(7) ` u(t+h) = u(t) - (mu*h)/(r(t)^3)*x(t)`
` v(t+h) = v(t) - (mu*h)/(r(t)^3)*y(t)`
` w(t+h) = w(t) - (mu*h)/(r(t)^3)*z(t)`
Widać, że znając składowe położenia i prędkości w chwili `t` znajdujemy z równań układu (7) wartości tych sześciu składowych na moment `(t+h)` - a to właśnie nazywa się całkowaniem numerycznym. Powtarzamy takie obliczenia n-krotnie aby z czasem "dojechać" do momentu `t_k`.
Metoda Leap Frog (wariant synchroniczny)
Jest to bardzo ciekawa i niemal równie prosta metoda całkowania numerycznego. Mimo iż na każdym kroku całkowania prawą stronę równań ruchu oblicza się tylko raz (jak w metodzie Eulera) to jest to metoda rzędu drugiego i co więcej "częściowo symplektyczna". Częściowo, bo jakkolwiek nie zachowuje ściśle całki energii ale powoduje, że wartość energii całkowitej oscyluje wokół wartości początkowej, zamiast "rozjeżdżać" się dowolnie jak to ma miejsce w metodzie Eulera.
Do zadania wybrałem wariant synchroniczny, najłatwiejszy chyba do zaprogramowania. Jeżeli mamy do czynienia z równaniem różniczkowym drugiego rzędu w ogólnej postaci (prawa strona nie może tu zależeć od prędkości):
(8) `ddot x(t) = f(t,x(t))`
to schemat obliczeniowy wygląda następująco:
(9) `x(t+h) = x(t) + h*dotx(t) + 1/2*h^2*f(t,x(t))`
(10) `dotx(t+h) = dotx(t) + 1/2*h*( (f(t,x(t)) + f(t+h,x(t+h)) )`
Podstawiając (8) do (9) i (10) możemy schemat zapisać nieco czytelniej:
(11) `x(t+h) = x(t) + h*dotx(t) + 1/2*h^2*ddotx(t)`
(12) `dotx(t+h) = dotx(t) + 1/2*h*( ddotx(t) + ddotx(t+h))`
Zastosowanie tego schematu do równań (3) daje, podobnie jak w metodzie Eulera, schemat obliczeniowy składający się z sześciu równań (używamy tu ponownie `u`, `v` i `w` na oznaczenie składowych prędkości):
`x(t+h) = x(t) + h*u(t) - (mu*h^2*x(t))/(2*r(t)^3)`
`y(t+h) = y(t) + h*v(t) - (mu*h^2*y(t))/(2*r(t)^3)`
`z(t+h) = z(t) + h*w(t) - (mu*h^2*z(t))/(2*r(t)^3)`
(13) `u(t+h) = u(t) - (mu*h)/2*( (x(t))/(r(t)^3) + (x(t+h))/(r(t+h)^3))`
`v(t+h) = v(t) - (mu*h)/2*( (y(t))/(r(t)^3) + (y(t+h))/(r(t+h)^3))`
`w(t+h) = w(t) - (mu*h)/2*( (z(t))/(r(t)^3) + (z(t+h))/(r(t+h)^3))`
Jak widać, w porównaniu z metodą Eulera w każdym wzorze jest dodatkowy człon, możliwy do skonstruowania dzięki temu, ze rozwiązujemy tu układ równań różniczkowych drugiego rzędu, w którym prawe strony nie zależą od prędkości.
Na koniec warto dodać, że obie te metody można oczywiście stosować do ruchu perturbowanego - prawe strony równań są wtedy bardziej skomplikowane ale nadal na podstawie schematów (6) oraz (11) i (12) można wyprowadzić odpowiednie wzory obliczeniowe.
Krok całkowania
Obie metody będziemy stosować ze stałym na całym przedziale całkowania, odpowiednio dobranym krokiem. W opisanym wyżej wariancie metody Leap Frog stały krok jest wręcz konieczny dla zachowania dokładności.