Informacje dla studentów - Stud - Badanie dokładności całkowania numerycznego

Autorem tego opracowania jest dr Piotr A. Dybczyński z Instytutu Obserwatorium Astronomiczne UAM w Poznaniu.

Najpierw prezent

Zanim przejdziemy do właściwego zagadnienia musimy stworzyć podstawowy warsztat rachunkowy dla ruchu keplerowskiego. Skorzystamy tu z formuł dla heliocentrycznego ruchu ciała o zaniedbywalnej masie (np. komety czy planetoidy). Oznacza to, że jednostką odległości będzie jednostka astronomiczna (AU), jednostką czasu doba a jednostką masy masa Słońca. W tych jednostkach stała grawitacji jest równa kwadratowi tzw. stałej Gaussa: G = k² = (0.01720209895)² . Zagubionym w obszarze mechaniki nieba polecam wykład Profesora Sławomira Breitera.

Musimy umieć wyliczyć pozycję ciała na zadany moment czasu z danych elementów orbity. Najtrudniejszym miejscem jest tu rozwiązanie równania Keplera. Sposobów praktycznego rozwiązania tego zagadnienia jest wiele, ja zaproponuję użycie jednej procedury, która wylicza anomalię prawdziwą po podaniu mimośrodu, odległości peryhelium i interwału czasu w dniach (może być ujemny) od momentu przejścia praz peryhelium. Pomysł i algorytm tej procedury pochodzi z pracy Profesora Grzegorza Sitarskiego (1968, Acta Astronomica,vol.18,no.2,strony 197-205).

Poniżej prezent, czyli gotowa i przetestowana procedura realizująca ten rachunek:

/****************************************************************************
******** ANOMALIAT.C *********
*****************************************************************************
******** **********
******** **********
*****************************************************************************
** **
** Pomysł i algorytm: **
** G.Sitarski (1968), Acta Astronomica,vol.18,no.2,pp 197-205, **
** **
*****************************************************************************
** Programmed by: Piotr A. Dybczyński, Last modified: 3.03.2012 **
****************************************************************************/

#define GAUSS 0.01720209895
#define E_SMALL (1.e-12)

#include<math.h>

double anomaliat(double dt,double q,double e)
{
double x,M,M1,E,E1,w,F,v,u,F1,z,n,mi;

x = (1.-e)/(1.+e);

mi=GAUSS*GAUSS; /* tylko dla ciał o zaniedbywalnej masie */

if(e<(1.-E_SMALL)) /* elipsa */
{
n = sqrt(mi*(1.-e)*(1.-e)*(1.-e)/(q*q*q));
M = n*dt;
M = fmod(M,2.0*M_PI);

M1 = fabs(M);
E1 = M_PI;
do {
E = E1;
E1 = E + (v=M1-E+e*sin(E))/(1.-e*cos(E));
} while( fabs(v) > fabs(M1-E1+e*sin(E1) ));

if(1.!=1.+cos(E1/2.)) w = sqrt(1./x)*tan(E1/2.);
else w = 1.0E+70; /* theta = 180 deg*/
}

else if(e>(1.+E_SMALL)) /*hiperbola */
{
M = sqrt(mi*(e-1.)*(e-1.)*(e-1.)/(q*q*q))*dt;

M1 = fabs(M);
F1 = 1;
while( (z = M1+F1-e*sinh(F1)) > 0.0 ) F1+=1.;
do {
F = F1;
F1 = F + (v=z)/(e*cosh(F)-1.);
} while( fabs(v) > fabs(z=M1+F1-e*sinh(F1)));
w = sqrt(-1./x)*tanh(F1/2.);
}
else /* bliskoparabola: e = 1 +/- E_SMALL */
{
M = sqrt(mi/(2.*q*q*q))*dt;

u = fabs(M)*1.5;
v = sqrt(1.+u*u);
w = pow((v+u),1./3.) - pow((v-u),1./3.);
}
return( M<0? -2.*atan(w): 2.*atan(w) );
}

Pobierz

Kilka słów wyjaśnienia: procedura definiuje stałą E_SMALL i jeśli mimośród różni się od jedynki o mniej niż E_SMALL korzysta ze wzorów dla paraboli, w pozostałych przypadkach rozwiązywane jest eliptyczne lub hiperboliczne równanie Keplera. Daje to optymalną dokładność numeryczną.

