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!

źródło