Numerics: Jak renormalizować następujące ODE

To pytanie dotyczy bardziej numerycznego rozwiązania problemu.

W małym projekcie chciałem zasymulować ruch koorbalny Janusa i Epimetheusa. Jest to w zasadzie problem trzech ciał. Niech wybiorę Saturna, który ma być ustawiony na początku $r_1$ i $r_2$ być odpowiednio wektorami lokalizacji Janusa i Epimetheusa. Ponieważ efekt występuje, gdy Janus i Epimetheus są bardzo blisko siebie, wybrałem współrzędne względne dla lepszej rozdzielczości, tj $r=r_1-r_2$ i $R=r_1+r_2$ . Teraz otrzymuję następujące równania ruchu:

\frac{{re}^{2)}}{re t^{2)}} (\binom{R}{r}) = - sol (m_{2)} \pm m_{1}) \frac{R}{R^{3)}} - 4 M. sol (\frac{r + R}{(r + R)^{3)}} \mp \frac{r - R}{(r - R)^{3)}})

$\frac {d^2}{dt^2} \binom{R}{r} = - G (m_2\pm m_1) \frac R {R^3} - 4 M G \left(\frac {r+R}{(r+R)^3} \mp \frac {r-R}{(r-R)^3}\right )$

gdzie odpowiada masom księżyców, jest masą Saturna, a stałą grawitacyjną. Problem pojawia się, gdy próbuję rozwiązać to numerycznie. Należy sobie radzić z wartościami o zupełnie różnych wielkościach, tj. i . I , są regionami 0 do 150.000. $m_i$ $M$ $G$ $M \sim e^{28}$ $m_i \sim e^{17}$ $r$ $R$

Szczerze mówiąc, nie jestem pewien, czy jest to forum dyskusyjne na temat takich problemów numerycznych.

Więcej informacji:

Kod jest napisany w Matlabie i do uzyskania wyniku używam standardowego solvera ODE. Jednak to się psuje, ponieważ wielkości kroku nie można zmniejszyć przy precyzji maszyny. (Uważam, że nie jest to zaskakujące, ponieważ trzeba poradzić sobie z już wspomnianymi rzędami wielkości).

ode

— GertVdE
źródło

Czy przeprowadzasz tę symulację w jednostkach SI? Jako minimum powinieneś podzielić wszystko przez pewien współczynnik , abyś mógł wyeliminować kilka rzędów wielkości.

G * m_{2}

$G*m_{2}$

Cześć, ja to, ale wciąż nie działa ... Te same problemy występują jak poprzednio. :(

Musisz ustawić swoją jednostkę masy na jedną z mas księżycowych, a jednostki długości / czasu, aby ustawić rzeczy na 1. Nic nie powinno być mniejsze niż 1/100, jeśli dobrze to napiszesz. Nie ma potrzeby korzystania z solvera bez recepty. Napisz kod, aby zrobić to sam, w którym kontrolujesz wielkość kroku. W przypadku kolizji mogą wystąpić awarie ze względu na tego rodzaju potencjały, w których solver będzie próbował zmniejszyć stopień aż do zbieżności, a przy zderzeniu nie ma zbieżności. Musisz upewnić się, że orbity nie są współliniowe, więc musisz zobaczyć symulację. Nie możesz uzyskać takiej odpowiedzi, jaka jest.

— Ron Maimon

Unikaj skrótów w tytule. DGL = Differentialgleichung?

Jakiego standardowego solvera ODE używasz?

— Geoff Oxberry

Odpowiedzi:

Twoje obecne podejście rujnuje stabilność liczbową; w rzeczywistości prawdopodobnie stracisz w ten sposób rozdzielczość.

Przyjmij jako współrzędne dla każdego satelity jego zmienne Keplera i kąt płaszczyzny zawierającej pozycję satelity, prędkość i początek. Równania różniczkowe przy braku interakcji między satelitami są wtedy banalnie proste i tylko interakcja staje się nieco skomplikowana. Ponieważ interakcja jest niewielka, jeśli satelity są daleko, wynikowa dynamika powinna być stabilna numerycznie.

— Arnold Neumaier
źródło

Zamiast używać „klasycznego” (sztywnego) solvera ODE, można użyć dedykowanych algorytmów do geometrycznej integracji numerycznej. Zobacz na przykład tę książkę i kody DNB, które można znaleźć na stronie internetowej Ernsta Hairera .

— GertVdE
źródło

Co powiesz, jeśli masz trzy kroki w swojej symulacji:

zaktualizuj pozycję Janusa, obliczając siłę Janusa - Saturna.
zaktualizuj pozycję Epimetheusa, obliczając siłę Epimetheus - Saturn.
zaktualizuj pozycję Janusa i Epimetheusa, obliczając siłę Janusa - Epimetheusa.

Być może przy użyciu dokładniejszych czasów dla # 3.

Nie jestem pewien, czy to pomoże. Przypuszczam, że prawdziwym problemem jest to, że wielkość siły jest różna w przypadku Księżyc - Księżyc i Księżyc - Saturn, chyba że księżyce są blisko?

Alternatywnie:

jeśli księżyce się zamkną, oblicz przybliżone księżyce - siłę Saturna, używając ich wektora środka masy i zaktualizuj obie pozycje tym samym wektorem.
jeśli są daleko od siebie, zaktualizuj je osobno.
jak wcześniej.

Powodzenia!