Gravitacijska simulacija

Gravitacijska sila med dvema telesoma je seveda dovolj, da preprosto simuliramo poljubno število teles (tu za ceno natančnosti lahko močno pridobimo na času), saj poljubno število teles lahko razdelimo v pare. Za le dve telesi obstajajo tudi analitične rešitve, za več teles pa se moramo poslužiti numeričnim metodam.

Toda 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 = 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. 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}

Primerjava integracijskih metod¶

Napišimo sedaj zelo preprost simulator, ki prikaže razliko med različnimi integracijskimi metodami. Vzamemo sistem Sonce-Zemlja in pogledamo kako so različne metode natančne in stabilne (pričakujemo namreč, da Zemlja kroži okoli Sonca in se ohranjajo vse količine).

from simulatorv2 import *
from numpy import sqrt

G = NewtonianGravity().G

b1 = Presets.sun(Vec2)

b2 = Presets.earth(Vec2)
b2.position = Vec2(1,0)
b2.velocity = Vec2(0,sqrt(G))

sim = Simulator(
    bodies=[b1,b2],
    timestep=0.001,
    steps=1000,
    diagnostics=[
        KineticEnergy(),
        PotentialEnergy(),
        TotalEnergy(),
        AngularMomentum()
    ]
)


leapfrog_result = sim.run("Leapfrog", "blue")
sim.reset()

sim.integrator = ExplicitEuler()
sim.reset()
explicit_euler_result = sim.run("Eksplicitni Euler", "red")

sim.integrator = SemiImplicitEuler()
sim.reset()
semi_euler_result = sim.run("Semi Implicitni Euler", "green")

sim.integrator = VelocityVerlet()
sim.reset()
velocity_verlet_result = sim.run("Hitrostna Verlet", "purple")

viewer = SimulationViewer(
    [
        leapfrog_result,
        explicit_euler_result,
        semi_euler_result,
        velocity_verlet_result
    ],
    ViewerConfig(
        diagnostics=[
            DiagnosticView(
                "Celotna Energija",
                ["vsi"]
            ),
            DiagnosticView(
                "Vrtilna Količina",
                ["vsi"]
            )
        ],
        camera_mode="follow",
        follow_body="Zemlja",
        trail_length=10,
        view_size=0.01,
        references = [
            ReferenceCircle()
        ]
    )
)

viewer.show()

Loading...

V zgornjem rezultatu vidimo, da je najbolj stabilna in natančna metoda Leapfrog, zato bomo od tu naprej uporabljali le to. Ostale metode niso tako slabe, toda vidimo, da energija za ostale metode ni periodična, torej ali bo v prihodnosti eksplodirala ali zadonila.

Aproksimacijske metode¶

Za velika števila teles se simulacija seveda zelo upočasni (program mora iti čez vse možne pare). Zato se moramo v takih primerih poslužiti tudi optimizacijam za ceno natančnosti.

Barnes-Hut aproksimacija¶

Ta sloni na sekanju prostora na podprostore, glede na kopičenje teles (torej več kot je teles na kupu, več podprostorov bo). Če je nek podprostor dovolj oddaljen od nekega telesa, ga to telo upošteva kot eno telo z eno zelo veliko skupno maso (kot mi na primer opisujemo mase kopic ali podobno).

Metoda mreže delcev¶

Še ena možnost je, da prostor diskretiziramo in naredimo mrežo delcev (kot v simulaciji tekočin). Tam lahko rešimo Poissonovo enačbo (druga oblika Newtonove gravitacijske enačbe) in s pomočjo potenciala premikamo maso po celicah (kot simulacija tekočin, le da tu delujejo drugi zakoni).

Testna simulacija¶

Prepričajmo se sedaj, da simulacija deluje pričakovano. To naredimo s preprosto simulacijo, za katero lahko predvidevamo rezultat - obravnavajmo hipotetični primer, da bi Zemlja bila kopija Sonca, ter da drugih teles v Osončju ne bi bilo. Po Keplerjevem zakonu lahko zelo preprosto napovemo:

T^2\propto\frac{a^3}{M}\to T\propto\frac{1}{\sqrt{M}}=\frac{1}{\sqrt{2}}\approx 0.707.

V tej enačbi smo upoštevali, da so ostale spremenljivke v formuli ostale enake. Pričakujemo torej, da bo perioda takih zvezd bila enaka približno 0.707 leta, kar je v resnici produkt:

P = T\cdot S,

kjer so $P$ perioda, $T$ časovni korak in $S$ število korakov. Poskrbeti moramo torej, da je produkt časovnega koraka in števila korakov enak željeni periodi, kar je prikazano v spodnji kodi.

from simulatorv2 import *
from numpy import sqrt

G = NewtonianGravity().G
v_mag = sqrt(G * 2 / 4)

b1 = Presets.sun(Vec2)
b1.name = "Sun 1"
b1.position = Vec2(-0.5, 0)
b1.velocity = Vec2(0, v_mag)

b2 = Presets.sun(Vec2)
b2.name = "Sun 2"
b2.position = Vec2(0.5, 0)
b2.velocity = Vec2(0, -v_mag)

sim = Simulator(
    bodies=[b1, b2],
    timestep=0.001,
    steps=707,
    diagnostics=[
        TotalEnergy(),
        AngularMomentum()
    ]
)

result = sim.run()

viewer = SimulationViewer(
    result,
    ViewerConfig(
        diagnostics=[
            DiagnosticView(
                "Celotna Energija",
                ["vsi"]
            ),
            DiagnosticView(
                "Vrtilna Količina",
                ["vsi"]
            )
        ],
        trail_length=100,
        view_size=0.6,
        references = [
            ReferenceCircle(r=0.5, opacity=0.2)
        ]
    )
)

viewer.show()

Loading...

Kot vidimo na zgornji animaciji, je simulacija zelo natančna, predvsem zaradi naslednjih stvari:

privzeto uporabljamo Leapfrog metodo, ki je najbolj natančna, toda tudi najbolj potratna iz navedenih,
perioda simulacije se sklada z napovedano periodo - telesi sta res po 0.707 leta v (približno) začetnih pozicijah,
celotna energija nakazuje periodično spreminjanje (ne eksplodira),
celotna energija in vrtilna količina imata obe zelo majhni amplitudi (razmerje med amplitudo in celotno vrednostjo je za energijo približno 10^-8 za vrtilno količino pa celo 10^-13).