Symulacja ruchu w polu grawitacyjnym

Ciao,
ostatnio próbowałem napisać prosty program który symulował by ruch planety po orbicie danej gwiazdy. Kiedy jednak zaadaptowałem wzór na przyspieszenie w polu grawitacyjnym, moja planeta zaczęła krążyć nie po krzywej stożkowej ale po spirali (sic).

Sam wzór sprawdzałem kilkanaście razy, problem stanowi najprawdopodobniej błąd przy manipulowaniu typem double - kilkakrotnie odejmuję liczby, potem dzielę przez pierwiastek z tej różnicy - innymi słowy, koprocesor więcej domyśla się niż liczy.

Czy ktoś wie, jak można by usprawnić algorytm żeby zminimalizować błąd? Nie chodzi mi o uzyskanie perfekcyjnie eliptycznej orbity, ale na razie dewiacja jest zbyt duża żeby użyć programu do czegokolwiek sensownego.

Tu zamieszczam kod źródłowy w C/Win32 API: http://taurendil.w.interia.pl/graw.zip

Z góry dzięki za odpowiedź. </url>

Akurat zacząłem się niedawno bawić w podobną rzecz, tyle że u mnie wokół Ziemi krąży satelita... :)

Jedna rzecz, którą widzę w Twoim kodzie to duża „odległość” od wyliczenia przyspieszenia do wpłynięcia jego na prędkość. Dokładnie mówiąc robisz tak: (1) dodajesz do pozycji planety jej prędkość (2) zwiększasz prędkość o przyspieszenie obliczone w poprzednim kroku (3) obliczasz przyspieszenie dla przesuniętej już pozycji.

Ja w pierwszym przybliżeniu miałem całe obliczenia zrobione dla kroku jednocześnie: obliczałem nową pozycję dodając (starą) prędkość do starej pozycji oraz obliczałem nową prędkość dodając przyspieszenie obliczone dla starej pozycji do starej prędkości. Niestety to rozwiązanie tez było niedokładne, całkowita energia satelity minimalnie, ale nieustannie rosła mi z czasem (a więc powolutku oddalał on się od Ziemi po spirali). Pomogło dopiero zaimplementowanie metody Cauchy'ego-Eulera:

static void NextStep()
{
	static const double h = 0.15;

	Vector f1 = Speed;
	Vector f2 = A(Position);

	Vector _y1 = Position + f1*h;
	Vector _y2 = Speed    + f2*h;

	Vector _f1 = _y2;
	Vector _f2 = A(_y1);

	Position += (f1 + _f1) * (h*0.5);
	Speed    += (f2 + _f2) * (h*0.5);

}

(Zdefiniowałem sobie klasę Vector z przeładowanymi operatorami dla wygody obliczeń, funkcja A zwraca wektor będący przyśpieszeniem satelity w danej pozycji)

Metoda C-E jest fajna, ale poczytaj jeszcze o metodzie Rungego-Kutty czterokrokowej (czwartego rzędu). Jest najprecyzyjniejszą w tego typu obliczeniach:
http://pl.wikipedia.org/wiki/Metoda_Rungego-Kutty

Dodam jeszcze, że dla dwóch ciał można rozwiązać równanie ruchu całkowicie analitycznie i odnaleźć wszystkie zależności. Tym samym można spokojnie wyeliminować błąd numeryczny za pomocą odpowiedniej jego kompensacji. Jeżeli masz trzy i więcej ciał nie można już w ten prosty sposób "naprawić" niedokładności komputera ponieważ nie istnieje rozwiązanie analityczne. Jedyną metodą jest przybliżone rozwiązanie numeryczne, dla którego w metodzie R-K zaczyna się "grzanie numeryczne" w okolicach 10-12 obrotu :)

Dziękuję za obydwie sugestie. Problem niedokładności zniknął (zmniejszył się do akceptowalnych rozmiarów) przy zastosowaniu metody Cauchy'ego-Eulera. Chwilowo nie mogę sprawdzić jak poprawi się precyzja przy metodzie R-K, jako że wiedza matematyczna na poziomie licealnym nie pozwala na obliczanie różniczki funkcji dwóch zmiennych ;-)

Koziołek napisał(a)

Metoda C-E jest fajna, ale poczytaj jeszcze o metodzie Rungego-Kutty czterokrokowej (czwartego rzędu). Jest najprecyzyjniejszą w tego typu obliczeniach:
http://pl.wikipedia.org/wiki/Metoda_Rungego-Kutty

Dodam jeszcze, że dla dwóch ciał można rozwiązać równanie ruchu całkowicie analitycznie i odnaleźć wszystkie zależności. Tym samym można spokojnie wyeliminować błąd numeryczny za pomocą odpowiedniej jego kompensacji. Jeżeli masz trzy i więcej ciał nie można już w ten prosty sposób "naprawić" niedokładności komputera ponieważ nie istnieje rozwiązanie analityczne. Jedyną metodą jest przybliżone rozwiązanie numeryczne, dla którego w metodzie R-K zaczyna się "grzanie numeryczne" w okolicach 10-12 obrotu :)

Do problemów fizycznych są inne specjalne metody numeryczne dla RR, które zachowują energię układu i wtedy rozwiązanie nie rozjeżdża się.

Ale te metody mają inną wadę: nie nadają się do analizy perpetuum mobile (użytecznego, czyli II-go rodzaju), bo tu energia z definicji nie jest zachowana.
Pewnie dlatego naukowcy (tacy marni, czyli zawsze wszystkowiedzący) nie potrafią odtworzyć starych wynalazków typu perpetuum mobile.
Takie układy niekiedy można obliczyć poprawnie, ale z równań ruchu Newtona, czyli z sił, a nie z energii.

Np. koło grawitacyjne Besslera. Testowane wielokrotnie oficjalnie 300 lat temu i ostatecznie zatwierdzone urzędowo jako PM.
Zaraz potem Leibniz wymyślił zasadę zachowania energii.
On znał osobiście Besslera i sam potwierdzał wielokrotnie użyteczność tej maszyny.

W grawitacji energia nie jest zachowana w ogólnym przypadku.
Samotna planeta na orbicie kołowej i o dużym promieniu zachowuje energię, w innych przypadkach jest inaczej.

Na tej stronie jest namiar na program i plik z działającą przykładową symulacja silnika PM: http://free.of.pl/m/magin

MAGIN

To co zastosowałeś to jest algorytm Eulera i on ma taką wadę, że popełnia systematyczne błędy (stąd spirala).
Inne proste algorytmy: które nie mają tej wady:
algorytm Verleta
algorytm LeapFrog