obrazek

Autor: mgr Tadeusz Sypek

obrazek

Numeryczne rozwiązanie problemu dwóch ciał

W niniejszej pracy, która ma charakter elementarny, propedeutyczny, ograniczymy się do problemu dwóch ciał. Zaniedbamy nieznaczne, zaburzające przyspieszenia grawitacyjne (zarówno dodatnie, jak i ujemne) działające na dwa ciała tworzące układ grawitacyjny, a wywołane działaniem trzeciej (lub kolejnej), odległej masy (np. oddziaływania perturbacyjne wywoływane przez Jowisza).

Jak wiadomo problem ten jest dokładnie rozwiązany analitycznie na gruncie newtonowskiej teorii grawitacji. Przy zadanych warunkach początkowych (energia potencjalna układu i moment pędu) możemy wyznaczyć elementy orbity: okres ruchu, półoś wielką, (w przypadku elipsy) mimośród. Rozwiązanie wspomnianego zagadnienia leży raczej poza możliwościami ucznia szkoły średniej: potrzebna jest znajomość całki energii lub wzoru Bineta, wymagane jest odpowiednie przygotowanie matematyczne (teoria krzywych stożkowych, współrzędne biegunowe). Korzystanie z gotowych wzorów – jak się wydaje – jest mało kształcące.

Rozsądną alternatywą dydaktyczną może być podejście numeryczne. Odpowiedni algorytm jest przedstawiony i dokładnie skomentowany np. w [3], [4]. Do jego zrozumienia całkowicie wystarczają wiadomości na poziomie niższych klas liceum (kinematyka i dynamika ruchów jednostajnie zmiennych, prawo powszechnego ciążenia). Przypomnijmy tutaj krótko i bardzo ogólnie, że krzywe stożkowe są przybliżane odcinkami obrazek (węzły główne algorytmu), a prędkości v_i+1 są liczone w momentach, gdy ciało znajduje się w węźle pomocniczym (i=1.2....,k+1), k – krok obliczeń). Przejdźmy do szczegółów.

Rozpatrzmy ciało centralne o dostatecznie dużej masie M (np. Słońce) oraz ciało o nieznacznej masie m (np. planeta). Oznacza to, że nie musimy uwzględniać ruchu ciała o masie M. Niech w chwili początkowej położenie planety jest określone jej promieniem wodzącym obrazek . Wiemy, że na planetę działa siła grawitacji skierowana ku Słońcu obrazek (p. Rys. 1a).

Rys. 1

obrazek

Szukamy składowych przyspieszenia obrazek (p . Rys. 1b).

Z II zasady dynamiki Newtona oraz prawa powszechnego ciążenia:

obrazek ,

obrazek (1)

W celu uproszczenia rachunków można przyjąć (jest to sprawa wyboru jednostek), GM ≡ 1; ostatecznie wzory (1) ulegają znacznemu uproszczeniu (2):

obrazek (2)

NUMERYCZNE OBLICZENIA BĘDĄ PRZEBIEGAĆ W NASTĘPUJĄCYCH ETAPACH:

wybieramy początkowe współrzędne położenie planety x(0) i y(0),
przyjmujemy początkowe wartości składowych wektora prędkości v_x(0) i v_y(0),
wybieramy krok rachunków k,

REKURENCJA

liczymy współrzędne planety x(k) i y(k) po pierwszym kroku (i) wykonywanego algorytmu,
wyznaczamy – dla wygody oznaczonego r (instrukcja podstawiania),
liczymy składowe wektora przyspieszenia i składowe wektora prędkości po pierwszym kroku wykonywanego programu,
zwiększamy wartość argumentu pętli i + 1, itd.

Listing 1

Tekst programu
Program

