Perspektiv er alt
Vi skal for det meste gj?re det veldig likt som i en tidligere post, Solsystemet p? dataform, og bruker samme antakelser og metode for ? g? l?s p? problemet.
Det er dog én vesentlig forskjell: vi m? skifte koordinatsystem! Grunnen til dette er fordi at i et koordinatsystem der stjerna er sentrum i origo, vil vi ikke registrere noen bevegelse for stjernen. Setter du deg p? stjernen og ser p? de andre planetene, ser det ut som at kun er de som beveger seg, og at stjerna vi sitter p? er helt i ro. For ? registrere bevegelsen til stjernen m? vi se p? bevegelsen til begge legemene i et annet system. Massesenteret mellom de to legemene er et veldig bra sted ? sette origo i det nye koordinatsystemet v?rt, fordi det er et inertialsystem. Hva vil det si? Et inertialsystem er et system der det ikke virker noen ytre krefter p? systemet, dvs. \(\Sigma \vec{F} = M\vec{A} = \vec{0}\) (der \(M\) er total masse for systemet, og \(\vec{A}\) er akselerasjonen til massesenteret). Det gj?r at vi ikke trenger ? tenke p? noen ekstra akselerasjoner og Newtons lover er oppfylt, som gj?r det enklere for oss ? regne. Du kan kanskje tenke: "men akselererer ikke hele systemet et sted i universet ogs??". Joda, men her antar vi at ingen ytre krefter virker p? systemet. Alts? neglisjerer vi alle andre krefter p? systemet og antar det kun er stjerna og planeten som eksisterer i dette rommet, og vi kan regne med Newtons lover. Det er veldig vanskelig ? finne et ekte inertialsystem i universet (igjen, livet er ikke lett for fysikere), fordi omtrent alle systemer beveger seg. For eksempel regner vi med Newtons lover p? jorda, uten ? tenke p? at jorda roterer. Men meteorologer kan ikke se bort i fra jordas rotasjon n?r de skal modellere v?ret, for da spiller Corioliskraften en rolle. Alts? er moralen her at st?rrelsen p? systemet vi modellerer spiller en rolle. Dersom vi plasserer en ball i l?se luften 10m overbakken og slipper den vil det ogs? virke en corioliskraft, men siden den er s? liten ser vi bort ifra den. Men dersom du slipper samme ball fra atmosf?ren blir den plutselig mye viktigere. Samme prinsipp gjelder ogs? andre systemer, som systemet vi skal ta en titt p? n?.
Hva blir forskjellen?
Forh?pentligvis har du lest Solsystemet p? dataform, fordi forskjellen blir ikke s?rlig stor. Metoden som blir brukt er kliss lik, s? jeg kommer ikke til ? ta en like n?ye gjennomgang (fordi det er helt samme metode). Idéen er det samme som f?r:
- Vi finner ut hvor lang tid vi vil simulere med Keplers 3. lov, \(P^2 = \frac{4\pi^2a^3}{G(m_1 + m_2)}\).
- Vi finner et uttrykk for akselerasjonen med Newtons 2. lov, \(\Sigma\vec{F} = m\vec{a}\) for begge legemene vi skal se p?.
- Vi setter initialbetingelser for hastighet og posisjon for begge legemene i det nye koordinatsystemet v?rt og bruker leapfrog til ? integrere opp hastighet og posisjon for begge legemene over hele tiden.
S? hva er forskjellen vi m? gj?re? Forskjellen er at vi m? flytte p? initialbetingelsene slik at de er i massesenterkoordinatsystemet \(CM\) og sim-salabim s? er problemet l?st. Det skal jeg vise dere hvordan man gj?r:
Jeg har enda ikke sagt hvordan vi finner massesenteret, men det er heldigvis ikke noe hokus-pokus. Med koordinatsystemet som har stjerna i origo vil posisjonen for massesenteret v?re: \(\vec{R} = \frac{m_1\vec{r}_1 + m_2\vec{r}_2}{m_1m_2}\), hvor \(m_1,\;m_2\) er massene til de to legemene, og \(\vec{r}_1,\;\vec{r}_2\) er posisjonsvektorene. Dette er et passende sted ? innf?re to st?rrelser som dukker opp mye som vi kommer til ? bruke, nemlig:
\(M=m_1 + m_2\) og \(\hat{\mu} = \frac{m_1m_2}{m_1 + m_2}\). \(M\) er alts? totalmassen av systemet, og \(\hat{\mu}\) er det vi kaller for den reduserte massen til systemet. Bit deg merke i disse. Det gj?r at vi kan skrive opp posisjonen for massesenteret med origo i stjerne som:
\(\vec{R} = \frac{\hat{\mu}}{m_2}\vec{r}_1 + \frac{\hat{\mu}}{m_1}\vec{r}_2 = \frac{\hat{\mu}}{m_1}\vec{r}_2\) (fordi \(\vec{r}_1\) er posisjonsvektoren til stjerna ut fra origo, og stjerna ligger i origo). Deriverer vi med hensyn p? tid kan vi f? hastigheten til massesenteret, \(\vec{V} = \frac{\hat{\mu}}{m_2}\vec{v}_1 + \frac{\hat{\mu}}{m_1}\vec{v}_2 = \frac{\hat{\mu}}{m_1}\vec{v}_2\)
Alle vektorer som g?r ut fra massesenterkoordinatsystemet vil ha en \(cm\) i toppen, som f.eks \(\vec{r}_1^{cm}\), mens \(\vec{r}_1\) betyr at den g?r ut fra det andre vilk?rlige koordinatsystemet. Det betyr at stjernas posisjon i massesenterkoordinatsystemet m? v?re \(\vec{r}_1^{cm} = -\vec{R} = -\frac{\hat{\mu}}{m_1}\vec{r}_2\) (skj?nner du hvorfor?). Jeg skal vise deg helt generelt hvordan vi oversetter dataen fra et koordinatsystem til et annet, og det er et par ting vi b?r legge merke til som kan hjelpe oss senere:
- \(\vec{R}^{cm} = \frac{\hat{\mu}}{m_2}\vec{r}_1^{cm} + \frac{\hat{\mu}}{m_1}\vec{r}_2^{cm} = \vec{0}\) (posisjonen til massesenteret fra massesenterkoordinatsystemet)
- \(\vec{r} = \vec{r}_2 - \vec{r}_1 = \vec{r}_2^{cm} - \vec{r}_1^{cm} \) (1)
MERK at med likning (1) s? har vi definert \(\vec{r}\) som vektoren fra legeme 1 (som kommer til ? v?re stjerna) til legeme 2 (som kommer til ? v?re planeten). Kanskje det kan hjelpe med en tegning for ? se litt n?yere p? det her:
Det kan hende det er litt vanskelig ? se helt p? figuren, men her har vi alts? to koordinatsystemer: et vilk?rlig et, og et med massesenteret \(CM\) i origo. Gjennom innlegget kommer alle ting som har indeks 1, som f.eks \(m_1,\;\vec{r}_1\) etc. h?re til stjernen, og alt som indekseres med 2, \(m_2, \vec{r}_2\) etc. h?re til planeten. Vi lar \(\vec{r}\) v?re avstandsvektoren som g?r fra sola til planeten. Denne trenger vi for ? regne ut akselerasjonen med Newtons 2. lov. For ? finne posisjonene til planetene i \(CM\)-systemet kan vi se p? tegningen at:
- \(\vec{r}_1^{cm} = \vec{r}_1 - \vec{R}\) og \(\vec{r}_2^{cm} = \vec{r}_2 - \vec{R}\) (2)
Vi kan finne hastighetene ved ? derivere med hensyn p? tid:
- \(\vec{v}_1^{cm} = \vec{v}_1 - \vec{V}\) og \(\vec{v}_2^{cm} = \vec{v}_2 - \vec{V}\) (3)
Siden vi allerede har initialverdiene \(\vec{r}_1,\;\vec{r}_2,\;\vec{v}_1,\;\vec{v}_2\) og massene \(m_1,\;m_2\) fra dataen i tabellen fra bloggpostene Solsystemet v?rt og Analytiske planetbaner finner vi n? lett initialverdiene for \(\vec{r}_1^{cm},\;\vec{r}_2^{cm},\;\vec{v}_1^{cm},\;\vec{v}_2^{cm}\) med likningene i (1), (2) og (3), som vi skal kaste inn i leapfrog. Da er vi faktisk klare til ? l?se dette problemet. Vi m? f?rst bare bestemme et par ting til: hvilken planet skal vi simulere, og hvor lenge skal vi simulere? Vi vil helst ha en planet med stor masse og n?r sola, slik at vi f?r st?rst mulig bevegelse p? sola. Sl?r vi opp i tabellene ser vi at planet nummer 7 (det som er oppf?rt som planet nr. 8 i tabellen i linkene ovenfor) har mye st?rre masse enn de andre planetene, og ligger relativt n?rme sola (se plott for planetbaner p? tidligere post, s? ser du planetbanen til planet nr. 7). Dermed er dette en veldig god kandidat. Vi velger ogs? ? simulere dette for 1 rotasjon for planet nr. 7 om stjernen, siden det er nok for ? registrere bevegelsen til stjernen. For ? n? l?se problemet, gj?r vi som jeg skrev tidligere i posten:
- Vi bruker Keplers 3. lov for ? finne tiden det tar planet 7 for 1 oml?p om stjerna: \(P = \sqrt{\frac{4\pi^2a^3}{G(m_1 + m_2)}}\approx 7.45yr \)
- Vi bruker n? Newtons 2. lov p? begge legemene og finner akselerasjonen. Fra N2L har vi: \(\vec{F}_1 = m_1\vec{a}_1 = G\frac{m_1m_2}{r^2}\frac{\vec{r}}{r}\) og \(\vec{F}_2 = m_2\vec{a}_2 = -G\frac{m_1m_2}{r^2}\frac{\vec{r}}{r}\). Dette gir oss: \(\vec{a}_1 = G\frac{m_2}{r^2}\frac{\vec{r}}{r}\) og \( \vec{a}_2 = G\frac{m_1}{r^2}\frac{\vec{r}}{r}\), der \(r=|\vec{r}|\). Fra likning (1) har vi at \(\vec{r} = \vec{r}_2^{cm} - \vec{r}_1^{cm}\), som vi kommer til ? bruke i l?kken for leapfrog n?r vi trenger \(\vec{r}\) i akselerasjonen. Husk at alt som er indeksert med 1 er for stjerna, og 2 for planet nr. 7.
- Vi bruker likningene fra (2) og (3) til ? finne initialbetingelsene, og setter dem i leapfrog. Jeg minner om at jeg ikke g?r i detalj p? leapfrog, siden det er gjort i tidligere bloggpost. Vi kj?rer med samme antall tiddsteg og \(dt\) i denne simuleringen ogs?, siden den er god nok og regneeffektiv.
P? tide ? kj?re kode!
Solas vals med planet 7
Etter ? ha kj?rt kode, plotter vi dataen av banene vi f?r og ender opp med dette resultatet:
Vi har f?tt bevegelse for stjerna, og det liker vi. Banen til stjerna er s? liten i forhold til banen for planet 7, s? i plottet til h?yre er kun banen for stjerna plottet slik at vi kan se bevegelsen for stjerna litt tydeligere ogs?. Vi ser at stjerna beveger seg i en bane som strekker seg ut omtrent 0.01AU, som er omtrentlig 1.5 million km om massesenteret i systemet! For ? teste n?yaktigheten av resultatene kunne vi brukt Keplers 2. lov en gang til, men det kan fort bli kjedelig. I stedet utleder jeg en likning for energien i systemet sett fra massesenterkoordinatsystemet, og bruker den istedet. Vi vet at den totale energien i et system kan skrives om:
\(E_{tot} = E_{tot,k} + E_{tot,p}\). Det betyr at vi m? finne uttrykk for den totale kinetiske energien, og den totale potensielle energien. Vi kan starte med f?rstnevnte:
\(E_{tot,k} = E_{1,k} + E_{2,k} = \frac{1}{2}m_1v_{1,cm}^2 + \frac{1}{2}m_2v_{2,cm}^2\).
Vi trenger uttrykk for farten \(v_{1,cm}\) og \(v_{2,cm}\). Husker du tidligere, s? fant vi at:
- \(\vec{R}^{cm} = \frac{\hat{\mu}}{m_2}\vec{r}_1^{cm} + \frac{\hat{\mu}}{m_1}\vec{r}_2^{cm} = \vec{0}\)
- \(\vec{r} = \vec{r}_2 - \vec{r}_1 = \vec{r}_2^{cm} - \vec{r}_1^{cm} \)
L?ser vi dette likningsystemet f?r vi at:
\(\vec{r}_1^{cm} = -\frac{\hat{\mu}}{m_1}\vec{r}\) og \(\vec{r}_2^{cm} = \frac{\hat{\mu}}{m_2}\vec{r}\). Deriverer vi disse uttrykkene med hensyn p? tid f?r vi:
\(\vec{v}_1^{cm} = -\frac{\hat{\mu}}{m_1}\vec{v}\) og \(\vec{v}_2^{cm} = \frac{\hat{\mu}}{m_2}\vec{v}\), som gir oss at \(v_{1,cm}^2 = |\vec{v}_1^{cm}|^2 = \left(\frac{\hat{\mu}}{m_1}v\right)^2\) og \(v_{2,cm}^2 = |\vec{v}_2^{cm}|^2 = \left(\frac{\hat{\mu}}{m_2}v\right)^2\). Setter vi dette inn i likningen for total kinetisk energi ovenfor og gj?r litt algebra ender vi opp med at \(E_{tot,k} = \frac{1}{2}\hat{\mu}v^2\). S? er det potensiell energi:
Total potensiell energi for et system er gitt som \(E_{tot,p} = -G\frac{m_1m_2}{r}\). Det som er viktig ? tenke p? her er at potensiell energi alltid er definert ut ifra et referanse-punkt. Det betyr at vi ikke kan splitte opp potensiell energi i et uttrykk for hver av planetene, siden potensiell energi kun gir mening dersom vi ser p? hele systemet samlet. Dermed er uttrykket for total potensiell energi ovenfor den totale potensiale energien i v?rt system ogs?. Vi kan gj?re en liten omskriving for \(m_1m_2 = \hat{\mu}(m_1 + m_2) = \hat{\mu}M\), og total potensiell energi blir \(E_{tot,p} = -\frac{GM\hat{\mu}}{r}\). Da blir den totale energien i systemet:
\(E_{tot} = \frac{1}{2}\hat{\mu}v^2 - \frac{GM\hat{\mu}}{r}\). Dette uttrykket bruker vi for ? sjekke den totale energien i systemet v?rt. Vi lager to nye arrays i leapfrog, og regner ut kinetisk, potensiell og total energi i hvert tidspunkt med integrasjonsl?kken, og f?r plottet energien:
Plottet ser ganske lovende ut. Den totale energien er konstant, som er godt ? se siden den skal v?re bevart. Kinetisk energi svinger opp n?r potensiell energi svinger ned, og vice versa. Det som kan v?re interessant ? se p? er at den totale energien i systemet er negativ. Jeg kommer ikke til ? utlede dette, men man kan vise at dersom den totale energien i et slikt system er negativt, g?r planetene i ellipsebaner! (og dersom total energi \(E_{tot} = 0\) har vi sirkelbaner). Som en bonus skal jeg ogs? vise hvordan man kan skrive det totale angul?re momentet til systemet. Angul?rt moment er definert som \(\vec{P} = \vec{r}\times\vec{p}\) der \(\vec{r}\) er posisjonsvektoren til legemet vi ser p? og \(\vec{p}\) er bevegelsesmengden. Vi regner totalt angul?rt moment av systemet sett fra massesenteret:
Vi kan skrive det totale angul?re momentet i systemet ved hjelp av superposisjonsprinsippet: \(\vec{P}_{tot} = \vec{P}_1 + \vec{P}_2\). Vi kan videre fomulere uttrykkene som:
\(\vec{P}_1 = \vec{r}_1^{cm}\times m_1\vec{v}_1^{cm}\) og \(\vec{P}_2 = \vec{r}_2^{cm}\times m_2\vec{v}_2^{cm}\). Tidligere fant vi at:
\(\vec{r}_1^{cm} = -\frac{\hat{\mu}}{m_1}\vec{r},\; \vec{r}_2^{cm} = \frac{\hat{\mu}}{m_2}\vec{r},\; \vec{v}_1^{cm} = -\frac{\hat{\mu}}{m_1}\vec{v},\; \vec{v}_2^{cm} = \frac{\hat{\mu}}{m_2}\vec{v}\). Kaster vi disse uttrykkene inn i \(\vec{P}_1\) og \(\vec{P}_2\) og gj?r litt algebra ender vi opp med dette uttrykket for total angul?r moment:
\(\vec{P}_{tot} = \vec{r}\times\hat{\mu}\vec{v}\). Det som er kult ? se er at uttrykkene for total energi og angul?r moment sett fra massesenteret ser nesten helt like ut som de gj?r til vanlig, med massen erstattet med redusert masse \(\hat{\mu}\) og total masse \(M\).
Diskusjon av resultater
N?yaktighet av resultater kommer av valgt tidssteg og metode. Ved ? velge flere tiddsteg og mindre \(dt\) kan vi f? mer n?yaktige svar (opp til et visst punkt), men det g?r utover regneeffektiviteten. Dette er allerede kommentert i posten Solsystemet p? dataform, og siden metoden vi bruker her er kliss lik med unntak av skifte i koordinatsystem vil jeg ikke gjenta det. Det ligger alltid en usikkerhet fra numeriske integrasjonsmetoder. Hvis du husker fra tidligere er det blitt nevnt at leapfrog er en energibevarende metode, og som vi s? var energien i systemet bevart. Eller var den det? La oss se p? et plott der jeg zoomer helt inn p? kurven for total energi:
Hva er det som foreg?r her? Er ikke leapfrog energibevarende likevel? Jeg er ingen ekspert p? numerikk (jeg er fysiker, ikke numeriker), men legg merke til st?rrelsesordenen p? y-aksen: merkelig, sant? N?r vi gj?r ting numerisk vil det alltid ligge en feil fra datamaskinen ogs?. Maskinen har begrenset med ressurser, og vil dermed alltid introdusere en slags feil. Her er feilen (uansett hvor den kommer ifra) s? liten at vi kan se bort fra den. I det st?rre bildet er energien bevart i systemet, og metoden er helt i orden. Dette er data jeg stoler p?, og dermed sender vi det av g?rde til forskergruppen som trengte dem. Samtidig mottar vi data fra et liknende system som vi har analysert n?, og ser hva slags spennende resultater vi kan komme frem til!