Planetbaner kommer i ulike former
Dersom du trodde at de analytiske banene vi har var nok, har jeg en d?rlig nyhet.. Universet er dessverre ikke s? enkelt, det hadde jo ikke v?rt noen utfordring for oss fysikere! Det m? kraftigere lut til, vi m? nemlig simulere planetbanene v?re tidssteg for tidssteg.
Jeg h?rer at du tenker "Men har vi ikke akkurat regnet ut planetbanene med de analytiske uttrykkene? Hvorfor m? vi simulere?", men frykt ikke, jeg skal oppklare: de analytiske uttrykkene gir oss riktig form p? ellipsebanen, men gir oss ikke riktig posisjon p? dem. Simuleringen skal gi oss posisjon av planetene som funksjon av tid, mens det analytiske uttrykket gir oss avstanden fra brennpunktet til ellipsebanen som en funksjon av vinkel. Dermed kan den analytiske banen og simulerte banen ligge rotert feil i forhold til hverandre, men fortsatt ha samme form p? ellipsebanen. Sitter dette litt l?st er det kanskje lurt ? g? tilbake til forrige post og sjekke det analytiske uttrykket for ellipsebanen. Det kan vi fikse ved ? rotere den analytiske ellipsebanen, men vi trenger fortsatt posisjonen av planeten som en funksjon av tid i stedet for vinkel. Vi er derfor n?dt til ? simulere det. Viktige antakelser vi ogs? kommer til ? gj?re her og i senere bloggpost er:
- Kun stjerna p?virker planetene med kraft, og planetene p?virker ikke hverandre med krefter.
- Vi ser bort ifra relativistiske effekter.
- I dette innlegget ser vi bort ifra kraften p? stjerna fra de andre planetene. I en senere bloggpost skal vi se p? stjernas bevegelse ogs?.
- Vi antar at alle planetbanene ligger i samme xy-plan, og ignorerer derfor z-aksen. Da jobber vi i 2 dimensjoner.
- Alle planeter roterer om sola og z-aksen (om seg selv) mot klokka.
- Andre legemer i solsystemet (som m?ner) blir sett bort ifra.
S?, hvordan skal vi g? frem for ? gj?re dette? F?rst vil jeg begynne med ? si noe om hvordan vi definerer koordinatsystemet n?r vi simulerer, som vil v?re et koordinatsystem med sentrum av sola i origo, alts? posisjon \((0,0)\) (slik som i Figur 1). Det betyr at simulerer vi planetbanene i et 2D-plan, xy-planet og ikke et 3D-rom. Det er kanskje greit ? avklare hvor lenge vi ?nsker ? kj?re simuleringen v?r. Vi ?nsker ? simulere planetbanene for minst 20 oml?p av hjemplaneten, (planet 1), om sola i systemet v?rt. Med andre ord m? vi regne ut hjemplanetens oml?pstid rundt sola i solsystemet (dette kan vi f.eks gj?re med Keplers 3. lov som vi kommer til ? utlede p? neste bloggpost). Dette er bare for ? s?rge for at vi simulerer alle planetbanene lenge nok, ogs? tar det jo tid ? reise gjennom solsystemet ogs? tross alt. Videre trenger vi et uttrykk for akselerasjonen til hver planet. Det finner vi lett ved ? bruke Newtons 2. lov. Neste steg blir ? bruke en kjent teknikk for ? finne posisjoner og hastigheter numerisk, numerisk integrasjon.
Som du kanskje husker, brukte vi dette da vi simulerte partikler i rakettmotoren tidligere. Det betyr at vi ogs? blir n?dt til ? finne et uttrykk for akselerasjonen til hver planet. Men denne gangen m?ter vi p? et problem: metoden vi brukte sist blir rett og slett ikke n?yaktig nok i simuleringen vi skal gj?re n?. Tidligere brukte vi noe kalt Eulers metode, som er god nok n?r vi ser p? korte intervaller slik som i rakettmotoren, men ikke god nok i v?rt tilfelle n?. I tillegg er ikke Eulers metode det vi kaller energibevarende, alts? f?r/mister systemet energi over tid n?r det ikke burde det (som i Figur 2). Hva som gj?r en metode energibevarende er jobb for en numeriker, s? jeg tar med ikke bryet med ? forklare det her, men det er likevel viktig ? vite om metoden vi skal bruke har denne egenskapen. Vi kunne kanskje brukt Eulers metode's storebror, Euler-Cromer! Han er jo energibevarende! Her sparer jeg meg selv tid ved ? ikke implementere det, fordi jeg vet allerede at denne metoden heller ikke er presis nok over tid..
Her trengs det andre astronomiske boller. Vi trenger mer presise, energibevarende astronomiske boller. Vi leter rundt i verkt?ykassa, og finner et redskap kalt leap-frog. Dette er en energibevarende og presis nok metode for numerisk integrasjon, som vi kommer til ? bruke i simuleringen v?r. Oppsummert + litt ekstra er dette overordnet hvordan vi l?ser problemet:
- Vi finner tiden til x antall oml?p for hjemplaneten om sola med Keplers 3. lov, \(P^2 = \frac{4\pi^2a^3}{G(m_1 + m_2)}\) og med x som minimum 20 oml?p.
- Vi finner uttrykket for akselerasjonen til planeten vi skal simulere med Newtons 2. lov, \(\Sigma\vec{F} = m\vec{a}\).
- Vi bruker data om initialposisjoner og initialhastigheter for planetene, ogs? integrerer vi opp posisjonen og hastigheten til alle planetene i solsystemet over tid med den numeriske metoden leapfrog over tiden for 20 oml?p om sola med hjemplaneten. Videre kan vi visualisere og teste dataen vi har f?tt.
Data for masser (planeter og stjerne), store halvakser, initialposisjoner og initialhastigheter for hver planet som vi kommer til ? bruke i denne posten kan du finne p? de forrige bloggpostene, Solsystemet v?rt og Analytiske planetbaner, s? jeg lager derfor ingen tabell i denne posten. MERK at vi i denne simuleringen kommer til ? bruke enheter som AU (astronomiske enheter) for lengde, \(\odot \) (solmasser) for masse, og yr (?r) for tid, som betyr at kjente konstanter som Newtons gravitasjonskonstant G eller \(\textbf{$\gamma$}\) (avhengig av hva man pleier ? bruke) ogs? vil se litt annerledes ut p? grunn av enhetene. f.eks er \(G = 4\pi^2 AU^3 yr^{-2}M^{-1}\) der \(M=M_{planet} + M_{sol}\) , generelt er \(M\) total masse i banesystemet vi ser p?. Vi bruker i denne oppgaven at \(M_{sol}=M\odot\) (planetmassen er liten i forhold, derav neglisjerbar), alts? ikke massen til stjerna i solsystemet v?rt, men en solmasse (massen til sola i jordas solsystem).
Hoppe frosk - numeriske datafrosker
Leapfrog er en metode for numerisk integrasjon, p? lik linje med metoden vi har brukttidligere. Ideen med numerisk integrasjon er ? regne ut hva den neste posisjonen og hastigheten er ved ? g? et bittelite tidssteg videre ved ? gj?re approksimasjoner. Derfor er det veldig viktig at tidssteget v?rt er lite nok. Ellers blir det en d?rlig approksimasjon. Andre eksempler p? numerisk integrasjon er Eulers-Cromer eller Eulers metode beskrevet i et tidligere innlegg. Vi vet hvor planetene befinner seg i en tid t=0, hvilken hastighet de har der, og akselerasjonen kan vi finne med Newtons lover. Da har vi alt vi trenger for ? integrere numerisk. Denne metoden er liiitt mer komplisert enn Eulers metode som vi har brukt tidligere, men kan fortsatt utledes p? en enkel m?te. Jeg kommer ikke til ? gj?re det her, men du kan finne det i et vedlegg linket her om du er interessert. Metoden for leapfrog tar flere former, s? jeg skal vise dere to av dem - og pr?ve ? gi litt intuisjon.
\((1) \begin{align} a_i &= A(x_i) \\ v_{i+\frac{1}{2}} &= v_{i-\frac{1}{2}} + a_i\Delta t \\ x_{i+1} &= x_i + v_{i+\frac{1}{2}}\Delta t \\ \end{align}\) \(????????????????(2) \begin{align} a_i &= A(x_i) \\ x_{i+1} &= x_i + v_i\Delta t + \frac{1}{2}a_i\Delta t^2 \\ a_{i+1} &= A(x_{i+1})\\v_{i+1} &= v_i + \frac{1}{2}(a_i + a_{i+1})\Delta t \end{align}\)
Hva er det som skjer i denne metoden? I stedet for ? g? hele steg om gangen, tar vi et mellomsteg p? hastigheten. For ? illustrere hva som skjer kan vi tegne opp en linje med indeks, posisjon og fart. Posisjon vil da representeres som gr?nn frosk, og hastighet som rosa frosk:
Det ser nesten ut som at leapfrog hopper bukk, og det er jo det den gj?r. Leapfrog oversatt direkte til norsk er jo hoppe bukk! Illustrert p? Figur 3 ser man hvordan leapfrog l?per p? indeksene. Nummereringen i illustrasjonen kan v?re litt misvisende, fordi rosa frosk er ferdig med hoppet sitt fra (1) til (3) f?r gr?nn frosk begynner ? hoppe i (2), mens i illustrasjonen kan det virke som om det er motsatt. Alts? blir rekkef?lgen p? regningen (1) til (3), (2) til (4) til (3) til (5), (4) til (6) osv. Det som er verdt ? merke seg er at vi ogs? m? regne ut akselerasjonen i det neste tidssteget underveis, som er annerledes enn fra tidligere metoder vi har sett p?.
I v?rt tilfelle vil posisjon og hastighet v?re tomme arrays (litt om arrays her om du ikke husker) som inneholder initialposisjon og initialhastighet for en tid t=0 til en planet i x og y-retning, som vi fyller opp med posisjoner og hastigheter planetene vil ha med tiden ved bruk av leapfrog. For numeriske beregninger egner formel (2) for leapfrog seg bedre (lettere ? implementere). Det er denne vi kommer til ? bruke senere for leapfrog.
Simulate orbits - 101
Fra punktlisten ovenfor vet vi n? akkurat hva vi m? gj?re. Men hvordan skal vi gj?re det? Fra punkt 1) st?r det at vi m? finne ut hvor lenge vi skal simulere, som er minimum 20 oml?p om sola for hjemplaneten v?r. Vi velger 21 oml?p s? vi er sikre p? at vi har ekstra god tid, uten at det skal ta for lang tid ? simulere. Vi finner s? tiden hjemplaneten bruker p? ett oml?p om sola, ved ? bruke Keplers 3. lov som er gitt som:
\(P = \sqrt{\frac{4\pi^2a^3}{G(m_1 + m_2)}}\approx 1.8yr\)
der \(a\) er store halvakse i ellipsebanen til hjemplaneten, \(G\) er Newtons gravitasjonskonstant, \(m_1\) er massen til solen i solsystemet og \(m_2\) er massen til hjemplaneten. \(P \) er da oml?pstiden til hjemplaneten v?r gitt i ?r, p? grunn av enhetene vi bruker. For ? videre f? tiden hjemplaneten bruker om sola 21 ganger er alt vi trenger ? gj?re ? gange \(P\) med 21, som gir oss total oml?pstid som
\(\text{total time} = 21\cdot P\)
N? som vi har total tid for simulering, kan vi g? videre til punkt 2), hva er akselerasjonen. Vi antar at eneste kraften som virker p? en planet er kraften fra sola, og vi neglisjerer krefter en planet f?r fra de andre planetene i solsystemet siden de blir veldig sm? i forhold (tenk sola har en masse med potens \(10^{30}\) mens jorda har en masse med potens \(10^{24}\), som betyr at kraften fra sola er rundt 1 million ganger st?rre). Dette gj?r det litt enklere for oss, og kraften fra sola p? en planet er gitt av Newtons gravitasjonslov:
\(\Sigma\vec{F} = m_{planet}\vec{a} = -G\frac{M_{sun}m_{planet}}{r^2}\frac{\vec{r}}{r}\)
der \(\vec{r}\)er posisjonsvektoren til planeten vi ser p? som peker fra sentrum av sola til sentrum av planeten, og \(r = |\vec{r}|\). De andre st?rrelsene b?r v?re trivielle. L?ser vi for akselerasjonen f?r vi uttrykket:
\(\vec{a} = -G\frac{M_{sun}}{r^2}\frac{\vec{r}}{r}\)
MERK at akselerasjonen blir det samme for alle planetene. Det som ogs? er verdt ? merke er at akselerasjonen er avhengig av posisjonen til planeten, som betyr vi blir n?dt til ? oppdatere uttrykket underveis og p? slutten av hver iterasjon med leapfrog.
Med data for initialbetingelser fra Solsystemet v?rt og Analytiske planetbaner er vi n? klare til ? anvende leapfrog og f? de faktiske planetbanene som funksjon av tid. Vi velger ? bruke omtrent 11000 tidssteg for hvert ?r vi simulerer (vi runder av til n?rmeste heltall fordi desimale tiddsteg ikke gir mening), som betyr at antall tidssteg blir \(\text{timesteps}=int(\text{total time} * 11000)\). Da f?r vi et tidssteg \(dt\) p? \(dt = \text{total time} / \text{timesteps}\), som kommer til ? bli noe rundt ca. \(\frac{1}{11000}\approx 9\cdot10^{-5}\) (ikke n?yaktig med noe rundt den st?rrelsen).
Solsystemet simulert
Dersom vi plotter x-posisjonen og y-posisjonen fra dataen leapfrog ga oss ender vi opp med dette plottet:
Tiden dette plottet gjelder for er 21 oml?p rundt sola for planeten v?r. Fra Keplers 3. lov fant vi oml?pstiden til hjemplaneten v?r, \(P = \sqrt{\frac{4\pi^2a^3}{G(m_1 + m_2)}} \approx1.8yr\), som betyr at tiden vi har simulert er ca. 37.8 ?r. Sola syntes ikke p? plottet, men ville befunnet seg i origo. Banene ser lovende ut, men det er ikke bra nok for ? v?re sikre p? om det stemmer, eller hva slags baner vi har. De ser f.eks veldig sirkul?re ut, men kan fortsatt v?re ellipser. Dette er ikke nok informasjon. Da er det kjekt ? teste resultatene, og det har vi tenkt til ? gj?re med Keplers lover. Men f?rst skal vi utlede dem!