'NUMERYCZNE OBLICZANIE ORBIT	KOMENTARZ LISTINGU Język programu FreeBASIC 0.4.5
Cls Screen 9 Color 4, 45 Locate 12, 64: Print "x" Locate 2, 42: Print "y" Window(-320, 175)-(319, -174) Color 3, 63 Line (-250, 0)-(250, 0), 40 Line (0, 170)-(0, -170), 40	Współrzędne ekranowe Osie układu
'WPROWADZANIE DANYCH Dim as double x, y, r, vx, vy, ax, ay, k	Aby nie mnożyć deklarowanych zmiennych wartości początkowe współrzędnych położenia i prędkości przyjęto zamiast x_o_,y_o, v_xo, v_yo odpowiednio x, y, vx, vy i wykorzystano instrukcję podstawiania
'WARUNKI POCZĄTKOWE x = 0.5: y = 0 vx = 0: vy = 1.6	Współrzędne położenia początkowego Składowe prędkości początkowej
'RACHUNKI ZASADNICZE (REKURENCJA) k = 1 / 150	Krok rachunków
For i = 1 To 10000 x = x + vx * k: y = y + vy * k r = (x * x + y * y) ^ (-3 / 2) ax = -x * r: ay = -y * r vx = vx + ax * k: vy = vy + ay * k Pset(180 * x, 180 * y) Next Sleep	Wyprowadzenie krzywej na ekran

Dokładność algorytmu

Jak wiadomo obliczenia komputerowe nie zawsze gwarantują otrzymanie dokładnych wyników. Ogólnie błędy wyników metod numerycznych związane są m. in z:

dokładnością i ilością danych wejściowych; jeżeli dane te będą mało dokładne lub będzie ich mało, wynik może być obarczony błędem,
błędami popełnionymi przez komputer w trakcie konwersji liczb pomiędzy układami liczbowymi jak i reprezentacją tych liczb w pamięci komputera,
krokiem rachunków i wielkością kroku,
błędami związanymi z zaokrągleniami liczb.

Dla przykładu zilustrujmy, jak wielkość kroków oraz ich ilość wpływa na dokładność wykonywania programu (Rys. 2a-d). Poniższa tabela posłuży do analizy rysunków.

Dane wejściowe (takie same dla Rys. 2a-d)

x₀ =0,5, y₀ = 0
v_xo = 0, v_yo = 1,6

Nr rysunku	Krok rachunków k	Ilość kroków i
Rys. 2a		10000
Rys. 2b		100000
Rys. 2c		100000
Rys. 2d		1000000

Rys. 2

obrazek

Algorytm wydaje się być bezpieczny, ponieważ nie występuje w nim w żadnym miejscu dzielenie przez 0. Krok rachunków zawsze można zadeklarować jako liczbę, a nie iloraz. Dla k=0 algorytm nie pracuje. Natomiast zwielokrotniając okres obiegu (Rys. 2d – 250 razy) upewniamy się o stabilności i efektywności algorytmu dla odpowiednio wybranych k.

Dokładność postępowania rekurencyjnego zależy od wielkość kroku rachunków k. Rysunki 2a-d zdają się potwierdzać to logiczne twierdzenie … błąd jest proporcjonalny do kwadratu kroku k, jeżeli tysiąckrotnie zmniejszymy krok rachunków, to wynik będzie tysiąc milion razy dokładniejszy… [4].

Krzywe, które są trajektoriami w problemie dwóch ciał są krzywymi stożkowymi: okrąg, elipsa, parabola, hiperbola, a w szczególnym przypadku dwie proste – krzywa stożkowa niewłaściwa.

Na początek zajmijmy się elipsą i na jej przykładzie pokażemy jak kształt elipsy oraz jej położenie zależy od warunków początkowych (Rys. 3a-d, dalej program Listing 1).

Rys. 3

a)	b)
x = 0,5, y = 0 v_x = 0,5, v_y = 1,6	x = 0, y = 0 v_x = -0,5, v_y = 1,6
c)	d)
x = 0, y = 0,5 v_x = 1,6, v_y = 0	x = -0,5, y = -0,5 v_x = 1,2, v_y = -0,5

Rys. 4

obrazek

1-sza prędkość kosmiczna:

obrazek ,

obrazek

(tutaj vyo = 1,82),

obrazek

2-ga prędkość kosmiczna

obrazek

Zwiększając ilość kroków obliczeń i oraz wyprowadzając na ekran wartość promienia wodzącego obrazek i długość wektora prędkości obrazek , stwierdzamy ławo, że gdy r→∞,v→0.

Aby uwidocznić drugą gałąź paraboli (podobnie z hiperbolą – p. Rys. 4c, d), zamiast linii:

Pset(50 * x, 180 * y) piszemy:

Pset(50 * x, 180 * y): Pset(50 * x, -180 * y).

obrazek

(tutaj v_yo = 2,5)

Wszystkie powyższe rysunki możemy poglądowo zestawić na diagramie:

Rys. 5

