Analiza numeryczna jest używana do obliczania przybliżeń rzeczy: wartości funkcji w pewnym punkcie, w którym znajduje się pierwiastek równania, lub rozwiązań układu równań różniczkowych. Jest to rozległy i ważny temat, ponieważ w praktyce większość rzeczywistych problemów w matematyce, nauce i technologii nie będzie miała jednoznacznego rozwiązania w formie zamkniętej (a nawet jeśli tak, obliczenie przy użyciu nieskończona precyzja - w końcu komputery reprezentują liczby ze skończoną precyzją). Ogólnie rzecz biorąc, istnieją kompromisy między dokładnością a szybkością obliczeniową.

W przypadku problemu trzech ciał mamy trzy masy punktowe w lokalizacjach początkowych $ \ mathbf {x} _i (0) $ z prędkościami $ \ mathbf {v} _i (0) $ , które chcemy obliczyć na późniejsze czasy $ t $ . Matematycznie chcemy znaleźć rozwiązanie dla systemu $$ \ mathbf {x} '_ i (t) = \ mathbf {v} _i (t), $$ $$ \ mathbf {v} '_ i (t) = \ mathbf {f} _i (t) / m_i, $$ $$ \ mathbf {f} _i (t) = Gm_i \ sum_ {j \ neq i} \ frac {m_j (\ mathbf {x} _j- \ mathbf {x} _i)} {|| \ mathbf {x} _i - \ mathbf {x} _j || ^ 3}. $$

Oczywistą metodą jest myślenie „jeśli zrobimy mały krok naprzód $ h $ w czasie, możemy przybliżyć wszystko tak, aby było liniowe”, więc tworzymy formuła, w której obliczamy stan w czasie $ t + h $ ze stanu w czasie $ t $ ( i tak dalej dla $ t + 2h $ i dalej): $$ \ mathbf {x} _i (t + h) = \ mathbf {x} _i (t) + h \ mathbf {v} _i (t), $$ $$ \ mathbf {v} _i (t + h ) = \ mathbf {v} _i (t) + h \ mathbf {f} _i (t). $$ To się nazywa metoda Eulera. Jest prosty, ale bywa niedokładny; błąd na krok to $ \ approx O (h ^ 2) $ i mają tendencję do narastania. Jeśli spróbujesz rozwiązać problem z dwoma ciałami, okrążające masy będą wykonywać precesyjną orbitę rozety z powodu nagromadzenia się błędów, zwłaszcza gdy zbliżają się do siebie.

Istnieje menażeria metod numerycznego rozwiązywania równań różniczkowych. Można użyć metod wyższego rzędu, które próbkują funkcje w większej liczbie punktów, a tym samym lepiej je aproksymują. Istnieją metody niejawne, które zamiast próbować znaleźć stan w późniejszym czasie tylko w oparciu o stan bieżący, szukają samoistnego stanu późnego i pośredniego. Najpoważniejsze metody rozwiązywania problemów z ODE zmniejszają także wielkość kroku $ h $ , gdy siły rosną podczas bliskich starć, aby zapewnić akceptowalną dokładność. Jak powiedziałem, to duży temat.

Jednak w przypadku symulacji mechanicznych możesz w szczególności przyjrzeć się metodom zaprojektowanym w celu zachowania energii i innych zachowanych ilości ( metody symplektyczne - są to metody używane przez profesjonalistów do długoterminowych obliczeń orbity) . Być może najprostsza jest częściowo ukryta metoda Eulera. Istnieje również metoda Verleta i skokowa integracja. Podoba mi się półjawna metoda Eulera, ponieważ jest super prosta (ale będąc metodą pierwszego rzędu, nadal nie jest zbyt dokładna): $$ \ mathbf {v} _i (t + h) = \ mathbf {v} _i (t) + h \ mathbf {f} _i (t), $$ $$ \ mathbf {x} _i ( t + h) = \ mathbf {x} _i (t) + h \ mathbf {v} _i (t + h). $$ Czy widzisz różnicę? Najpierw obliczasz zaktualizowaną prędkość, a następnie używasz jej do aktualizacji pozycji - drobna sztuczka, ale nagle orbity 2-ciałowe zachowują się dobrze.

