Resultater, testing og diskusjon av numeriske baner

Resultat, analyse (og litt diskusjon her ogs?)

Som vist i tidligere bloggpost fikk vi dette plottet som resultat fra simulering av planetbaner:

Vi kan aldri v?re helt sikre p? resultatene v?re uten ? teste dem. Vi har dermed tenkt til ? teste planetbanene v?re med Keplers lover. Vi bruker Keplers 2. lov for ? sjekke om arealet mellom 2 ulike tidsintervaller (et n?rt aphelium, et n?rt perihelium) og sjekker om de blir like, som de skal if?lge K2L (MERK at i plottet ovenfor er de numerert fra 0-7, men noen ganger er blir forskj?vet til 1-8 n?r vi skriver). Vi har at \(P\approx 1.8yr\) er oml?pstiden for hjemplaneten (planet 0) 1 gang om solen. Da blir tidsintervallene vi har valgt \(T1 \in \left(0, \frac{P}{10}\right)\) for areal 1, og \(T2 \in \left(\frac{P}{2}, \frac{P}{2} + \frac{P}{10}\right)\), for areal 2. Tiden er valgt slik at de er ca. p? motsatt side av banen, vilk?rlig sm? men sm? nok til at de ikke n?rmer seg hverandre s? vi f?r effekten av at det ene arealet er n?rmere aphelium og vice versa perihelium. Vi vet heller ikke hvilket av intervallene som er p? aphelium eller perihelium. Regner vi ut K2L f?r vi f?lgende resultat:

Det vi ser her er at arealene vi har f?tt er utrolig like. Tar vi en titt p? tallene fra arealene ser vi at tallene er helt like opp til syvende desimal! Den relative feilen forteller oss ogs? at feilen mellom de to tallene relativt til et av dem er bitte-lite! Om du ikke vet hvordan man regner relativ feil, s? har vi gjort det slik: \(\frac{|A2 - A1|}{A1}\), som betyr hva forskjellen mellom de to tallene er relativt til st?rrelsen p? et av tallene. Vi ser ogs? at den minste av arealene er \(99.9\%\) lik den andre.

Det er ogs? interessant ? se hvor langt planeten reiste i l?pet av disse to intervallene, det kan gi oss informasjon om hvor i banen vi er. Dette kan vi regne ut ved ? bruke definisjonen av radianer, \(\theta = \frac{b}{r}\), der \(b\) er buelengden, \(\theta\) er vinkelen mellom vektoren som peker p? startposisjonen til planeten og der den befinner seg, og \(r\) er avstanden ut til planeten. Tanken her er ? gj?re \(\theta\) infinitesimalt liten, slik at vi har en \(d\theta\). Vi f?r dermed at \(db = rd\theta\) (Figur 3), ogs? summerer vi opp \(db \) over hele intervallet slik at vi til slutt finner den totale lengden \(b\) planeten reiste. \(d\theta \) har vi ikke, men den kan vi finne ved at \(d\theta = \omega dt\) der \(\omega \) er vinkelfarten til planetbanen om sola, og at \(\omega = \frac{v}{r}\) der \(v\) er farten tangensialt p? banen, og \(r \) er avstanden ut. 

Med denne informasjonen regner vi ut hvor langt planeten reiste i hvert tidsintervall:

Her ser vi at planeten reiser lenger i intervallet \(T2\) enn i intervallet \(T1\) (omtrent 1.8millioner km). Siden tidsintervallene er like store, betyr det at gjennomsnittsfarten i intervallet i \(T2\) m? v?re h?yere enn snittfarten i \(T1\), og derfor b?r omr?det ved \(T2\) v?re n?r planetens aphelium. Da blir \(T1\) intervallet som ligger n?r perihelium. Hang du med p? det? Dette er rett og slett \(s = v\cdot \Delta t\), siden \(\Delta t \) er likt for begge intervallene m? alts? snittfarten v?re h?yere i \(T2\) enn i \(T1\). Siden planeten akselererer fortere n?r den n?rmer seg sola, vil det v?re der den har h?yest fart i banen, og dermed tilh?rer \(T2 \) intervallet n?rmest aphelium. Vi sjekker om dette stemmer med dataene ved ? regne ut snittfarten i intervallene (vi summerer opp farten i hvert punkt og deler p? antall m?lepunkter): 