Pozostaje napisać prosty program testowy, który poprosi o podanie q, e oraz dt i wylicza anomalię prawdziwą f. W cytowanej pracy Profesora Sitarskiego jest tabelka z danymi testowymi ale prościej jest chyba porównać działanie takiego programu z moim, zamieszczonym na pracowni w katalogu /home/COMMON/varia programem anomalia.

Położenie i prędkość ciała na orbicie keplerowskiej

Orbitę keplerowską opisuje się przy pomocy tzw. elementów orbity. Zależnie od zagadnienia czy typu orbity stosuje się różne, równoważne zestawy sześciu elementów. Dla komet najwygodniejszym zestawem jest:

odległość peryhelium w AU q
mimośród e
moment przejścia peryhelium (data lub dzień juliański) T
długość węzła wstępującego Ω
argument peryhelium ω
nachylenie (najczęściej do płaszczyzny ekliptyki) i

W innych zagadnieniach (ale i w ruchu komet) używane bywają również inne elementy:

parametr orbity p = q(1 + e)
wielka półoś orbity a, dla elipsy: a = q / (1 - e), dla hiperboli: a = q / (e - 1)
odległość aphelium (dla elipsy) Q = q (1 + e) / (1 - e)
ruch średni n = k / a³
anomalia średnia M_o na jakiś ustalony moment czasu

Położenie ciała na orbicie keplerowskiej opisujemy tradycyjnie przy pomocy jednej z anomalii: prawdziwej f, mimośrodowej E lub średniej M. Ponieważ będziemy korzystać z załączonej wyżej procedury użyjemy anomalii prawdziwej.

W płaskim (tzw. orbitalnym) układzie odniesienia Oξη związanym z płaszczyzną orbity, ze środkiem w masie centralnej (Słońcu), z osią Oξ skierowaną do peryhelium a osią Oη prostopadle do niej i zgodnie z kierunkiem ruchu ciała, składowe położenia i prędkości wyrażają się prostymi wzorami:

`r = p/(1+e*cos(f))`

`xi = r*cos(f)`

`eta = r*sin(f)`

`dotxi = -sqrt(mu/p)*sin(f)`

`doteta = sqrt(mu/p)*(cos(f)+e)`

Aby uzyskać prostokątne składowe położenia i prędkości w układzie równikowym (w takim będziemy pracować) przechodzimy po drodze przez układ ekliptyczny (względem którego mamy określone elementy) a cała transformacja wymaga w sumie czterech obrotów (trzecia składowa położenia i prędkości w układzie orbitalnym jest równa zero):

`r_(alphadelta) = R_x(-epsilon)R_z(-Omega)R_x(-i)R_z(-omega)r_(xieta)`

gdzie `epsilon` jest nachyleniem ekliptyki do równika niebieskiego. Proszę używać wartości dla epoki J2000: `epsilon` = 23° 26' 21.406" .

Teraz zadanie

Rys. 1 Zależność dokładności i czasu obliczeń od
użytego kroku całkowania.
Krok jest podany jako ułamek okresu.
rt - wektor położenia wyliczony ze wzorów
ruchu keplerowkiego, r - wektor położenia
z całkowania numerycznego.
Autor: mgr F.Berski.

Napisać program, który prosi o podanie elementów dowolnej (niezdegenerowanej) orbity keplerowskiej oraz dwa dowolne momenty czasu `t_1` i `t_2`. Następnie wylicza z powyższych wzorów położenia i prędkości na moment pierwszy (jako dane startowe) i przeprowadza badanie dokładności i czasochłonności całkowania numerycznego zagadnienia keplerowskiego dwoma metodami na `(t_1,t_2)`, w zależności od długości kroku całkowania.

Ale jak całkować numerycznie?!

Opis dwóch prostych metod całkowania znajdziecie Państwo w tym materiale.

Program realizujący zadanie musi więc zapytać użytkownika którą metodę ma badać a następnie wykonać serię kilkudziesięciu (może więcej) całkowań numerycznych i przedstawić wynik badania w postaci dwóch nałożonych na siebie wykresów. Na osi poziomej ma być krok całkowania a na pionowej błąd całkowania w chwili `t_2` oraz czas obliczeń. Zakres zmian kroku całkowania należy dobierać do danych podanych przez użytkownika, w szczególności do wielkiej półosi całkowanej orbity i jej mimośrodu.

Czas obliczeń proszę mierzyć funkcją clock() a dokładność wyniku całkowania długością wektora opisującego różnicę przestrzennych położeń naszego ciała w chwili `t_2` uzyskanych z całkowania numerycznego i z dokładnych wzorów podanych wyżej.