Problem trzech ciał jest chaotyczny w prawdziwym matematycznym sensie. Wiemy, że istnieją sytuacje, w których drobne różnice w warunkach początkowych są skalowane w górę do dowolnie dużych różnic w późniejszych pozycjach (nawet jeśli wykluczymy bardzo bliskie podania między masami). Tak więc nawet przy arbitralnie dużej precyzji liczbowej nastąpi moment, w którym nasze obliczone orbity będą całkowicie błędne. Ogólny „styl” trajektorii może nadal być poprawny, dlatego też można bawić się częściowo ukrytym Eulerem, o ile nie planuje się żadnej misji kosmicznej na podstawie wyników.

„Rozwiązanie” problemu trzech ciał można zapisać jako parę równań różniczkowych, \ begin {align} \ vec {v} & = \ frac {\ mathrm d \ vec {x}} {\ mathrm dt} \\ \ vec {a} & = \ frac {\ mathrm d \ vec {v}} {\ mathrm dt} \ end {align} gdzie to drugie jest zwykle zapisywane w kategoriach siły, $ m \ vec {a} = \ vec {F} $ . Następnie używając definicji pochodnej, $$ \ frac {\ mathrm df} {\ mathrm dx} = \ lim_ {h \ to0} \ frac {f (x + h) -f (x)} {h}, $$ równania różniczkowe, od których zaczęliśmy, można zapisać jako równania różnicy skończonej , \ begin {align} \ vec x (t + \ mathrm dt) & \ simeq \ vec x (t) + \ vec v \, \ mathrm dt \\ \ vec v (t + \ mathrm dt) & \ simeq \ vec v (t) + \ vec a \, \ mathrm dt \ end {align} co jest wykonywane przy założeniu, że $ \ mathrm dt $ jest „tylko mały”, a nie nieskończenie mały.

To tych równań używa komputer do "rozwiązania" problemu trzech ciał: mając warunek początkowy dla każdego z ciał, komputer iteracyjnie przesuwa się naprzód w czasie, używając powyższych równań różnic skończonych gdzie siła jest $ n $ -body force.

Jak wspomniano w komentarzach, uproszczony model powyżej nosi nazwę integracja Eulera , która wcale nie jest dobrze dopasowana do tego problemu, ponieważ nie oszczędza energii. Lepszą opcją jest prędkość Verlet , która jest omawiana w innych postach na temat fizyki.SE (polecam ten mój post, ponieważ podaje przyzwoite, ale zwięzłe szczegóły implementacji).

Rozsądnie dostępny artykuł na temat problemu trzech ciał, obejmujący dyskusję na temat algorytmów numerycznych, to Musielak & Quarles (2014).Szersze omówienie bardziej przydatnych technik wyższego rzędu, zgodnie z sugestią homocomputeris, patrz Farrés et al (2012)

Kilka lat temu stworzyłem model wewnętrznego układu słonecznego, słońca, księżyca Ziemi i planet od Merkurego do Marsa.

Użyłem wariantu Leapfrog Rk4 i uzyskałem przyzwoite wyniki, okresy orbitalne zgodne z początkowymi konfiguracjami, które pozostawały stabilne w czasie. Nawet niektóre złe rzeczy, które otrzymałem, sugerowały dobre praktyki liczbowe.

Inicjowałbym moje planety z początkowym pędem, ze słońcem mającym zerowy początkowy pęd. Pęd planet nie wyłączył się, więc nastąpił ten dziwny dryf sieci w mniej więcej prostym kierunku, którego zrozumienie zajęło trochę czasu.

Oprócz uporządkowania tego, jedną z ważnych funkcji, które uznałem za ważną, była kolejność wprowadzania zmian w pozycji i szybkości.

Nie chcesz iterować po swoich ciałach, obliczać siły wypadkowej w jednym przejściu ciało po ciele, a następnie stosować swoje dynamiczne zmiany ciało po ciele. W przeciwnym razie, powiedzmy, zaktualizuje pozycję Wenus, zachowując niezmienioną pozycję Marsa, podczas gdy Ty aktualizujesz pozycję Ziemi. Twoje kroki czasu tracą synchronizację.

Przechowuj gdzieś w pamięci nowe przemieszczenia i prędkości podczas iteracji po ciałach i dopiero po znalezieniu swojej dynamiki w obecnym tiku, przeprowadź i zastosuj swoje zmiany. Przepraszamy i gratulacje, jeśli złapałeś to za pierwszym razem. Trochę mnie to zaskoczyło.

