INOM
EXAMENSARBETE TEKNIK, GRUNDNIVÅ, 15 HP
STOCKHOLM SVERIGE 2021,
Numerisk simulering av solsystemet
AGNES BJÖRNFOT SANDRA FJELKESTAM
KTH
SKOLAN FÖR TEKNIKVETENSKAP
Sammanfattning
Den h¨ar rapporten unders¨oker olika numeriska metoder och hur v¨al de kan simulera solsystemet.
Tv˚a modeller simuleras: ett tv˚akroppsproblem f¨or att ta fram jordens planetbana runt solen samt ett niokroppsproblem f¨or att ta fram alla ˚atta planeters samt solens bana. Tv˚a- och niokropps- problemet ¨ar formulerade som en matematisk modell med hj¨alp av differentialekvationer fr˚an Newtons lagar. Systemet av andra ordningens ordin¨ara differentialekvationer som framkommer
¨ar ett Hamiltonianskt system d¨ar den totala energin i systemet ¨ar konstant. D¨arf¨or unders¨oker rapporten symplektiska metoder som konserverar energin f¨or systemet. Dessa symplektiska me- toder j¨amf¨ors med numeriska standardmetoder. De numeriska metoderna som analyseras ¨ar fram˚at Euler, bak˚at Euler, Runge-Kutta 4, symplektisk Euler, St¨ormer-Verlet och symplektiska metoder av h¨ogre ordning, mer specifikt av ordning fyra och ordning tio. Metoderna analyseras utifr˚an noggrannhet och f¨orm˚aga att konservera energi.
Resultatet fr˚an simuleringarna blir bland annat att l¨osningen av tv˚akroppsproblemet f¨or ett tidsintervall p˚a 180 000 ˚ar ger ett betydligt mindre fel f¨or fj¨arde ordningens symplektiska me- tod j¨amf¨ort med Runge-Kutta 4. Runge-Kutta 4 ger ¨aven ett st¨orre fel ¨an en andra ordningens symplektisk metod (St¨ormer-Verlet ) och n¨astan lika stort fel som f¨orsta ordningens symplek- tiska metod (symplektisk Euler ). Slutsatsen av rapporten ¨ar att numeriska metoder som ¨ar symplektiska ¨ar betydligt mer ¨onskv¨arda vid simuleringen av solsystemet, speciellt under l˚anga tidsintervall. Symplektiska metoder ¨ar dessutom ¨overl¨agsna med avseende p˚a energikonserve- ring.
Abstract
This paper examines different numerical methods and how well they can simulate the solar sy- stem. Two models are simulated: a two-body problem to obtain earth’s orbit around the sun and then a nine-body problem to include all eight planets’ and the sun’s orbits. The two- and nine-body problem are formulated as mathematical models with the help of differential equa- tions originating from Newton’s laws. The obtained system of second order ordinary differential equations is a Hamiltonian system, where the total energy of the system is constant. Therefore this paper examines symplectic methods that conserve the energy in the system, and compare these to standard numerical methods to review if energy conservation is important to obtain a small error in the solution. The numerical methods that are analysed are forward Euler, back- ward Euler, Runge-Kutta 4, symplectic Euler, St¨ormer-Verlet and symplectic methods of higher order, namely of order four and ten. The methods are analysed based on accuracy and ability to conserve energy.
One of the results is that the solution of the two-body problem, with a duration of 180 000 years, has a considerable smaller error when using a symplectic method of fourth order than Runge-Kutta 4. Runge-Kutta 4 also has a larger error than a symplectic method of second order (St¨ormer-Verlet ) and almost the same error as a symplectic method of first order (symplectic Euler ). The conclusion of this paper is that symplectic numerical methods are definitely more eligible for a simulation of the solar system, especially with large time intervals. Furthermore symplectic methods are superior based on ability to conserve energy.
Inneh˚ allsf¨ orteckning
1 Inledning 1
1.1 Bakgrund . . . 1
1.1.1 Newtons lagar och Keplers problem . . . 1
1.1.2 Hamiltonianska system . . . 2
1.2 Matematisk modell och problemformulering . . . 2
1.2.1 Tv˚akroppsproblemet . . . 2
1.2.2 Niokroppsproblemet . . . 3
1.2.3 Problemformulering . . . 3
2 Metod 3 2.1 Numeriska standardmetoder . . . 3
2.1.1 F¨orsta ordningens system . . . 3
2.1.2 Eulers stegmetoder . . . 4
2.1.3 Runge-Kutta 4 . . . 5
2.2 Symplektiska metoder . . . 5
2.2.1 F¨orsta ordningens symplektisk metod . . . 6
2.2.2 Andra ordningens symplektisk metod . . . 6
2.2.3 Generell symplektisk metod f¨or separabla Hamiltonianer . . . 7
2.2.4 Fj¨arde ordningens symplektisk metod . . . 7
2.2.5 Tionde ordningens symplektisk metod . . . 8
2.3 Felanalys . . . 8
2.4 Applicering . . . 9
3 Resultat 10 3.1 Numeriska standardmetoder . . . 10
3.2 Symplektiska metoder . . . 12
3.2.1 Tv˚akroppsproblemet . . . 12
3.2.2 Niokroppsproblemet . . . 14
3.3 Energikonservering . . . 18
4 Analys 20 4.1 Noggrannhetsanalys . . . 20
4.2 Energikonservering . . . 22
5 Diskussion 23 5.1 Slutsatser . . . 24
6 Tack till 26 7 Referenser 27 8 Bilagor 28 8.1 Bilaga 1: Koefficienter f¨or tionde ordningens metod . . . 28
8.2 Bilaga 2: Begynnelsev¨arden p˚a positionen . . . 28
8.3 Bilaga 3: Begynnelsev¨arden p˚a hastigheten . . . 29
8.4 Bilaga 4: Fysikaliska data . . . 29
8.5 Bilaga 5: 3D-plottar av planetbanorna . . . 30
1 Inledning
1.1 Bakgrund
Solsystemet med solen och dess himlakroppar har intresserat m¨anniskor under m˚anga ˚ar. F¨orem˚al i r¨orelse i rymden p˚averkas av gravitation som g¨or att dessa kan h˚allas kvar i omloppsbanor. Det- ta ¨ar anledningen till att solsystemets ˚atta planeter (Merkurius, Venus, Jorden, Mars, Jupiter, Saturnus, Uranus och Neptunus) roterar i planetbanor kring solen. Astronomer, matematiker och fysiker har historiskt sett velat unders¨oka hur planeternas r¨orelse utvecklas ¨over tid. Detta kan g¨oras med hj¨alp av differentialekvationer, som d˚a beskriver dessa fysikaliska fenomen med en matematisk modell. Ett system av f¨orsta ordningens differentialekvationer (ODE) kan h¨arledas fr˚an Newtons andra lag och Newtons gravitationslag och utifr˚an dessa g˚ar det att ber¨akna pla- neternas r¨orelse i solsystemet.
Differentialekvationerna kan l¨osas med hj¨alp av numeriska metoder som tar fram approximativa v¨arden p˚a himlakropparnas position och hastighet med hj¨alp av initialv¨arden. De numeriska metoderna har olika egenskaper, som till exempel implicit, explicit och symplektisk, samt olika noggrannhetsordningar som syftar p˚a hur exakta metoderna kan vara. Metoderna karakt¨ariseras av noggrannheten i varje tidssteg, d˚a den numeriska till¨ampningen g˚ar ut p˚a att l¨osa ekvationerna genom att skapa diskreta stegintervall med en specifik stegl¨angd.
Det ¨ar naturligtvis ¨onskv¨art med en numerisk metod som har h¨og noggrannhet eftersom den ger en exaktare l¨osning, men den h¨ar rapporten unders¨oker ¨aven n¨ar det ¨ar l¨ampligt att anv¨anda numeriska metoder som ¨ar symplektiska. N¨ar en numerisk metod ¨ar symplektisk inneb¨ar det att den ¨ar energikonserverande. Planetbanor ¨ar ett typiskt exempel d¨ar den totala energin alltid ¨ar konstant. Numeriska metoder som inte ¨ar symplektiska kommer inte att bevara energin och ger ett globalt fel i energin som v¨axer med tiden.
1.1.1 Newtons lagar och Keplers problem
Figur 1: Kraften fr˚an gravitationen mellan tv˚a kroppar a och b.
Figur 1 visar Newtons gravitationslag som beskriver kraften som tv˚a kroppar uts¨atter varandra f¨or enligt
F = Gmamb
r2 , (1)
d¨ar F ¨ar beloppet av kraften mellan kropp a och b, G ¨ar gravitationskonstanten, ma¨ar massan f¨or kropp a, mb ¨ar massan f¨or kropp b och r ¨ar avst˚andet mellan kropparna. Om ekvation (1) skrivs om till vektorform f˚as d˚a
F = −Gmamb
|rba|2ˆrba= −Gmamb
|rba|3rba, (2)
d¨ar ˆrba=|rrba
ba|¨ar enhetsvektorn och rba=
xb− xa
yb− ya
zb− za
¨ar avst˚andet mellan de tv˚a kropparna och (xi, yi, zi) ¨ar koordinaterna f¨or positionen hos kropp i.
Fr˚an Newtons andra lag f˚as ¨aven f¨or en kropp i r¨orelse att
F = m · a, (3)
d¨ar F ¨ar kraften, m ¨ar massan hos kroppen i r¨orelse och a ¨ar accelerationen hos kroppen. Efter ins¨attning av ekvation (2) i ekvation (3) f˚as f¨oljande:
maa = −Gmamb
|rba|3rba, (4)
d¨ar a|= ¨xa, ¨ya, ¨za. Detta ¨ar ¨aven kallat Keplers problem. L¨osning av ODE:en i ekvation (4) ger kropparnas r¨orelse och hastighet relativt varandra.
1.1.2 Hamiltonianska system
Funktionen H(p, q) ¨ar Hamiltonianen f¨or ODE:er som kallas Hamiltonianska system vilket ¨ar ett system som styrs av Hamiltons ekvationer [5]. Hamiltonianen, H(p, q), representerar den totala energin i systemet som alltid ¨ar konstant n¨ar den inte explicit beror p˚a tiden. I en planetbana
¨ar energin konstant och planeternas banenergi kan beskrivas enligt Hamiltonianen H(p, q) = p2
2 · ma −Gmamb
q , (5)
d¨ar q ¨ar positionen och p ¨ar r¨orelsem¨angden hos planeten. B˚ade ekvation (4) och (5) beskriver en himlakropps r¨orelsebana och allts˚a ¨ar ekvation (5) en Hamiltonian till ODE:en i ekvation (4). Vissa Hamiltonianer ¨ar separabla vilket betyder att de kan skrivas p˚a formen
H(p, q) = H1(p) + H2(q). (6)
H1 beror allts˚a bara p˚a p och H2 beror bara p˚a q. Ekvation (5) ¨ar separabel d¨ar H1= 2mp2
a och H2 = −Gmqamb. Detta inneb¨ar att numeriska metoder anpassade f¨or separabla Hamiltonianer kan till¨ampas f¨or detta problem.
1.2 Matematisk modell och problemformulering
I detta arbete studeras tv˚a olika problem, tv˚akroppsproblemet och niokroppsproblemet. Dessa ska l¨osas med olika numeriska metoder f¨or att f˚a fram hur planeterna r¨or sig i banor i solsyste- met.
1.2.1 Tv˚akroppsproblemet
Tv˚akroppsproblemet handlar om att ta fram hur tv˚a himlakroppar r¨or sig runt ett masscentrum i en elliptisk bana. Fr˚an Newtons lagar framkom ett andra ordningens system av olinj¨ara ODE:er.
F¨or att simplifiera problemet s˚a betraktas nu endast jordens planetbana runt solen, d¨ar solen d˚a blir masscentrum. Solen antas vara stillast˚aende och placerad i origo vilket fr˚an ekvation (4) leder till systemet:
¨ x
¨ y
¨ z
= −G ms
kr3kr, (7)
d¨ar ms¨ar solens massa och r|=x, y, z som ¨ar jordens position relativt solen.
Tv˚akroppsproblemet ¨ar ett Hamiltonianskt system d¨ar energin representeras med Hamiltonia- nen
H(p, q) = kpk2 2mj
−Gmjms
kqk , (8)
d¨ar p ¨ar r¨orelsem¨angden i x-, y- och z-led hos jorden, q ¨ar positionen i x-, y- och z-led hos jorden och mj ¨ar jordens massa.
1.2.2 Niokroppsproblemet
N¨ar fler ¨an tv˚a himlakroppar betraktas m˚aste det ocks˚a tas h¨ansyn till att de olika kropparna p˚averkar varandra. F¨or solsystemet blir det d˚a ett niokroppsproblem f¨or ˚atta planeter och solen.
Detta bildar ett stort system av ODE:er med sex ekvationer f¨or varje himlakropp. Det ¨ar tre ekvationer i position och tre ekvationer i hastighet. Planetbanorna kan l¨osas numeriskt med hj¨alp av ODE:ernas Hamiltonian. En generell Hamiltonian f¨or N + 1 antal kroppar vars banor p˚averkar varandra ¨ar:
H(p, q) = 1 2
N
X
i=0
kpik2 mi
− G
N
X
i=1 i−1
X
j=0
mimj
1
kqi− qjk, (9)
d¨ar pi och qi¨ar r¨orelsem¨angden i x-, y- och z-led f¨or vardera (N + 1)-himlakroppar och mi ¨ar himlakropps i :s massa. F¨or simulering av solsystemet med de ˚atta planeterna och solen s¨atts d˚a N = 8.
1.2.3 Problemformulering
Tv˚akroppsproblemet ska l¨osas f¨or jordens bana runt solen med hj¨alp av olika sorters numeriska metoder f¨or begynnelsev¨ardesproblem. ¨Aven niokroppsproblemet ska l¨osas med hj¨alp av symplek- tiska metoder och systemets Hamiltonian. Metodernas noggrannhet och hur stort fel de ger ska unders¨okas med erk¨anda metoder. D˚a energin f¨or systemet ska vara konstant ska det unders¨okas hur v¨al energin har konserverats efter de numeriska ber¨akningarna. Slutligen ska metoderna j¨amf¨oras utifr˚an hur v¨al de lyckats simulera solsystemet utifr˚an kriterierna ovan.
2 Metod
2.1 Numeriska standardmetoder
2.1.1 F¨orsta ordningens system
Tv˚akroppsproblemet ger ett andra ordningens system av ODE:er enligt ekvation (7). F¨or att im- plementera en l¨osning av problemet med numeriska standardmetoder beh¨ovs ett f¨orsta ordning- ens system. F¨or att g¨ora en omskrivning till f¨orsta ordningens system ans¨atts en vektor:
Q =
q1 q2 q3 q4 q5
q6
=
x y z
˙ x
˙ y
˙ z
. (10)
Sedan deriveras vektorn och eftersom andraderivatorna ¨x, ¨y, och ¨z ¨ar accelerationerna, ¨ar de givna fr˚an ekvation (7). Accelerationerna s¨atts in systemet och ger f¨oljande f¨orsta ordningens system.
Q =˙
˙ x
˙ y
˙ z
¨ x
¨ y
¨ z
=
q4
q5
q6
−Gms q1 (q21+q22+q23)3/2
−Gms q2
(q21+q22+q23)3/2
−Gms q3 (q21+q22+q23)3/2
(11)
Detta f¨orsta ordningens system kan l¨osas med standardmetoder f¨or begynnelsev¨ardesproblem f¨or givna startpositioner x0, y0och z0och starthastigheter ˙x0, ˙y0, och ˙z0.
2.1.2 Eulers stegmetoder
F¨or att l¨osa begynnelsev¨ardesproblemet numeriskt delas tidsaxeln in i diskreta stegintervall och differentialekvationen diskretiseras. Tv˚a standardmetoder ¨ar fram˚at Euler som ¨ar en explicit stegmetod och bak˚at Euler som ¨ar en implicit stegmetod [1]. Explicit inneb¨ar att n¨astkommande iteration ber¨aknas fr˚an nuvarande tidssteg medan implicit tar fram n¨asta tidssteg genom att l¨osa en ekvation som b˚ade inneh˚aller nuvarande och n¨astkommande iteration. Fram˚at Euler till¨ampas genom att derivatan i differentialekvationerna approximeras med fram˚atdifferens. F¨or att implementera l¨osningen av ODE:n s¨atts stegl¨angd h, numerisk l¨osningsvektor Qnum och derivatan fr˚an ekvation (11), ˙Q = f (Q). Index i anv¨ands f¨or aktuellt tidsteg. Fram˚at Euler blir d˚a:
Qnum,i+1= Qnum,i+ h · f (Qnum,i), i = 0, 1, 2, ..., n, (12) d¨ar n ¨ar det totala antalet tidssteg. Bak˚at Euler till¨ampas ist¨allet genom att derivatan approx- imeras med bak˚atdifferens. Implementationen f¨or varje tidssteg blir d˚a:
Qnum,i+1= Qnum,i+ h · f (Qnum,i+1), i = 0, 1, 2, ..., n. (13) D˚a bak˚at Euler ¨ar en implicit metod beh¨over ett olinj¨art ekvationssystem l¨osas i varje steg. F¨or detta anv¨ands Newtons metod f¨or system. Fr˚an ekvation (13) flyttas h¨ogerledet till v¨ansterledet och bildar den vektorv¨arda funktionen g(Qi+1) som anv¨ands f¨or att approximera l¨osningen i det nya tidssteget:
g(Qnum,i+1) = Qnum,i+1− Qnum,i− h · f (Qnum,i+1), i = 0, 1, 2, ..., n. (14) Newtons metod f¨or system ger en numerisk l¨osning vid det nya tidssteget efter att iterationen ger ett fel mindre ¨an 10−15. Indexeringen k + 1 h¨anvisar till roten vid den nuvarande iterationen och indexeringen k ¨ar roten fr˚an en tidigare iteration.
Qnum,i+1,k+1= Qnum,i+1,k−
Jg(Qnum,i+1,k)−1· g(Qnum,i+1,k)
, k = 0, 1, 2, ..., n. (15) Jg(Qi+1) ¨ar jacobianen av ekvation (14) vilket i det till¨ampade fallet med tv˚akroppsproblemet blir:
Jg(Qi+1) = Gms
0 0 0 Gm1
s 0 0
0 0 0 0 Gm1
s 0
0 0 0 0 0 Gm1
s
−2q21+q22+q32 (q12+q22+q23)5/2
−3q1q2
(q21+q22+q23)5/2
−3q1q3
(q21+q22+q32)5/2 0 0 0
−3q1q2
(q12+q22+q23)5/2
q21−2q22+q23 (q21+q22+q23)5/2
−3q2q3
(q21+q22+q32)5/2 0 0 0
−3q1q3 (q12+q22+q23)5/2
−3q2q3 (q21+q22+q23)5/2
q21+q22−2q23
(q21+q22+q32)5/2 0 0 0
. (16)
2.1.3 Runge-Kutta 4
En stegmetod f¨or att l¨osa begynnelsev¨ardesproblem med h¨ogre noggrannhet ¨an Eulers stegme- toder ¨ar fj¨arde ordningens Runge-Kuttametod. Metoden g˚ar ut p˚a att ta ett tidssteg genom att ta ett viktat medelv¨ardet av fyra olika lutningar [1].
k1= f (Qnum,i) k2= f (Qnum,i+h
2 · k1) k3= f (Qnum,i+h
2 · k2) k4= f (Qnum,i+ h · k3) Qnum,i+1= Qnum,i+h
6 · (k1+ 2k2+ 2k3+ k4), i = 0, 1, 2, ..., n.
(17)
Ekvation (17) ger Runge-Kutta 4 med noggrannhetsordning fyra.
2.2 Symplektiska metoder
Symplektiska metoder kan anv¨andas f¨or att l¨osa ODE:er som ¨ar Hamiltonianska system f¨or att konservera energin i l¨osningen. Metoderna som anv¨ands i rapporten utg˚ar fr˚an Hamiltons ekva- tioner som visar att tidsderivatan av r¨orelsem¨angden samt tidsderivatan av positionen motsvarar komponenter av gradienten till Hamiltonianen H(p, q) [5]:
˙
p = −∇qH(p, q)
˙
q = ∇pH(p, q), (18)
d¨ar ˙p ¨ar tidsderivatan av r¨orelsem¨angden, ˙q ¨ar tidsderivatan av position och indexering q och p h¨anvisar till de olika komponenterna av gradienten. Dessa komponenter av gradienten best¨ams enligt:
∇H(p, q) =
"
∂H
∂q (p, q),∂H
∂p(p, q)
#
=
"
∇qH(p, q), ∇pH(p, q)
#
(19) Ekvation (18) anv¨ands f¨or att implementera l¨osningen till ODE:en med ett diskret tidsintervall d¨ar det i varje tidssteg stegas i b˚ade position och r¨orelsem¨angd. F¨or att till¨ampa de symplek- tiska metoderna f¨or tv˚a- och niokroppsproblemet anv¨ands vardera systems Hamiltonian. F¨or tv˚akroppsproblemet g¨aller Hamiltonianen fr˚an ekvation (8). Den ¨ar separabel vilket g¨or att de tv˚a komponenterna av gradienten bara beror p˚a en variabel vardera och dessa blir
∇pH(p) = ∂H(p)
∂p = p
mj ∇qH(q) = ∂H(q)
∂q = Gmjms
kqk3 q (20)
d¨ar p och q ¨ar r¨orelsem¨angd och position i x-, y- och z-koordinater hos jorden.
F¨or niokroppsproblemet g¨aller fortfarande att Hamiltonianen fr˚an ekvation (9) ¨ar separabel.
De partiella derivatorna av Hamiltonianen i gradienten blir dock annorlunda ¨an ekvation (20) eftersom den nu best˚ar av 18 vektorer (position och r¨orelsem¨angd f¨or 9 kroppar). Komponen- terna av gradienten blir d˚a tv˚a vektorv¨arda funktioner enligt:
∇pH(p, q) = ∂H
∂p = pk
mk
∇qH(p, q) = ∂H
∂q = −
k−1
X
i=1 i6=k
Gmkmi kqk− qik3qk+
8
X
i=k+1
Gmimk kqi− qkk3qk,
(21)
d¨ar pk ¨ar r¨orelsem¨angden och qk ¨ar positionen hos himlakropp k, k = 0, 1, 2, ..., 8, i x-, y- och z-koordinater. Massan hos kropp k ges av mk och hos kropp i av mi. Positionen qi ¨ar position hos himlakropp i i x-, y- och z-koordinater. D¨armed f˚as en ny kolumn i gradienten f¨or varje ny himlakropp k enligt:
∇pH(p, q) =p0 m0, mp1
1, mp2
2, mp3
3, mp4
4, mp5
5, mp6
6, mp7
7, mp8
8
∇qH(p, q)|=
GP8 i=1
m0·mi
kq0−qik3q0|
−Gkqm1·m0
1−q0k3q1|+ GP8 i=2
m1·mi
kq1−qik3q1|
−GP1 i=0
m2·mi
kq2−qik3q2|+ GP8 i=3
m2·mi
kq2−qik3q2|
−GP2 i=0
m3·mi
kq3−qik3q3|+ GP8 i=4
m3·mi
kq3−qik3q3|
−GP3 i=0
m4·mi
kq4−qik3q4|+ GP8 i=5
m4·mi
kq4−qik3q4|
−GP4 i=0
m5·mi
kq5−qik3q5|+ GP8 i=6
m5·mi kq5−qik3q5|
−GP5 i=0
m6·mi
kq6−qik3q6|+ GP8 i=7
m6·mi kq6−qik3q6|
−GP6 i=0
m7·mi
kq7−qik3q7|+ G m7·m8
kq7−q8k3q7|
−GP7 i=0
m8·mi
kq8−qik3q8|
(22)
2.2.1 F¨orsta ordningens symplektisk metod
En f¨orsta ordningens symplektisk metod ¨ar symplektisk Euler [3]. Hamiltonianen f¨or systemet anv¨ands f¨or att ber¨akna n¨asta position efter en stegl¨angd h p˚a liknande sett som Eulers steg- metoder fr˚an tidigare kapitel. En semi-implicit version av symplektisk Euler ¨ar:
pi+1= pi− h · ∇qH(pi+1, qi)
qi+1= qn+ h · ∇pH(pi+1, qi), i = 0, 1, 2, ..., n, (23) d¨ar i ¨ar aktuella tidssteget, n ¨ar antalet steg, h ¨ar stegl¨angden, q ¨ar positionen och p ¨ar r¨orelsem¨angden. Metoden ¨ar semi-implicit eftersom steget i r¨orelsem¨angden p ¨ar implicit medan steget i position q ¨ar explicit. Metoden blir dock explicit d˚a Hamiltonianen ¨ar separabel enligt ekvation (6) eftersom ∇qH(p, q) inte kommer bero av p. F¨or system med separabla Hamiltoni- aner blir d˚a symplektisk Euler
pi+1= pi− h · ∇qH(qi)
qi+1= qi+ h · ∇pH(pi+1), i = 0, 1, 2, ..., n. (24) 2.2.2 Andra ordningens symplektisk metod
En andra ordnings symplektisk metod ¨ar till exempel St¨ormer-Verlet [3]. Metoden g˚ar ut p˚a att det tas tv˚a steg i r¨orelsem¨angd och position f¨or varje tidssteg, till skillnad fr˚an symplektisk Euler d¨ar det tas ett steg per tidsteg. Om stegen f¨orst tas i position blir metoden
qi+1/2= qi+h
2 · ∇pH(pi, qi+1/2) pi+1= pi−h
2 ·
∇qH(pi, qi+1/2) + ∇qH(pi+1, qi+1/2) qi+1= qi+1/2+h
2 · ∇pH(pi+1, qi+1/2), i = 0, 1, 2, ..., n.
(25)
Notera att ¨aven denna metod ¨ar semi-implicit f¨or en generell Hamiltonian och explicit f¨or en separabel Hamiltonian. F¨or system med separabla Hamiltonianer blir d˚a St¨ormer-Verlet
qi+1/2= qi+h
2 · ∇pH(pi) pi+1= pi− h · ∇qH(qi+1/2) qi+1= qi+1/2+h
2 · ∇pH(pi+1), i = 0, 1, 2, ..., n.
(26)
2.2.3 Generell symplektisk metod f¨or separabla Hamiltonianer
F¨or att hitta h¨ogre ordningens symplektiska metoder f¨or separabla Hamiltonianer har Yoshida [11] tagit fram en generell s-steg integrator f¨or varje diskret tidssteg:
qk+1= qk+ hck+1∂H
∂p(pk) pk+1= pk− hdk+1
∂H
∂q (qk+1), k = 0, 1, 2, ..., s − 1,
(27)
d¨ar k ¨ar det aktuella steget. Allts˚a blir s totala antalet steg f¨or varje tidssteg. Efter att s steg har tagits tas ett nytt diskret tidssteg f¨or metoden genom att upprepa ekvation (27), vilket g¨ors tills totala antalet tidssteg har n˚atts. Koefficienterna c och d summeras vardera till ett eftersom det totalt kommer att tas ett diskret tidssteg per s steg:
c1+ c2+ ... + cs= 1
d1+ d2+ ... + ds= 1. (28)
Observera att f¨or s = 1 f˚as c1 = 1 och d1 = 1, ger ekvation (27) samma ekvation som f¨or symplektisk Euler f¨or separabla Hamiltonianer enligt ekvation (24). En l¨osning f¨or s = 2 ¨ar d˚a
¨aven St¨ormer Verlet f¨or separabla Hamiltonianer. Denna generella form kan ¨aven anv¨andas f¨or st¨orre v¨arden p˚a s f¨or att skapa symplektiska metoder av h¨ogre ordning.
Figur 2: Tidslinje som visar hur stegen i position och r¨orelsem¨angd f¨orh˚aller sig till varandra.
Figur 2 visar hur den generella s-steg integratorn kan skapa h¨og noggrannhet enligt ekvation (27) genom att varje steg i positionen utnyttjar det f¨oreg˚aende steget i r¨orelsem¨angden och vice versa.
2.2.4 Fj¨arde ordningens symplektisk metod
En fj¨arde ordningens metod, d˚a s = 4, uppt¨acktes av Forest och Ruth 1990 [2]. Koefficienterna f¨or denna metod ¨ar:
c1= c4= 1
2(2 − 21/3), c2= c3= 1 − 21/3 2(2 − 21/3), d1= d3= 1
2 − 21/3, d2= − 21/3
2 − 21/3, d4= 0.
(29)
Metoden implementeras enligt ekvation (27) d¨ar det allts˚a tas fyra steg i position och r¨orelsem¨angd f¨or varje tidssteg h.
2.2.5 Tionde ordningens symplektisk metod
F¨or att hitta h¨ogre ordningens l¨osningar kan koefficienter w0, w1, ..., wr, r =s−22 anv¨andas. Dessa anv¨ands f¨or att ta fram koefficienterna c1, c2..., csoch d1, d2..., dsi ekvation (28). L¨osningen blir d˚a att:
d2r+2= 0,
d1= d2r+1= wr, d2= d2r= wr−1, ..., dr= dr+2= w1, dr+1= w0
c1= c2r+2=1 2wr, c2= c2r+1=1
2(wr+ wr−1), ..., cr+1= cr+2=1
2(w1+ w0).
(30)
En tionde ordningens metod, d˚a s = 34, togs fram av Tsitouras 1999 [10]. Eftersom s = 34 s˚a
¨
ar r = 16. Det blev d˚a 34 stycken koefficienter ck och dk varderda enligt ekvation (30):
d34= 0,
d1= d33= w16, d2= d32= w15, d3= d31= w14, ..., d16= d18= w1, d17= w0
c1= c34= 1 2w16, c2= c33= 1
2(w16+ w15), c3= c32= 1
2(w15+ w14), ..., c17= c18=1
2(w1+ w0).
(31)
De 17 olika koefficenterna w0, w1..., w16togs fram av Tsitouras [10] med 45 decimalers noggrann- het och finns i Bilaga 1.
2.3 Felanalys
F¨or att unders¨oka hur noggranna de olika numeriska metoderna ¨ar g¨ors olika noggrannhetskon- troller f¨or att se att metoderna verkligen har korrekt noggrannhetsordning enligt teorin och ¨aven uppskatta hur stort felet blir. Ett s¨att att uppskatta felet ¨ar att ta slutpositionen med en viss stegl¨angd och sedan upprepa ber¨akningen med en halverad stegl¨angd. Tas normen av differensen mellan dessa kommer det att ges ett approximativ fel enligt:
˜ eh=
˜uh− ˜uh/2
(32)
d¨ar ˜eh¨ar det approximativa felet och ˜uh,h/2¨ar l¨osningen i sista steget f¨or tv˚a olika stegl¨angder h [9]. Om de uppskattade felen plottas mot stegl¨angden i en logaritmisk skala g˚ar det ocks˚a att best¨amma vilken noggrannhet metoden har. Detta g¨ors genom att plotta en referenslinje som har samma lutning som metodens ordning [9].
Ett till experimentellt s¨att f¨or att best¨amma vilken noggrannhetsordning en viss metod har, ¨ar genom en kvot av felen vid olika stegl¨angder enligt:
log2
˜eh
˜ eh/2
= log2
˜uh− ˜uh/2 ˜uh/2− ˜uh/4
= p + O(h), (33)
d¨ar p ¨ar ordningen hos metoden och O(h) ¨ar resten [9]. Det betyder att f¨or en f¨orsta ordning- ens metod avtar differenserna med en faktor tv˚a och f¨or en andra ordningens metod kommer differenserna att avta med en faktor fyra.
D˚a systemet av ODE:er som unders¨oks ¨ar ett Hamiltonianskt system s˚a unders¨oks det ¨aven om metoderna ¨ar symplektiska genom att betrakta konserveringen av energi. Eftersom energin i him- lakropparnas planetbana ska vara konstant i tiden och energin representeras av Hamiltonianen, b¨or dess v¨arde vara konstant i tiden enligt:
H(p0, q0) = H(p(t), q(t)), ∀ t > 0. (34)
Hur v¨al de olika metoderna lyckas med ekvation (34) studeras genom att ber¨akna Hamiltonianen f¨or varje tidssteg och plotta den f¨or att se hur konstant energin h˚aller sig. F¨or tv˚akroppsproblemet
¨ar det d˚a Hamiltonianen fr˚an ekvation (8) som ska plottas och f¨or niokroppsproblemet ¨ar det Hamiltonianen fr˚an ekvation (9) [5].
2.4 Applicering
F¨or att till¨ampa de numeriska metoderna skrivs ber¨akningarna i MATLAB, d¨ar ¨aven alla grafer skapas. Som indata l¨ases planeternas startposition och starthastighet in i x -,y- och z -koordinater via MATLABs inbyggda funktion planetEphemeris [4]. Funktionen planetEphemeris ger indatan av planeterna fr˚an NASA. Position och hastighet hos planterna och solen ¨ar givna relativt ett masscentrum som i tv˚akroppsproblemet ¨ar solen och i niokroppsproblemet ¨ar solsystemets mass- centrum. Startpositionerna och starthastigheterna ¨ar h¨amtade med startdatum fr˚an 1/1-2021, se v¨arden i Bilaga 2 och 3. Fysikaliska data om planeternas massa och gravitationskonstanten som anv¨ands g˚ar att se i Bilaga 4 vilka ¨aven dessa ¨ar h¨amtade fr˚an NASA, gravitationskonstanten fr˚an [6], solens massa fr˚an [8] och planetmassor fr˚an [7].
Alla planetbanor kommer att presenteras i tv˚a dimensioner med x - och y-koordinater f¨or att f˚a tydliga figurer av banorna. Hur planetbanorna ser ut i tre dimensioner g˚ar att se i Bilaga 5.
Allt resultat fr˚an niokroppsproblemet visas separat i de inre och yttre himlakropparna eftersom l¨angden av planetbanorna hos de olika himlakropparna varierar. Merkurius har den kortaste planetbanan p˚a 176 dagar och Neptunus har den l¨angsta planetbanan som tar 165 ˚ar vilket motsvarar 60 190 dagar.
Tv˚akroppsproblemet l¨oses med alla metoder. F¨or att j¨amf¨ora metoderna under en l¨angre tid- period valdes ett tidsintervall p˚a T = 180 000 ˚ar f¨or alla metoder f¨orutom fram˚at Euler, bak˚at Euler och tionde ordningens metod. F¨or det tidsintervallet s¨atts ¨aven stegl¨angden till h = 1 dag eftersom alla metoderna uppvisar sin teoretiska noggrannhet f¨or den stegl¨angden och det gav en tillr¨ackligt noggrann planetbana. F¨or niokroppsproblemet s˚a studeras endast de symplektis- ka metoderna och det problemet l¨oses med ett kortare tidsintervall p˚a T = 200 ˚ar som ungef¨ar motsvarar ett varv p˚a den l¨angsta planetbanan som ¨ar Neptunus. ¨Aven niokroppsproblemet l¨oses med h = 1 dag f¨orutom f¨or tionde ordningens metod som l¨oses med h = 8 dagar eftersom den metoden ¨ar noggrannare och mer ber¨akningstung.
MATLABs standardinst¨allning p˚a 16 decimalers noggrannhet begr¨ansar till en viss del noggrann- heten av utr¨akningarna med metoderna. Problemet l¨oses med att anv¨anda skalning p˚a all fy- sikalisk data med AU ist¨allet f¨or km s˚a att ett mindre antal siffror beh¨ovs f¨or varje utr¨akning.
F¨or niokroppsproblemet fungerar inte detta helt med sm˚a stegl¨angder f¨or de yttre planeterna.
Avst˚anden de r¨or sig ¨ar d˚a f¨or korta vilket g¨or att felet ackumuleras. Avrundningsfelet f¨orst¨or noggrannheten f¨or tionde ordningens metod d˚a utr¨akningarna sker i v¨aldigt sm˚a steg (34 steg per tidssteg) och felet ofta ¨ar mindre ¨an 10−15. Dessutom ¨ar konstanterna fr˚an Tabell 1, som anv¨ands f¨or att formulera metoden, givna med 45 decimalers noggrannhet vilket betyder att ¨aven MATLAB m˚aste klara 45 decimaler. F¨or att kunna till¨ampa tionde ordningens metod i MATLAB m˚aste allts˚a programmets inbygga noggrannhet h¨ojas. Detta l¨oses med MATLABs inbyggda funktioner vpa och digits. Problemet blir d˚a att ber¨akningskostnaderna f¨or programmet blir betydligt h¨ogre och l¨osningen kan inte tas fram i ett lika l˚angt tidsintervall eller med lika liten stegl¨angd som f¨or l¨agre ordningens metoder. Tv˚akroppsproblemet l¨oses d¨arf¨or endast med T = 10 ˚ar med tionde ordningens metod. Med tv˚akroppsproblemet ¨ar det mer intressant att j¨amf¨ora de numeriska stan- dardmetoderna med de symplektiska metoderna, d¨arav anv¨andes endast ett kort tidsintervall f¨or tionde ordningens metod. Dock f¨or niokroppsproblemet l¨ostes tionde ordningens metod med T = 200 ˚ar eftersom att de yttre planeterna m˚aste hinna ˚aka minst ett varv runt solen om de olika symplektiska metoderna ska kunna j¨amf¨oras med varandra. Jordens planetbana ber¨aknas med fram˚at Euler och bak˚at Euler, i endast T = 5 ˚ar j¨amf¨ort med de andra metoderna som f¨or tv˚akroppsproblemet ber¨aknas i T = 180 000 ˚ar. Dessutom k¨ors fram˚at Euler och bak˚at Euler med mindre stegl¨angd j¨amf¨ort med de andra metoderna. Detta ¨ar dels f¨or att metoderna ger ett stort fel redan av korta tidsintervaller och ocks˚a p˚a grund av att bak˚at Euler inte klarar av ett l¨angre tidsintervall, d˚a jordens bana tillslut krockar med solen och det g˚ar inte forts¨atta att l¨osa planetbanan.
3 Resultat
3.1 Numeriska standardmetoder
Med fram˚at och bak˚at Euler har endast tv˚akroppsproblemet studerats. I tiden T = 5 ˚ar och med stegl¨angden h = 1/4 dag f˚as f¨oljande resultat f¨or jordens planetbana runt solen.
(a) Fram˚at Euler (b) Bak˚at Euler
Figur 3: L¨osning av tv˚akroppsproblemet med fram˚at och bak˚at Euler i T = 5 ˚ar och h = 1/4 dag.
Tv˚akroppsproblemet l¨ostes ¨aven med Runge-Kutta 4 och i tiden T = 180 000 ˚ar och med stegl¨angden h = 1 dag f˚as f¨oljande resultat:
Figur 4: L¨osning av tv˚akroppsproblemet med Runge-Kutta 4 i T = 180 000 ˚ar och med h = 1 dag.
Felet f¨or de numeriska standardmetoderna, f¨or l¨osningarna som ses i Figurerna 3 och 4, har uppskattats enligt avsnitt 2.3 och ekvation (32).
Tabell 1: Uppskattade felet f¨or tv˚akroppsproblemet med de numeriska standardmetoderna.
Metod T [˚ar] h [dag] Uppskattade felet [AU]
Fram˚at Euler 5 1/4 1,9258
Bak˚at Euler 5 1/4 1,2207
Runge-Kutta 4 180 000 1 1,5693
Den uppskattade noggrannheterna fr˚an numeriska metoderna vid en viss tid och stegl¨angder har ber¨aknats enligt ekvation (33) och j¨amf¨ors med den teoretiska noggrannheten f¨or att se att de st¨ammer ¨overens.
Tabell 2: Ber¨aknad noggrannhet f¨or de numeriska standardmetoderna.
Metod T [dag] h [dag] Uppskattad
noggrannhet
Teoretisk noggrannhet
Fram˚at Euler 64 1, 1/2, 1/4 0,9607 1
Bak˚at Euler 64 1, 1/2, 1/4 1,0444 1
Runge-Kutta 4 256 4, 2, 1 4,1281 4
De logaritmiska plottarna av de uppskattade felen togs fram f¨or att unders¨oka metodernas noggrannhet och se hur felet ¨andrats. Samma tidsintervall som n¨ar noggrannheten ber¨aknades i Tabell 2 anv¨andes f¨or de olika metoderna.
(a) Fram˚at Euler (b) Bak˚at Euler
Figur 5: De uppskattade felen f¨or fram˚at Euler och bak˚at Euler i logaritmisk skala.
Figur 6: De uppskattade felen f¨or Runge-Kutta 4 i logaritmisk skala.
3.2 Symplektiska metoder
Med de symplektiska metoderna har b˚ade tv˚akroppsproblemet och niokroppsproblemet analyse- rats, med solen och de ˚atta planeterna.
3.2.1 Tv˚akroppsproblemet
Med symplektiskt Euler och St¨ormer-Verlet i tiden T = 180 000 ˚ar och stegl¨angden h = 1 dag blev det f¨oljande resultat f¨or solen och jorden:
(a) Symplektisk Euler (b) St¨ormer-Verlet
Figur 7: L¨osning av tv˚akroppsproblemet med symplektisk Euler och St¨ormer-Verlet med h = 1 dag och T = 180 000 ˚ar.
Med fj¨arde ordningens metod i tiden T = 180 000 ˚ar och med stegl¨angden h = 1 dag och tionde ordningens metod i tiden T = 10 ˚ar och med stegl¨angden h = 1 dag blev det f¨oljande resultat f¨or solen och jorden:
(a) Fj¨arde ordningens metod (b) Tionde ordningens metod
Figur 8: L¨osning av tv˚akroppsproblemet med fj¨arde ordningens metod och tionde ordningens metod med h och T enligt ovan.
Felet f¨or de symplektiska metoderna, vid l¨osningarna fr˚an Figurer 7 och 8 har uppskattats enligt avsnittet 2.3 f¨or tv˚akroppsproblemet och ekvation (32).
Tabell 3: Uppskattade felet f¨or tv˚akroppsproblemet med de symplektiska metoderna.
Metod T [˚ar] h [AU] Uppskattade felet [AU]
Symplektisk Euler 180 000 1 1,996788406166508 St¨ormer-Velert 180 000 1 0,829913306713017 Fj¨arde ordningen 180 000 1 0,075983579034666 Tionde ordningen 10 1 1,715507801925319 ·10−18
Noggrannheten testades enligt ekvation (33) f¨or de olika metoderna med olika tidsintervall och stegl¨angder.
Tabell 4: Ber¨aknad noggrannhet f¨or metoderna.
Metod T [dag] h [dag] Uppskattad noggrannhet
Teoretisk noggrannhet
Symplektisk Euler 128 1, 2, 1/2 0,9967 1
St¨ormer-Verlet 256 4, 2, 1 1,9999 2
Fj¨arde ordningen 256 4, 2, 1 3,9993 4
Tionde ordningen 256 4, 2, 1 9,9735 10
De uppskattade felen plottades logaritmiskt f¨or att unders¨oka metodernas noggrannhet och hur felet ¨andrades. Samma tidsintervall som i Tabell 4 anv¨andes f¨or de olika metoderna.
(a) Symplektiskt Euler (b) St¨ormer-Verlet
Figur 9: De uppskattade felen f¨or f¨orsta och andra ordningens metod i logaritmisk skala.
(a) Fj¨arde ordningen (b) Tionde ordningen
Figur 10: De uppskattade felen f¨or fj¨arde och tionde ordningens metod i logaritmisk skala.
3.2.2 Niokroppsproblemet
Niokroppsproblemet har l¨osts med de symplektiska metoderna i T = 200 ˚ar och med stegl¨angden h = 1 dag f¨or symplektisk Euler, St¨ormer-Verlet och fj¨arde ordningens metod och med h = 8 dagar f¨or tionde ordningens metod.
(a) Hela solsystemet (b) De fem inre himlakropparna
Figur 11: L¨osningen av niokroppsproblemet med symplektisk Euler i T = 200 ˚ar och h = 1 dag.
(a) Hela solsystemet (b) De fem inre himlakropparna Figur 12: L¨osningen av niokroppsproblemet med St¨ormer-Verlet i T = 200 ˚ar och h = 1 dag.
(a) Hela solsystemet (b) De fem inre himlakropparna
Figur 13: L¨osningen av niokroppsproblemet med fj¨arde ordningens metod i T = 200 ˚ar och h = 1 dag.
(a) Hela solsystemet (b) De fem inre himlakropparna
Figur 14: L¨osningen av niokroppsproblemet med tionde ordningens metod i T = 200 ˚ar och h = 8 dagar.
De uppskattade felen f¨or de symplektiska metoderna vid l¨osningen av planetbanorna med tiden T = 200 ˚ar och h = 1 dag f¨or alla metoder utom tionde ordningens metod d¨ar T = 200 ˚ar och h = 8 dagar, f˚as fr˚an de numeriska l¨osningarna som visas i Figurerna 11, 12, 13 och 14 till:
Tabell 5: Uppskattade felet f¨or de inre himlakropparna.
Metod
Planet
Solen Merkurius Venus Jorden Mars
Symplektisk Euler 2,6223·10−5 0,1848 0,3947 0,2200 0,5791 St¨ormer-Verlet 1,0312·10−6 0,0612 0,2858 0,0943 0,0232 Fj¨arde ordningen 1,1211·10−8 0,0647 6,8324·10−4 8,5089·10−5 6,8567·10−6 Tionde ordningen 1,8864·10−9 0,0114 5,5235 ·10−7 7.6357 ·10−8 2,4768 ·10−8
Tabell 6: Uppskattade felet f¨or de yttre himlakropparna.
Metod
Felet [AU]
Jupiter Saturnus Uranus Neptunus
Symplektisk Euler 0,0267 0,0028 0,0020 0,0014
St¨ormer-Verlet 3,0278·10−4 3,7273·10−5 2,7996·10−6 5,5630·10−7 Fj¨arde ordningen 2,1851·10−9 1,9567·10−10 3,3607·10−11 3,8061·10−11 Tionde ordningen 1,4223·10−10 3,0856 ·10−11 5,4271·10−12 6,3528·10−12
Noggrannhetsordningen ber¨aknades med hj¨alp av ekvation (33) f¨or de inre planeterna med stegl¨angden h = [1, 1/2, 1/4] dag och tiden T = 256 dagar f¨or symplektiska Euler samt St¨ormer- Verlet Med fj¨arde ordningens metod anv¨andes stegl¨angderna h = [2, 1, 1/2] och tiden T = 8192 dagar. Med tionde ordningens metod anv¨andes stegl¨angderna h = [4, 2, 1] och tiden T = 256 dagar.
Tabell 7: Ber¨aknad noggrannhet f¨or metoderna f¨or de inre planeterna.
Metod
Planet
Solen Merkurius Venus Jorden Mars Teoretisk noggrannhet Symplektisk Euler 1,0016 1,2003 0,8132 0,9444 0,9976 1 St¨ormer-Verlet 1,9948 1,9935 1,9993 1,9997 1,9999 2 Fj¨arde ordningen 3,9817 3,9803 3,9981 3,9993 3,9998 4 Tionde ordningen 9,7673 9,9408 9,9295 9,7574 9,0697 10
Samma sak gjordes f¨or de yttre planeterna med stegl¨angderna h = [4, 2, 1] dag och tiden T = 8192 dagar f¨or symplektisk Euler samt St¨ormer-Verlet Samma tid men stegl¨angderna h = [2, 1, 1/2] dag anv¨andes f¨or fj¨arde ordningens metod. F¨or tionde ordningens metod anv¨andes stegl¨angderna h = [16, 8, 4] dag och tiden T = 256 dagar.
Tabell 8: Ber¨aknad noggrannhet f¨or metoderna f¨or de yttre planeterna.
Metod
Planet
Jupiter Saturnus Uranus Neptunus Teoretisk noggrannhet Symplektisk Euler 1,1247 0,9812 0,9998 1,0000 1 St¨ormer-Verlet 1,9999 2,0000 2,0000 2,0000 2 Fj¨arde ordningen 4,0010 3,9779 3,9899 3,9977 4 Tionde ordningen 9,0797 6,6401 6,0800 6,7422 10
Felen och en referenslinje med f¨orv¨antad noggrannhet plottades f¨or olika stegl¨angder och ti- der.
(a) Inre planeterna med T = 256 dagar. (b) Yttre planeterna med T = 8192 dagar.
Figur 15: Noggrannheten f¨or alla planetbanorna f¨or niokroppsproblemet med symplektisk Euler.
(a) Inre planeterna med T = 256 dagar. (b) Yttre planeterna med T = 8192 dagar.
Figur 16: Noggrannheten f¨or alla planetbanorna f¨or niokroppsproblemet med St¨ormer-Verlet.
(a) Inre planeterna med T = 8192 dagar. (b) Yttre planeterna med T = 8192 dagar.
Figur 17: Noggrannheten f¨or alla planetbanorna f¨or niokroppsproblemet med fj¨arde ordningens metod.
(a) Inre planeterna med T = 256 dagar. (b) Yttre planeterna med T = 256 dagar.
Figur 18: Noggrannheten f¨or alla planetbanorna f¨or niokroppsproblemet med tionde ordningens metod.
I figurerna ovan ¨ar ˜e0, ˜e1, ˜e2, ˜e3, ˜e4, ˜e5, ˜e6, ˜e7 och ˜e8 felet f¨or varje himlakropp, allts˚a felet f¨or Solen, Merkurius, Venus, Jorden, Mars, Jupiter, Saturnus, Uranus respektive Neptunus.
3.3 Energikonservering
Resultatet om hur energin ¨andrades f¨or alla numeriska metoderna plottades och ber¨aknades med hj¨alp av ekvation (8) f¨or tv˚akroppsproblemet.
(a) Fram˚at Euler (b) Bak˚at Euler
Figur 19: V¨ardet av Hamiltonianen f¨or fram˚at Euler och bak˚at Euler med h = 1/2 dag och i T
= 5 ˚ar.
(a) Runge-Kutta 4 (b) Sympletisk Euler
Figur 20: V¨ardet av Hamiltonianen var 100 000:e dag f¨or Runge-Kutta 4 och symplektisk Euler med h = 1 dag och i T = 180 000 ˚ar.
(a) St¨ormer-Verlet (b) Fj¨arde ordningens metod
Figur 21: V¨ardet av Hamiltonianen var 100 000:e dag f¨or St¨ormer-Verlet och fj¨arde ordningens metod med h = 1 dag och i T = 180 000 ˚ar.
Hur energin f¨or¨andrades ¨over tiden T = 10 ˚ar f¨or fj¨arde respektive tionde ordningens symplek- tiska metod j¨amf¨ordes f¨or att se om de skiljde sig ˚at.
(a) Fj¨arde ordningens metod (b) Tionde ordningens metod
Figur 22: V¨ardet av Hamiltonianen f¨or fj¨arde ordningens metod och tionde ordningens metod med h = 1 dag och i T = 10 ˚ar.
F¨or niokroppsproblemet togs energin fram och plottades med hj¨alp av Hamiltonianen fr˚an ek- vation (9) f¨or de symplektiska metoderna.
(a) Symplektisk Euler (b) St¨ormer-Verlet
Figur 23: V¨ardet av Hamiltonianen f¨or symplektisk Euler och St¨ormer-Verlet med h = 1 dag och i T = 200 ˚ar.
(a) Fj¨arde ordningens metod (b) Tionde ordningens metod
Figur 24: V¨ardet av Hamiltonianen f¨or fj¨arde ordningens metod och tionde ordningens metod med h = 1 dag respektive h = 8 dagar och i T = 200 ˚ar.
F¨or att kunna j¨amf¨ora fj¨arde respektive tionde ordningens metod har v¨ardet av Hamiltonoianen tagits fram f¨or tiden T = 10 ˚ar.
(a) Fj¨arde ordningens metod (b) Tionde ordningens metod
Figur 25: V¨ardet av Hamiltonianen f¨or fj¨arde ordningens metod och tionde ordningens metod med h = 1 dag respektive h = 8 dagar och i T = 10 ˚ar.
4 Analys
4.1 Noggrannhetsanalys
Noggrannheten ber¨aknades f¨or alla numeriska metoder och i Tabeller 2, 4, 7 och 8 g˚ar det att se att alla metoder har noggrannhetsordning enligt teorin, f¨orutom noggrannhetensordningen i de yttre planetbanorna med tionde ordningens metod. Detta eftersom noggrannhetsordningen egentligen beh¨ovs testas i ett l¨angre tidsintervall.
N¨ar noggrannhetsordningen testades f¨or de olika metoderna och planetbanorna anv¨andes olika stegl¨angder och tidsintervall. Detta eftersom vissa metoder ¨ar mer ber¨akningstunga ¨an andra och de yttre planeterna har i regel l¨angre tidsintervall ¨an de inre himlakropparna eftersom att de skulle hinna f¨orflytta sig mer d˚a de har l¨angre omloppsbana och l¨agre hastigehet.
De numeriska metoderna kommer inte att ge r¨att noggrannhetsordning f¨or alla stegl¨angder och f¨or olika tidsintervall bekr¨aftas det nedan vilka stegl¨angder som gav r¨att noggrannhetsordning f¨or de olika numeriska metoderna.
Fram˚at Euler
• F¨or tv˚akroppsproblemet med tiden T = 64 dagar visas i Figur 5a att metoden har r¨att noggrannhetsordning f¨or stegl¨angder mellan 1/32 till 4 dagar.
Bak˚at Euler
• F¨or tv˚akroppsproblemet med tiden T = 64 dagar visas i Figur 5b att metoden har r¨att noggrannhetsordning f¨or stegl¨angder mellan 1/32 till 8 dagar.
Runge-Kutta 4
• F¨or tv˚akroppsproblemet med tiden T = 256 dagar visas i Figur 6 att metoden har r¨att noggrannhetsordning f¨or stegl¨angder mellan 1/16 till 64 dagar.
Symplektisk Euler
• F¨or tv˚akroppsproblemet med tiden T = 128 dagar visas i Figur 9a att metoden har r¨att noggrannhetsordning f¨or stegl¨angder mellan 1/8 till 8 dagar.
• F¨or niokroppsproblemet med tiden T = 256 dagar visas i Figur 15a att metoden har r¨att noggrannhet f¨or de inre planetbanorna vid stegl¨angder upp till 2 dagar.
• F¨or niokroppsproblemet med tiden T = 8192 dagar visas i Figur 15b att metoden har r¨att noggrannhet f¨or de yttre planetbanorna vid stegl¨angder upp till 8 dagar.
St¨ormer-Verlet
• F¨or tv˚akroppsproblemet med tiden T = 128 dagar visas i Figur 9b att metoden har r¨att noggrannhetsordning f¨or stegl¨angder mellan 1/8 till 64 dagar.
• F¨or niokroppsproblemet med tiden T = 256 dagar visas i Figur 16a att metoden har r¨att noggrannhet f¨or de inre planetbanorna vid stegl¨angder upp till 8 dagar.
• F¨or niokroppsproblemet med tiden T = 8192 dagar visas i Figur 16b att metoden har r¨att noggrannhet f¨or de yttre planetbanorna vid stegl¨angder upp till 64 dagar.
Fj¨arde ordningens metod
• F¨or tv˚akroppsproblemet med tiden T = 256 dagar visas i Figur 10a att metoden har r¨att noggrannhetsordning f¨or stegl¨angder mellan 1/8 till 16 dagar.
• F¨or niokroppsproblemet med tiden T = 256 dagar visas i Figur 17a att metoden har r¨att noggrannhet f¨or de inre planetbanorna vid stegl¨angder upp till 4 dagar.
• F¨or niokroppsproblemet med tiden T = 8192 dagar visas i Figur 17b att metoden har r¨att noggrannhet f¨or de yttre planetbanorna vid stegl¨angder upp till 2 dagar.
Tionde ordningens metod
• F¨or tv˚akroppsproblemet med tiden T = 256 dagar visas i Figur 10b att metoden har r¨att noggrannhetsordning f¨or stegl¨angder mellan 2 till 32 dagar.
• F¨or niokroppsproblemet med tiden T = 256 dagar visas i Figur 18a att metoden har r¨att noggrannhet f¨or de inre planetbanorna vid stegl¨angder upp till 8 dagar.
• F¨or niokroppsproblemet med tiden T = 256 dagar visas i Figur 18b att metoden inte gav r¨att noggrannhetsordning. Detta beror antagligen p˚a det korta tidsintervallet d˚a de yttre planeterna inte hunnit f¨orflytta sig lika mycket som de inre.
Som f¨orv¨antat gav de numeriska standardmetoderna fram˚at Euler och bak˚at Euler ett stort fel i l¨osningen. Detta g˚ar att betrakta i Figur 3 d¨ar det tydligt syns att metoderna inte lyckas med en noggrann planetbana f¨or jorden. Dessutom kan det uppskattade felet betraktas i Tabell 1 som visar att redan efter T = 5 ˚ar och med en l˚ag stegl¨angd p˚a h = 1/2 dag s˚a har fram˚at Euler och bak˚at Euler gett ett stort uppskattat fel p˚a st¨orre ¨an 1 AU.
Runge-Kutta 4 som ocks˚a ¨ar en numerisk standardmetod men med h¨ogre noggrannhetsordning var betydligt b¨attre p˚a att l¨osa en planetbana till jorden vilket g˚ar att se i Figur 4. Trots det gav metoden, efter att ha k¨ort i T = 180 000 ˚ar med h = 1 dag ett stort uppskattat fel p˚a 1,569 AU vilket g˚ar att se i Tabell 1.
De sympektiska metoderna symplektisk Euler och St¨ormer-Verlet anv¨andes ocks˚a f¨or att l¨osa problemet med samma tidsintervall och stegl¨angd som Runge-Kutta 4. Felet hos symplektisk Euler g˚ar att se i Tabell 3 d¨ar det uppskattade felet ¨ar 1,997 AU vilket ¨ar ett st¨orre fel ¨an felet fr˚an Runge-Kutta 4. D¨aremot g˚ar det uppskattade felet fr˚an St¨ormer-Verlet att l¨asa av i samma tabell till 0,8299 AU vilket ¨ar ett mindre fel ¨an fr˚an Runge-Kutta 4 trots att metoden har en l¨agre noggrannhetsordning. Skillnaden ¨ar att St¨ormer-Verlet ¨ar en symplektisk metod. I Figur 7 plottas l¨osningen i xy-planet f¨or symplektisk Euler och St¨ormer-Verlet. Planetbanorna ser ut att ha tjocka streck vilket ¨ar l¨osningskurvan som under tidsintervallet skiftar mellan olika radier vid samma punkt i banan. Detta st¨ammer inte ¨overens med jordens riktiga planetbana och f¨or korta tidsintervall s˚a ger de symplektiska metoderna av ordning ett och tv˚a ett st¨orre fel ¨an en standardmetod av ordning fyra som Runge-Kutta 4. D¨aremot ger Runge-Kutta 4 ett globalt fel som v¨axer med tiden och efter T = 180 000 ˚ar har detta fel blivit st¨orre ¨an felet fr˚an St¨ormer-Verlet, en symplektisk metod av ordning tv˚a.
Aven fj¨¨ arde ordningens metod anv¨andes med samma tidsintervall och stegl¨angd som Runge- Kutta 4 f¨or tv˚akroppsproblemet. D˚a dessa metoder har samma noggrannhetsordning var det speciellt intressant att se resultatet. Fr˚an Figur 8a g˚ar det att se jordens planetbana ber¨aknad
av fj¨arde ordningens metod. Denna planetbana har inte lika stort problem med de skiftande radierna, som symplektisk Euler och St¨ormer-Verlet hade, utan planetbanan ser enhetlig ut precis som i l¨osningen av Runge-Kutta 4. Det uppskattade felet fr˚an fj¨arde ordningens metod finns i Tabell 3 och blev 0,0760 AU vilket ¨ar ett betydligt mindre fel ¨an felen fr˚an Runge-Kutta 4, St¨ormer-Verlet och symplektisk Euler. S˚aklart ¨ar h¨og noggrannhet ¨onskv¨art f¨or att f˚a ett litet fel men detta visar ¨aven p˚a att den symplektiska egenskapen ¨ar viktigt f¨or att f˚a ett litet fel f¨or l˚anga tidsintervall vid l¨osning av ett Hamiltonianskt system.
Tionde ordningens metod l¨ostes i endast T = 10 ˚ar och med stegl¨angd h = 1 dag. Vid detta tidsintervall blev det uppskattade felet 1.7155 · 10−18 AU. Felet ¨ar betydligt mindre ¨an fr˚an de f¨oreg˚aende metoderna vilket f¨orklaras av att tidsintervallet var kortare och att metoden ¨ar symplektisk och har h¨og noggrannhet.
F¨or niokroppsproblemet anv¨andes ist¨allet tidsintervallet T = 200 ˚ar f¨or alla de symplektiska metoderna. Felen f¨or alla himlakropparna fr˚an symplektisk Euler, St¨ormer-Verlet och fj¨arde ordningens metod med stegl¨angd h = 1 dag och f¨or tionde ordningens metod med stegl¨angden h = 8 dagar, visas i Tabell 5 och 6. Felen f¨or de olika himlakroppara minskar ju h¨ogre nog- grannhetsordning metoden har. Dock f¨or Merkurius g˚ar det att se att felet ¨okar lite f¨or fj¨arde ordningens metod j¨amf¨ort med St¨ormer-Verlet. Detta var inte ett f¨orv¨antat resultat men troli- gen s˚a har St¨ormer-Verlet gett ett l¨agre fel ¨an f¨orv¨antat f¨or Merkurius vid just tidsintervallet T = 200 ˚ar. Felen f¨or planetbanorna blev mindre f¨or de yttre planeterna j¨amf¨ort med de inre planetbanorna. Till exempel med symplektisk Euler blev felet f¨orSaturnus planetbana 0,0028 AU medan felet f¨or Jordens planetbana blev 0,22 AU. Detta beror p˚a att de yttre planeterna har en mycket l¨angre planetbana d¨ar Jorden har hunnit fler varv ¨an Saturnus under tidsintervallet och d˚a f˚ar ett st¨orre fel.
4.2 Energikonservering
De numeriska standardmetoderna ¨ar inte symplektiska vilka inneb¨ar att energin i l¨osningen inte konserverats. Detta kan ses tydligt f¨or fram˚at Euler och bak˚at Euler i Figur 19 d¨ar energin ¨okar drastiskt f¨or fram˚at Euler medan den minskar f¨or bak˚at Euler. Detta kan ocks˚a ses tydligt i Figur 3, d¨ar jorden inte stannar i sin planetbana.
Runge-Kutta 4 ¨ar en explicit osymplektisk metod och kommer allts˚a inte att konservera energin i l¨osningen. Detta ¨ar sv˚art att se genom att endast studera l¨osningen av planetbanan i Figur 4 d¨ar det ser ut som att jorden lyckas h˚alla sin planetbana. D¨aremot g˚ar det att se i Figur 20a att energin i l¨osningen hela tiden avtar med tiden. Dessutom g˚ar det att studera felet fr˚an Tabell 1 och se att Runge-Kutta 4 ger ett betydligt st¨orre fel ¨an en symplektisk metod med ordning fyra.
Vid inzoomning av Runge-Kutta 4 med stegl¨angden h = 1 dag och i tiden T = 5 ˚ar syns det ocks˚a tydligt att jorden r¨or sig in˚at vilket kan ses i figuren nedan.
(a) Hela planetbanan (b) Inzoomad Runge-Kutta 4 Figur 26: Runge-Kutta 4 med T = 5 ˚ar och h = 1 dag
