Problem N teles - Sim-Astro

V splošnem imamo seveda več kot le dve telesi. Gibanje vseh teles ne moremo izračunati analitično, zato se moramo poslužiti numerični integraciji. Zaradi popolnega nadzora in izgleda smo se odločili, da bomo napisali svoj integrator in prikaz rezultata v HTML s pomočjo javascript-a, za čim hitrejši in nemoteni prikaz.

Astronomski pojavi običajno trajajo dolgo časa (na primer perioda orbite Zemlje okoli Sonca je po definiciji enaka enemu letu). Da v mednarodnem sistemu enot lahko simuliramo astronomske pojave torej potrebujemo zelo veliko korakov. Zato uporabimo astronomski sistem enot, kjer ena časovna enota predstavlja eno leto (po definiciji $86400\cdot 365.25\ \text{s}$ ), masna enota predstavlja maso Sonca (po definiciji $1.98892\cdot 10^{30}\ \text{kg}$ ) in dolžinska enota predstavlja astronomsko enoto (po definiciji $149.597.870.700\ \text{m}$ ). V takem sistemu zaradi tretjega Keplerjevega zakona v primeru Zemlje velja še:

T^2 = \frac{4\pi^2}{GM}a^3\to T = 1,\ a = 1,\ M\approx 1\to G\approx 4\pi^2.

Integracijske metode¶

Izraz za silo poznamo točno, iz izraza za silo pa točno poznamo pospešek telesa. Toda tu sedaj nastane problem v temu, kako upoštevamo ta pospešek - čas namreč simuliramo po korakih, in te koraki so običajno veliki, zato da rezultat simulacije dobimo v zmernem času. Poslužiti se moramo torej aproksimacijam, ki jih je veliko.

Eksplicitna Eulerjeva metoda¶

Je najbolj enostavna in posledično najbolj nezanesljiva metoda. Ta pravi, da upoštevamo klasične enačbe za gibanje telesa in upoštevamo diskretni korak:

\begin{aligned} \vec{v}_{n+1} &= \vec{v}_n + \vec{a}_n\Delta t,\\ \vec{r}_{n+1} &= \vec{r}_n + \vec{v}_n\Delta t. \end{aligned}

Takoj lahko uvidimo zakaj je ta način neprimeren. Predstavljajmo si satelit, ki kroži okoli Zemlje - njegova hitrost je vedno tangencialna na krožnico. Če ga premaknemo za vektor hitrosti, bo satelit že po prvi iteraciji padel dol iz krožnice.

Semi-implicitna Eulerjeva metoda¶

Ta je zelo podobna eksplicitni metodi, le da za premik upošteva hitrost v naslednjem koraku:

\begin{aligned} \vec{v}_{n+1} &= \vec{v}_n + \vec{a}_n\Delta t,\\ \vec{r}_{n+1} &= \vec{r}_n + \vec{v}_{n+1}\Delta t. \end{aligned}

Če pomislimo na enakem primeru kroženja satelita okoli Zemlje, je to res nekoliko bolje. Satelit bo v naslednji točki res imel hitrost, ki kaže bolj usmerjeno proti središču Zemlje (odmik od krožnice bo torej manjši).

Pozicijska Verlet metoda¶

Ta je že nekoliko bolj neintuitivna, toda jo je zelo lahko izpeljati. Ideja je da novo pozicijo izračunamo kar neposredno iz pospeška, ki ga poznamo točno:

\frac{\Delta^2}{\Delta t^2}\vec{r}_n^2 = \frac{\frac{\vec{r}_{n+1} - \vec{r}_n}{\Delta t} - \frac{\vec{r}_n - \vec{r}_{n-1}}{\Delta t}}{\Delta t}\to\vec{r}_{n+1} = 2\vec{r}_n - \vec{r}_{n-1} + \vec{a}_n\Delta t^2.