Jednym ze sposobów symulacji trzech ciał jest użycie przyrostów czasu, w których udaje się, że każde ciało przyspiesza w stałym tempie przez krótki okres czasu. Zaczynasz od wektorów pozycji początkowej, wektorów prędkości i mas ciał, a następnie używasz ich do obliczenia wektora przyspieszenia każdego ciała. W fizyce newtonowskiej wektor przyspieszenia grawitacyjnego ciała A z ciała B jest określony przez równanie $$ \ vec {a_ {GA}} = - \ frac {G \ left (\ vec { A} - \ vec {B} \ right) M_B} {|| \ vec {A} - \ vec {B} || ^ 3} $$ z $ G $ to stała grawitacyjna, $ \ vec {A} $ to wektor położenia ciała A, $ \ vec {B} $ to wektor położenia ciała B, $ M_B $ to masa ciała B, a $ \ vec {a_ {GA}} $ to wektor przyspieszenia grawitacyjnego ciała A. Aby znaleźć całkowity wektor przyspieszenia grawitacyjnego ciała A, należy dodać każdy z wektorów przyspieszenia grawitacyjnego ciała A z każdego z pozostałych organów.

Wektory położenia i wektory prędkości można znaleźć na końcu kroku czasowego $ n $ , korzystając z ogólnych formuł $ $ \ vec {x_n} = \ frac {1} {2} \ vec {a_ {n-1}} s ^ 2 + \ vec {v_ {n-1}} s + \ vec {x_ {n-1}} $$ $$ \ vec {v_n} = \ vec {a_ {n-1}} s + \ vec {v_ {n-1}} $$ span > gdzie $ \ vec {a_ {n-1}} $ jest wektorem przyspieszenia dla kroku czasowego $ n-1 $ , $ \ vec {v_ {n-1}} $ to wektor prędkości dla kroku czasowego $ n -1 $ , $ \ vec {x_ {n-1}} $ to wektor pozycji $ n -1 $ , $ s $ to długość przyrostów kroku czasowego, $ \ vec {x_n} $ jest wektorem pozycji dla kroku czasowego $ n $ i $ \ vec {v_n} $ Zakres > będący wektorem prędkości dla kroku czasowego $ n $ . Na każdym kroku czasowym przeliczasz również wektory przyspieszenia, ponieważ nie są one stałe.

Ta metoda symulowania ciał $ n $ działa dla dowolnego prawa sił, co oznacza, że ​​możesz jej używać do symulowania $ n $ treści dla prawa siłowego $ F = f (r) $ , z $ f (r) $ to pewna funkcja odległości, chociaż w przypadku prawa siły $ F = f (r) $ równanie, które określiłem wcześniej dla obliczenia przyspieszenia treści A z ciała B staje się $$ \ vec {a_ {GA}} = - \ frac {\ vec {A} - \ vec {B}} {|| \ vec { A} - \ vec {B} ||} GM_Bf \ left (|| \ vec {A} - \ vec {B} || \ right) $$ ale inne formuły, o których wspomniałem, w tym część o dodawanie wektorów przyspieszenia pozostaje bez zmian.

