Planeetan liike kiinteän Auringon ympäri

Planeetan liike Auringon ympäri voidaan yksinkertaisimmassa muodossa mallintaa pitämällä Aurinkoa kiinteänä ja ajattelemalla ainoaksi vaikuttavaksi voimaksi Auringon ja planeetan välinen gravitaatio. Todellisuudessa aivan näin ei ole: Planeetta ja Aurinko muodostavat kappaleparin, jossa molemmat vaikuttavat toisiinsa ja kumpikin liikkuu systeemin yhteisen massakeskipisteen suhteen. Liikkuvia kappaleita on myös useampia, Auringon lisäksi monia planeettoja, ja kaikki vaikuttavat toisiinsa. Yksinkertainen malli on suhteellisen hyvä approksimaatio, jos voidaan ajatella, että Auringon massa on merkittävästi planeetan massaa suurempi ja muut planeetat ovat niin etäällä, että niiden vaikutus voidaan jättää huomiotta.

Gravitaatiolain mukaan Auringon ja planeetan välinen vetovoima on

F = ((γ (Subscript[m, 1]) (Subscript[m, 2]))/(r^2)),

missä (Subscript[m, 1]) on Auringon massa ja (Subscript[m, 2]) planeetan massa, γ on gravitaatiovakio ja r on kappaleiden välinen etäisyys. Oletetaan, että Aurinko sijaitsee origossa ja planeetan paikkavektori on (OverBar[r]). Tällöin planeettaan kohdistuva voima suuntautuu origoon. Newtonin lain mukaiseksi vektorimuotoiseksi liikeyhtälöksi saadaan

(Subscript[m, 2]) (OverBar[r])''=-((γ (Subscript[m, 1]) (Subscript[m, 2]))/(r^2)) ((OverBar[r])/r),

Sievennettynä tämä on

(OverBar[r])''=-(α/((r^3))) (OverBar[r]),

missä α  =  γ (Subscript[m, 1]). Kyseessä on toisen kertaluvun vektorimuotoinen yhtälö.

Jos oletetaan, että planeetan rata sijaitsee tasossa (mikä ei ole itsestään selvää, mutta seuraa kyllä liikeyhtälöstä), voidaan tason koordinaatteja merkitä x ja y, jolloin paikkavektori on (OverBar[r])(t) = x(t) (OverBar[i]) + y(t) (OverBar[j]). Differentiaaliyhtälö hajoaa tällöin kahdeksi komponenttiyhtälöksi

x''=-((α x)/(((((x^2) + (y^2)))^(3/2)))),

y''=-((α y)/(((((x^2) + (y^2)))^(3/2)))).

Merkitsemällä x' = u ja y' = v voidaan tämä palauttaa ensimmäistä kertalukua olevaksi neljän yhtälön normaaliryhmäksi, joka seuraavassa syötetään laskentaa varten Mathematicaan.

Aluksi kuitenkin on syytä hävittää mahdollisista aiemmista laskuista jääneet muuttujat:

Remove["Global`*"]

Normaaliryhmä on

ryhma =  {x'[t]  ==  u[t],   y'[t]  ==  
          v[t],  u'[
            t]  ==  -α x[
                t]/(x[t]^2  +  y[t]^2)^(3/2),   v'[
            t]  ==  -α y[
                t]/(x[t]^2  +  y[t]^2)^(3/2)}

ja probleeman tuntemattomat funktiot

tuntemattomat =  {x[t], y[t], u[t], v[t]}

Vakio α on Auringon massan ja gravitaatiovakion tulo (yksikkönä N (m^2) / kg = (m^3)/(s^2)):

α =  (1.989*10^30)*(6.673*10^(-11))

Alkuehdot voidaan antaa vaikkapa jonkin planeetan mukaisesti. Maan keskietäisyys Auringosta on noin 150 miljoonaa kilometriä ja ratanopeus keskimäärin noin 30 kilometriä sekunnissa. Seuraavat alkuehdot ovat tämän mukaisia. Yksikköinä on käytettävä metrejä ja sekunteja kuten eo. vakiossakin:

alkuehto = {x[0] == 0, y[0] == 150*10^9, u[0] == 30*10^3, 
        v[0] == 0}

Luontevaa on laskea rata yhden Maan vuoden ajalta yksikkönä sekunti:

tmax =  365*24*60*60

Differentiaaliyhtälöryhmän ratkaisu numeerisesti:

rtk =  
      NDSolve[Flatten[{ryhma, alkuehto}], 
        tuntemattomat, {t, 0, tmax}]

Planeetan rata voidaan piirtää parametrikäyränä; kuvioon lisätään Aurinkoa osoittava täplä:

rata =  {x[t], y[t]} /. First[rtk]

planeetta =  
      ParametricPlot[rata, {t, 0, tmax}, 
        DisplayFunction -> Identity]

aurinko =  
      Graphics[{PointSize[0.1], Point[{0, 0}]}]

Show[planeetta, aurinko, AspectRatio -> Automatic, 
      DisplayFunction -> $DisplayFunction]

Ratanopeuden vaihtelusta voidaan myös piirtää kuvaaja:

ratanopeus =  
      Sqrt[u[t]^2 + v[t]^2] /. First[rtk]

Plot[ratanopeus, {t, 0, tmax}, 
      PlotRange -> {0, Automatic}]

Tehtävä

Tutki, miten alkuehdossa annettu nopeus vaikuttaa rataan. Miten käy, jos nopeutta a) pienennetään, b) suurennetaan, c) sen suuntaa muutetaan? Mitä tapahtuu, jos nopeuden molemmat komponentit ovat = 0? Taivaankappaleiden radat tässä yksinkertaisessa mallissa voivat olla ellipsejä (erikoistapauksena ympyröitä), paraabeleja ja hyperbelejä. Saatko alkuehtoja muuttamalla esiin eri ratatyypit?

SKK & JP 03.05.2001