W szkole średniej rozpatrujemy rzut poziomy, czyli ruch w jednorodnym polu grawitacyjnym z prędkością początkową w płaszczyźnie prostopadłej do kierunku pola. Torem ruchu jest parabola o wierzchołku w punkcie rzutu. Odpowiada on ruchowi ciała rzuconego poziomo w polu grawitacyjnym Ziemi, z pewnej wysokości, przy założeniu braku oporu ruchu i prędkości znacznie mniejszej od I prędkości kosmicznej. Wówczas pole grawitacyjne Ziemi można uznać w przybliżeniu za jednorodne.

Ale w rzeczywistości ruch odbywa się po elipsie o bardzo małym mimośrodzie, a środek Ziemi pokrywa się z jednym ognisk elipsy.

Elementy orbity

Pozostańmy przy tematyce astronomicznej. Dostosujmy Listing 1 do potrzeb dydaktycznych; przyjrzyjmy się dokładnie parametrom (zmiennym) występującym w programie. Prześledzimy jak zmieniają się współrzędne (x,y) planety w czasie, zbadamy wartości i znaki składowych prędkości [v_x,v_y], wyznaczymy odległość planety od Słońca (r), w tym odległość maksymalną (x max), mimośród orbity (e), półoś wielką (a), okres obiegu (T). Dlatego też Listing 2 jest bogato zaopatrzony w instrukcje wyprowadzające wartości zmiennych na ekran.

Rys. 6

obrazek

Listing 2

Tekst programu
Program

'NUMERYCZNE OBLICZANIE ORBIT
‘ELEMENTY ORBITY

Cls: Screen 9

Color 4, 45
Locate 12, 64: Print "x"
Locate 2, 42: Print "y"

Window(-320, 175)-(319, -174)

Color 3, 63
Line (-250, 0) - (250, 0), 40: Line (0, 170) - (0, -170), 40

'WPROWADZANIE DANYCH  

Dim as double x, y, r, a, vx, vy, v, ax, ay, k
Dim as double xmax, i, T, e

Color 3, 63


'WARUNKI POCZĄTKOWE

Locate 1, 2: Print "DANE WEJSCIOWE"

x = 0.5: y = 0

Locate 1, 18: Print "xo="; x
Locate 2, 18: Print "yo="; y

vx = 0: vy = 1.95

Locate 1, 28: Print "vxo="; vx
Locate 2, 28: Print "vyo="; vy

Color 5, 63

xmax = x


'RACHUNKI ZASADNICZE (REKURENCJA)

Color 44, 63

k = 1 / 10000

Locate 24, 2: Print "k="; k

For i = 1 To 62200
  Locate 24, 30: Print "i="; i
  
  Color 4, 63

  x = x + vx * k: y = y + vy * k

  Locate 17, 52: Print "x="; x
  Locate 18, 52: Print "y="; y

  If Abs(x) > Abs(xmax) Then xmax = Abs(x)

  r = (x * x + y * y) ^ (-3/2)
  
  Locate 19, 52: Print "r="; r ^ (-1 / 3)

  Color 5, 63

  ax = -x * r: ay = -y * r

  Locate 21, 52: Print "vx="; vx

  vx = vx + ax * k: vy = vy + ay * k

  Locate 22, 52: Print "vy="; vy

  v = Sqr(vx * vx + vy * vy)

  Locate 23, 52: Print "v="; v

  Color 1,63

  Pset(150 * x, 150 * y)

  Locate 2, 52: Print "xmax="; xmax

Next

Color 34, 63

a = 0.5 * (xmax + x)

Locate 4, 52: Print "a="; a

e = (xmax - a) / a

Locate 5, 52: Print "e="; e

T = k * i

Locate 6, 52: Print "T="; T

Sleep

W razie potrzeby dodatkowo możemy prześledzić składowe i wartość całkowitego przyspieszenia, kąt między wektorami przyspieszenia i promienia wodzącego, prędkość ucieczki, itp. Otwiera to możliwości różnych zastosowań.

Dokument ten rozpowszechniany jest zgodnie z zasadami licencji
GNU Free Documentation License.

Pytania proszę przesyłać na adres email: i-lo@eduinf.waw.pl

W artykułach serwisu są używane cookies. Jeśli nie chcesz ich otrzymywać,
zablokuj je w swojej przeglądarce.
Informacje dodatkowe