Końcowy wykres będący wynikiem działania programu może wyglądać np. jak na Rys.1 (to wynik badania innej, dokładniejszej metody).

≡


Strona Wydziałowa Informacje ogólne Materiały dla studentów, z którymi mam zajęcia. Co nowego? ≡ Powered by	Autorem tego opracowania jest dr Piotr A. Dybczyński z Instytutu Obserwatorium Astronomiczne UAM w Poznaniu. Najpierw prezent Zanim przejdziemy do właściwego zagadnienia musimy stworzyć podstawowy warsztat rachunkowy dla ruchu keplerowskiego. Skorzystamy tu z formuł dla heliocentrycznego ruchu ciała o zaniedbywalnej masie (np. komety czy planetoidy). Oznacza to, że jednostką odległości będzie jednostka astronomiczna (AU), jednostką czasu doba a jednostką masy masa Słońca. W tych jednostkach stała grawitacji jest równa kwadratowi tzw. stałej Gaussa: G = k² = (0.01720209895)² . Zagubionym w obszarze mechaniki nieba polecam wykład Profesora Sławomira Breitera. Musimy umieć wyliczyć pozycję ciała na zadany moment czasu z danych elementów orbity. Najtrudniejszym miejscem jest tu rozwiązanie równania Keplera. Sposobów praktycznego rozwiązania tego zagadnienia jest wiele, ja zaproponuję użycie jednej procedury, która wylicza anomalię prawdziwą po podaniu mimośrodu, odległości peryhelium i interwału czasu w dniach (może być ujemny) od momentu przejścia praz peryhelium. Pomysł i algorytm tej procedury pochodzi z pracy Profesora Grzegorza Sitarskiego (1968, Acta Astronomica,vol.18,no.2,strony 197-205). Poniżej prezent, czyli gotowa i przetestowana procedura realizująca ten rachunek: /************************************************************************** **** ANOMALIAT.C ***** ************************************************************************* **** ****** **** ****** ************************************************************************* Pomysł i algorytm: G.Sitarski (1968), Acta Astronomica,vol.18,no.2,pp 197-205, ************************************************************************* Programmed by: Piotr A. Dybczyński, Last modified: 3.03.2012 *************************************************************************/ #define GAUSS 0.01720209895 #define E_SMALL (1.e-12) #include<math.h> double anomaliat(double dt,double q,double e) { double x,M,M1,E,E1,w,F,v,u,F1,z,n,mi; x = (1.-e)/(1.+e); mi=GAUSSGAUSS; /* tylko dla ciał o zaniedbywalnej masie / if(e<(1.-E_SMALL)) / elipsa / { n = sqrt(mi(1.-e)(1.-e)(1.-e)/(qqq)); M = ndt; M = fmod(M,2.0M_PI); M1 = fabs(M); E1 = M_PI; do { E = E1; E1 = E + (v=M1-E+esin(E))/(1.-ecos(E)); } while( fabs(v) > fabs(M1-E1+esin(E1) )); if(1.!=1.+cos(E1/2.)) w = sqrt(1./x)tan(E1/2.); else w = 1.0E+70; /* theta = 180 deg/ } else if(e>(1.+E_SMALL)) /hiperbola / { M = sqrt(mi(e-1.)(e-1.)(e-1.)/(qqq))dt; M1 = fabs(M); F1 = 1; while( (z = M1+F1-esinh(F1)) > 0.0 ) F1+=1.; do { F = F1; F1 = F + (v=z)/(ecosh(F)-1.); } while( fabs(v) > fabs(z=M1+F1-esinh(F1))); w = sqrt(-1./x)tanh(F1/2.); } else / bliskoparabola: e = 1 +/- E_SMALL / { M = sqrt(mi/(2.qqq))dt; u = fabs(M)1.5; v = sqrt(1.+uu); w = pow((v+u),1./3.) - pow((v-u),1./3.); } return( M<0? -2.atan(w): 2.atan(w) ); } Pobierz Kilka słów wyjaśnienia: procedura definiuje stałą E_SMALL* i jeśli mimośród różni się od jedynki o mniej niż E_SMALL korzysta ze wzorów dla paraboli, w pozostałych przypadkach rozwiązywane jest eliptyczne lub hiperboliczne równanie Keplera. Daje to optymalną dokładność numeryczną. Pozostaje napisać prosty program testowy, który poprosi o podanie q, e oraz dt i wylicza anomalię prawdziwą f. W cytowanej pracy Profesora Sitarskiego jest tabelka z danymi testowymi ale prościej jest chyba porównać działanie takiego programu z moim, zamieszczonym na pracowni w katalogu /home/COMMON/varia programem anomalia. Położenie i prędkość ciała na orbicie keplerowskiej Orbitę keplerowską opisuje się przy pomocy tzw. elementów orbity. Zależnie od zagadnienia czy typu orbity stosuje się różne, równoważne zestawy sześciu elementów. Dla komet najwygodniejszym zestawem jest: odległość peryhelium w AU q mimośród e moment przejścia peryhelium (data lub dzień juliański) T długość węzła wstępującego Ω argument peryhelium ω nachylenie (najczęściej do płaszczyzny ekliptyki) i W innych zagadnieniach (ale i w ruchu komet) używane bywają również inne elementy: parametr orbity p = q(1 + e) wielka półoś orbity a, dla elipsy: a = q / (1 - e), dla hiperboli: a = q / (e - 1) odległość aphelium (dla elipsy) Q = q (1 + e) / (1 - e) ruch średni n = k / a³ anomalia średnia *M_o* na jakiś ustalony moment czasu Położenie ciała na orbicie keplerowskiej opisujemy tradycyjnie przy pomocy jednej z anomalii: prawdziwej f, mimośrodowej E lub średniej M. Ponieważ będziemy korzystać z załączonej wyżej procedury użyjemy anomalii prawdziwej. W płaskim (tzw. orbitalnym) układzie odniesienia *Oξη* związanym z płaszczyzną orbity, ze środkiem w masie centralnej (Słońcu), z osią Oξ skierowaną do peryhelium a osią Oη prostopadle do niej i zgodnie z kierunkiem ruchu ciała, składowe położenia i prędkości wyrażają się prostymi wzorami: `r = p/(1+ecos(f))` `xi = rcos(f)` `eta = rsin(f)` `dotxi = -sqrt(mu/p)sin(f)` `doteta = sqrt(mu/p)(cos(f)+e)` Aby uzyskać prostokątne składowe położenia i prędkości w układzie równikowym (w takim będziemy pracować) przechodzimy po drodze przez układ ekliptyczny (względem którego mamy określone elementy) a cała transformacja wymaga w sumie czterech obrotów (trzecia składowa położenia i prędkości w układzie orbitalnym jest równa zero): `r_(alphadelta) = R_x(-epsilon)R_z(-Omega)R_x(-i)R_z(-omega)r_(xieta)` gdzie `epsilon` jest nachyleniem ekliptyki do równika niebieskiego. Proszę używać wartości dla epoki J2000: `epsilon` = 23° 26' 21.406"* . Teraz zadanie Rys. 1 Zależność dokładności i czasu obliczeń od użytego kroku całkowania. Krok jest podany jako ułamek okresu. rt - wektor położenia wyliczony ze wzorów ruchu keplerowkiego, r - wektor położenia z całkowania numerycznego. Autor: mgr F.Berski. Napisać program, który prosi o podanie elementów dowolnej (niezdegenerowanej) orbity keplerowskiej oraz dwa dowolne momenty czasu `t_1` i `t_2`. Następnie wylicza z powyższych wzorów położenia i prędkości na moment pierwszy (jako dane startowe) i przeprowadza badanie dokładności i czasochłonności całkowania numerycznego zagadnienia keplerowskiego dwoma metodami na `(t_1,t_2)`, w zależności od długości kroku całkowania. Ale jak całkować numerycznie?! Opis dwóch prostych metod całkowania znajdziecie Państwo w tym materiale. Program realizujący zadanie musi więc zapytać użytkownika którą metodę ma badać a następnie wykonać serię kilkudziesięciu (może więcej) całkowań numerycznych i przedstawić wynik badania w postaci dwóch nałożonych na siebie wykresów. Na osi poziomej ma być krok całkowania a na pionowej błąd całkowania w chwili `t_2` oraz czas obliczeń. Zakres zmian kroku całkowania należy dobierać do danych podanych przez użytkownika, w szczególności do wielkiej półosi całkowanej orbity i jej mimośrodu. Czas obliczeń proszę mierzyć funkcją clock() a dokładność wyniku całkowania długością wektora opisującego różnicę przestrzennych położeń naszego ciała w chwili `t_2` uzyskanych z całkowania numerycznego i z dokładnych wzorów podanych wyżej. Końcowy wykres będący wynikiem działania programu może wyglądać np. jak na Rys.1 (to wynik badania innej, dokładniejszej metody). ≡