Okazuje się, że znając trajektorie parami można wyznaczyć równanie ruchu ciał. Równania ruchu są zapisywane jako $$ \ frac {d ^ 2 \ vec r ^ k} {ds ^ 2} = - G \ sum_ {n = 1, n \ ne k} ^ N \ frac {m_n (\ vec r ^ k- \ vec r ^ n)} {| \ vec r ^ k- \ vec r ^ n | ^ 3} \ tag {1} $$ Rozwiązujemy pomocniczy problem parowania ciał o zredukowanej masie bezwładnościowej $$ \ frac {m_n m_k} {\ sum_ {k = 1} ^ N m_k} \ frac {d ^ 2 \ vec R ^ {kn}} {ds ^ 2} = -G \ frac {m_n m_k \ vec R ^ {kn}} {| \ vec R ^ {kn} | ^ 3} \ tag {2} $$ Równania te różnią się w różnych warunkach początkowych. Odejmij od równania (2) równanie (1), otrzymujemy $$ \ frac {d ^ 2 R_0 ^ k} {ds ^ 2} - \ frac {d ^ 2 r ^ k} {ds ^ 2} = G \ sum_ {n = 1, n \ ne k} ^ N \ frac {m_n (\ vec r ^ k- \ vec r ^ n)} {| \ vec r ^ k- \ vec r ^ n | ^ 3} - G \ sum_ {n = 1, n \ ne k} ^ N \ frac {m_n \ vec R ^ {kn}} {| \ vec R ^ {kn} | ^ 3} = 0, \ vec R ^ {kn} = \ vec r ^ k- \ vec r ^ n $$ Gdzie $ \ vec R_0 ^ k = \ sum_ {n = 1 n \ ne k} ^ N \ frac {m_n \ vec R ^ {kn}} {\ sum_ {p = 1 } ^ N m_p} $ Zdobądź wzór $$ \ frac {d ^ 2 \ vec R_0 ^ k} {ds ^ 2} = \ frac {d ^ 2 \ vec r ^ k} {ds ^ 2} $$ Całkując tę ​​równość, otrzymujemy równanie ruchu dla każdego z N ciał. $$ \ vec r_0 ^ k = \ frac {d \ vec r ^ k (0)} {ds} s + \ vec r ^ k (0) + \ vec R_0 ^ k ( s) - \ frac {d \ vec R ^ k (0)} {ds} s- \ vec R ^ k (0) = \ vec R_0 ^ k (s) $$ Ruch każdego ciała jest określany przez ruch środka bezwładności układu ciała pary. Warunki początkowe są określane na podstawie warunku $$ \ frac {d \ vec r ^ k (0)} {ds} s + \ vec r ^ k (0) = \ frac {d \ vec R ^ k (0)} {ds} s + \ vec R ^ k (0) \ tag {3} $$ Wybieramy układ współrzędnych $ \ sum_ {n = 1} ^ N m_n \ vec r ^ n (0) = 0, \ sum_ {n = 1} ^ N m_n \ frac { d \ vec r ^ n (0)} {ds} = 0 $ , a następnie warunek (3) jest spełniony identycznie. Otrzyma $ \ vec R ^ {kn} (0) = \ vec r ^ k (0) - \ vec r ^ n (0), \ frac {d \ vec R ^ {kn}} {ds} (0) = \ frac {d \ vec r ^ k} {ds} (0) - \ frac {d \ vec r ^ n} {ds} (0) $ energia pary trajektorii układu jest obliczana według wzoru $$ E_ {kn} = \ frac {m_k m_n} {\ sum_ {q = 1} ^ N m_q} (\ frac {d \ vec R ^ {kn}} {ds} (0)) ^ 2 / 2- \ frac {G m_k m_n} {| \ vec R ^ {kn} (0) |} + \ frac {M_ { kn} ^ 2 \ sum_ {q = 1} ^ N m_q} {2 m_k m_n | R ^ {kn} | ^ 2} = \ frac {m_k m_n} {\ sum_ {q = 1} ^ N m_q} (\ frac {d \ vec R ^ {kn}} {ds} (0)) ^ 2 / 2- \ frac {G m_k m_n} {| \ vec R ^ {kn} (0) |} + \ frac {m_k m_n } {2 \ sum_ {q = 1} ^ N m_q | R ^ {kn} (0) | ^ 2} (\ frac {d \ vec R ^ {kn}} {ds} (0) \ times \ vec R ^ {kn} (0)) ^ 2 $$ Moment impulsu jest określany przez iloczyn wektorowy wzoru $$ M_ {kn} = \ frac {m_k m_n } {\ sum_ {q = 1} ^ N m_q} \ frac {d \ vec R ^ {kn}} {ds} (0) \ times \ vec R ^ {kn} (0) $$ W tym przypadku skonstruowane rozwiązanie jest podstawiane do obliczonych trajektorii $ \ vec r_0 ^ k- \ vec r_0 ^ p $ planet i używane $ \ vec R ^ {kp} = r ^ kr ^ p $ , wtedy otrzymujemy $ \ vec R ^ {kp} = \ vec r_0 ^ k- \ vec r_0 ^ p $ . W rzeczy samej $$ r_0 ^ k-r_0 ^ p = \ sum_ {n = 1} ^ N \ frac {m_n \ vec R ^ {kn}} {\ sum_ {q = 1} ^ N m_q} - \ sum_ {n = 1} ^ N \ frac {m_n \ vec R ^ {pn}} {\ sum_ {q = 1} ^ N m_q} = \ sum_ {n = 1} ^ N \ frac { m_n (\ vec r ^ k- \ vec r ^ n)} {\ sum_ {q = 1} ^ N m_q} - \ sum_ {n = 1} ^ N \ frac {m_n (\ vec r ^ p- \ vec r ^ n)} {\ sum_ {q = 1} ^ N m_q} = \ vec r ^ k- \ vec r ^ p = R_ {kp} $$ Formuła uproszczona $$ \ vec r_0 ^ u = \ sum_ {n = 1} ^ N \ frac {m_n \ vec R ^ {un} (s)} {\ sum_ {q = 1} ^ N m_q | R ^ {un} (s) |} \ frac {p_ {un}} {1 + e_ {un} cos (\ varphi_ {un} - \ varphi_ {un0})} $$ Kąt jest obliczany według wzoru $ \ varphi_ {un} - \ varphi_ {un0} = \ int_0 ^ s \ frac {M_ {un}} {| \ vec R ^ {pn } (u) | ^ 2} du $ Wielkość $$ \ frac {R ^ {un} (\ varphi)} {| R ^ {un} (\ varphi) |} = \ frac {\ vec R ^ {un} (0)} {| \ vec R ^ {un} (0) |} cos \ varphi + \ frac {\ vec R ^ {un} (0) \ razy \ vec M_ {un}} {| \ vec R ^ {un} (0) \ times \ vec M_ {un} |} sin \ varphi $$ kąta $ \ varphi $ wspólna dla wszystkich trajektorii par odpowiada kierunkowi $ \ frac {R ^ {un} (\ varphi)} {| R ^ {un} (\ varphi) |} $ z pojedynczymi ortami $$ \ frac {\ vec R ^ {un} (0)} {| \ vec R ^ {un} ( 0) |}, \ frac {\ vec R ^ {un} (0) \ times \ vec M_ {un}} {| \ vec R ^ {un} (0) \ times \ vec M_ {un} |} $ $ . Mamy ostateczny wzór na trajektorię ciał $$ \ vec r_0 ^ u (\ varphi) = \ sum_ {n = 1} ^ N \ frac {m_n} {\ sum_ { q = 1} ^ N m_q} [\ frac {\ vec R ^ {un} (0)} {| \ vec R ^ {un} (0) |} cos \ varphi + \ frac {\ vec R ^ {un} (0) \ times \ vec M_ {nie}} {| \ vec R ^ {un} (0) \ times \ vec M_ {nie} |} sin \ varphi] \ frac {p_ {un}} {1 + e_ {un} cos \ varphi} $$ Ale z powodu braku nieskończonego rozwiązania nawet przez krótki czas konieczne jest założenie $$ | \ vec R ^ { un} (\ varphi) | = \ frac {p_ {un}} {1 + e_ {un} cos \ varphi} >0, e_ {un} <1 $$ To nakłada ograniczenie na wartość parametrów ; energia interakcji par musi być ujemna. W ten sposób wykluczone są ciała o małej masie, których wkład w energię układu nie jest znaczący. Zgodnie ze wzorem interakcji par, promień zmienia się zgodnie z prawem, patrz Landau Lifshits tom 1 $$ r = a (ecosh \ xi-1); t = \ sqrt {\ frac {a ^ 3} {G \ sum_ {q = 1} ^ N m_q}} ( esinh \ xi- \ xi), r \ ok c \ sqrt {\ frac {r_g} {a}} t $$ Po prostych przekształceniach ta formuła jest redukowana do $$ r = \ sqrt {V ^ 2 (0) + \ frac {(\ vec V (0) \ times \ vec R (0))) ^ 2} {R (0) ^ 2} - \ frac {2G \ sum_ {q = 1} ^ N m_q} {R (0)}} t $$ Ten liniowy wzrost promienia przesuwaśrodek ciężkości systemu.Cząstka poruszająca się w nieskończoność prowadzi do przenoszenia całego układu i przesuwania środka ciężkości.Druga prędkość kosmiczna odpowiada dodatniej wartości radicandy, czyli dodatniej energii układu.