Pri tej metodi torej sploh ne upoštevamo hitrosti, zato je ne bomo simulirali (metoda ni uporabna za natančnost, morali pa bi zelo spremeniti strukturo simulatorja, kar se ne splača za samo eno primerjavo). Začetni pogoj sta sedaj dve točki, kjer drugo lahko posredno izračunamo iz hitrosti:

\vec{r}_1 = \vec{r}_0 + \vec{v}_0\Delta t + \frac{1}{2}\vec{a}_0\Delta t^2.

Hitrostna Verlet metoda¶

Ta je še nekoliko boljša kot pozicijska Verlet metoda, saj izračuna še vmesni korak, s pomočjo katerega lahko boljše oceni nove parametre telesa:

\begin{aligned} \vec{r}_{n+1} &= \vec{r}_n + \vec{v}_n\Delta t + \frac{1}{2}\vec{a}_n\Delta t^2, \\ \vec{v}_{n+1} &= \vec{v}_n + \frac{1}{2}(\vec{a}_n + \vec{a}_{n+1})\Delta t. \end{aligned}

Toda kot je razvidno iz formule, hitrost nastavimo glede na pospešek, ki ga moramo ponovno izračunati v novi poziciji (torej gre za bolj časovno potratno metodo).

Leapfrog metoda¶

Ta je zelo podobna hitrostni Verlet metodi, le da zanjo ni potrebno shraniti prejšnjih pospeškov, ampak vmesne hitrosti:

\begin{aligned} \vec{v}_{n+\frac{1}{2}} &= \vec{v}_n + \frac{1}{2}\vec{a}_n\Delta t, \\ \vec{r}_{n+1} &= \vec{r}_n + \vec{v}_{n+\frac{1}{2}}\Delta t, \\ \vec{v}_{n+1} &= \vec{v}_{n+\frac{1}{2}} + \frac{1}{2}\vec{a}_{n+1}\Delta t. \end{aligned}

Runge-Kutta metoda¶

Ta je nekoliko podobna Leapfrog in hitrostni Verlet metodi, le da ima več vmesnih izračunov (zato je tudi najbolj potratna). Metoda sloni na izračunu naslednjega stanja s pomočjo odvoda in nekega poljubnega koraka $h>0$ :

\frac{\text{d}y}{\text{d}t} = f(t, y),\ y(t_0) = y_0.

Naslednjo vrednost lahko izračunamo z več vmesnimi izračuni. Na primer metoda RK4 ima 4 vmesne izračune:

\begin{aligned} y_{n+1} &= y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4),\\ k_1 &= f(t_n, y_n),\\ k_2 &= f\left(t_n + \frac{h}{2}, y_n + k_1\frac{h}{2}\right),\\ k_3 &= f\left(t_n + \frac{h}{2}, y_n + k_2\frac{h}{2}\right),\\ k_4 &= f(t_n + h, y_n + hk_3). \end{aligned}

Da izračunamo neko vrednost v naslednjem trenutku torej potrebujemo še njen odvod. Za vektorje to seveda velja po komponentah, zato vse skupaj lahko zapišemo v vektorski obliki, zanimata pa nas tako pozicija kot hitrost:

\vec{y} = \begin{pmatrix}\vec{r} \\ \vec{v}\end{pmatrix}\to f(t, \vec{y}) = \begin{pmatrix}\vec{v} \\ \vec{a}\end{pmatrix}.

import sys
from pathlib import Path

BOOK_DIR = Path.cwd().parent

if str(BOOK_DIR) not in sys.path:
    sys.path.insert(0, str(BOOK_DIR))

from simulatorv3.plotter import show_simulation
from simulatorv3.bodies import Presets
from simulatorv3.vectors import Vec2

earth = Presets.earth()
earth.position = Vec2(1.0, 0.0)
earth.velocity = Vec2(0.0, 1.0)

sun = Presets.sun()
sun.position = Vec2()
sun.velocity = Vec2()

show_simulation(
    bodies=[sun, earth],
    filename="N_body_problem.html"
)

PosixPath('/srv/simastro/book/simulations/N_body_problem.html')