Det viser seg ogs? her som postulert tidligere at snittfarten i intervallet \(T2\) er st?rre enn i \(T1\). Til sammenlikning har jorda en gjennomsnittlig fart i banen sin p? ca. 107208km/h. Til slutt sjekker vi Keplers 3.lov og Newtons forbedret versjon.  Om du ikke husker hvordan de ser ut, s? kan jeg minne deg om:

\(\text{Kepler:} \;\;\;\;\; \text{Newton:} \\ P^2 \propto a^3 \;\;\;\;\; P^2 = \frac{4\pi^2a^3}{G(m_1 + m_2)}\). Vi f?r denne tabellen for alle planetene i solsystemet:

Her ser vi at Newtons og Keplers lov stemmer godt overens. Det ser kanskje ikke s?nn ut ved f?rste ?yekast, men siden konstanten \(G=4\pi^2AU^3yr^{-2}M\odot^{-1}\) s? vil \(4\pi^2\) i Newtons likning forsvinne. Neglisjerer vi ogs? planetens masse (siden den er liten ift. stjerna), og stjernas masse er ca. \(1M\odot \), s? blir disse like. Her er stjerna i solsystemet omtrent  p? \(1.955M\odot \), s? vi vil f? \(P^2 = \frac{a^3}{1.955}\). Siden jeg har brukt \(P\) i tabellen i stedet for \(P^2\) kan vi for eksempel regne ut \(\frac{1}{0.715^2}\approx 1.956\), som er solmassen til sola i systemet v?rt! Siden alle planetmassene er s? sm? ift. stjernas masse blir den omtrentlig lik for alle planetene. Dermed funker Keplers tredje lov slik Kepler beskrev det, \(P^2 = a^3 \), veldig bra i et system der stjerna har masse \(1M\odot\) (MERK at om vi hadde definert masser utifra massen til stjernen v?r istedet for sola s? hadde Keplers lov holdt for \(P^2=a^3\) i stedet for \(P^2\propto a^3\)).

Litt mer diskusjon av resultater

1) Keplers 2. lov

Vi s? tidligere i innlegget at resultatene ble gode, men ikke prikk like. Hvorfor ikke det? Jeg tror avviket i all hovedsak kommer av metoden vi har brukt, og valgt tidssteg i simuleringen. Her har vi brukt approksimasjoner som \(dA = \frac{1}{2}r^2d\theta\), der \(d\theta = \omega dt\), \(\omega = \frac{v}{r}\). Dermed kommer avviket av \(dt\) vi velger, og feilen i farten \(v\) og posisjonen \(r\) fra leapfrog-metoden. Avviket blir lite nok til at det p?virker resultatene, som jeg kommenterer nederst.

2) Lengde tilbakelagt og snittfart

P? samme m?te som forrige punkt tror jeg avvikene her kommer av metode og tidssteg. I tillegg til approksimasjonen i leapfrog-metoden bruker vi her approksimasjonen \(d\theta = \omega dt\) likt som f?r, som er bestemt av det valgte tiddsteget vi bruker, og fart og posisjon \(v,\; r\) fra den numeriske integrasjonen for \(\omega = \frac{v}{r}\).

Tidssteget vi brukte 

S? lenge vi regner numerisk blir det mye tiln?rminger, og derfor er det viktig med sm? nok intervaller slik at resultatet blir godt nok. Samtidig m? vi balansere det meg regnekraft, slik at beregningstiden ikke blir for stor. Vi valgte ? bruke 11 000 tidssteg pr. ?r vi skulle simulere, som betyr at vi hadde en \(dt\) p? ca. \(\frac{1}{11000} \approx 9\cdot 10^{-5}\). Tester vi simulering med 20 000 steg og 50 000 steg tar plutselig beregningene lenger tid uten ? bli noe s?rlig mer n?yaktig enn for simulering med 11 000 steg, og dermed bestemte vi oss for at 11 000 steg var godt nok i de beregningene vi skulle gj?re. Tidssteget er spesielt viktig fordi den inng?r i n?yaktigheten av leapfrog, og 2 ganger i \(d\theta = \omega dt \) (fordi \(\omega \) bruker \(v,\;r\) som ogs? har n?yaktighet fra metoden og \(dt\)).

Forrige innlegg <<                                                                             Neste innlegg >>