1. No category

Simulering av solsystemet
Datorlab med MATLAB
Daniel Vågberg
Institutionen för fysik
Umeå Universitet
17 april 2013
Innehåll
Introduktion
3
Redovisning
3
Simulering av Newtons rörelseekvationer
4
Gravitation
6
Uppgifter
Uppgift 1: Simulering av satellit .
Uppgift 2: Omloppstiden . . . . .
Uppgift 3: Tvåkroppsproblemet .
Uppgift 4: Solsystemet . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
8
. 8
. 9
. 10
. 11
Appendix
Verlet-integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Velocity-Verlet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Exempelkod . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
12
12
12
14
Introduktion
Er uppgift är att simulera ett planetsystem i MATLAB och mäta olika egenskaper på systemet. Instruktionerna börjar med en genomgång av de numeriska
metoder som ska användas och fortsätter sedan med hur man tillämpar dem på
den typ av gravitationsproblem vi är intresserade av. Efter det kommer själva
uppgifterna.
Syftet med uppgifterna är dels att ni ska bli bättre på att använda MATLAB
till fysikaliska simuleringar och dels att ni ska få undersöka ett fysikaliskt system
och se hur olika kvantiteter (tex energi och rörelsemängd) bevaras över tid. Till
vissa uppgifter kommer ni behöva data för planeterna tex jordens och solens
massa, dessa kan tas ifrån tex Physics Handbook, ange alltid vilka värden ni
använder och varifrån de kommer.
Redovisning
Ni jobbar två och två och lämna in en skriftlig redovisning tillsammans. Rapporten behöver inte följa en strikt rapportmall, men det måste tydligt gå att
följa och förstå vad ni gjort.
• Rapporten ska innehålla svar på samtliga frågor. Svaren ska vara tillräck-
ligt utförliga så man kan förstå dem utan att ha tillgång till själva frågan.
• Figurer som efterfrågas i uppgifterna ska nnas med i rapporten. Ni får
naturligtvis även lägga till egna gurer där ni tycker det passar.
• Den MATLAB-kod ni använder för att lösa problemen ska nnas med i
rapporten. Längre kod kan med fördel läggas i ett appendix i slutet av
rapporten.
• Glöm inte att kommentera MATLAB-koden.
• Ange vilka data ni använder för planeterna och varifrån de kommer.
• Rapporten skickas in som pdf-l
• Glöm inte att skriva namn och email på rapporten!
3
Simulering av Newtons rörelseekvationer
Vi börjar med att studera rörelse i en dimension, när vi löst det endimensionella
fallet är det enkelt att generalisera lösningen till två eller tre dimensioner. Antag
att vi har ett föremål med position x, hastighet v och acceleration a. Vi vill
beräkna hur föremålets position förändras över tiden, vi vet föremålets position
och hastighet vid tiden t = 0.
Accelerationen hos ett föremål beskrivs av Newtons andra lag
F = ma → a =
F
,
m
(1)
där F är kraften som påverkar föremålet, m är föremålets massa och a är
föremålets acceleration. Accelerationen denierats som andraderivatan av positionen
d2 x
(2)
a= 2.
dt
Om vi känner till systemets position och hastighet vid tiden t = 0 och accelerationen (eller kraften) är en känd funktion kan vi lösa ovanstående dierentialekvation och se hur positionen x varierar med tiden. För vissa problem går
ekvationen att lösa för hand, men för mer avancerade problem eller problem
med många föremål som rör sig samtidigt är ofta enda lösningen att använda
datorhjälpmedel.
Ekvation 2 är en andra ordningens dierentialekvation. Den går att skriva
om som två kopplade första ordningens dierentialekvationer genom att införa
hastigheten v som en hjälpvariabel. Problemet som ska lösas blir då



dx(t)
dt
dv(t)
dt
= v(t)
= a(t)
(3)
där initialvärden x(0) = x0 och v(0) = v0 är kända, och a(t) är en känd funktion.
Både ekvation 2 och 3 beskriver samma fysik. Numeriskt är det ofta lättare att
hantera första ordningens dierentialekvationer så vi kommer använda ekvation
3 som grund för våra simuleringar. Notera att den numeriska metoden vi ska
använda förutsätter att kraften F är konservativ.
För att lösa problemet numeriskt måste vi diskretisera problemet. Vi börjar
med att dela in tiden i diskreta punkter t0 , t1 , t2 ..., tn , tn+1 , ... mellan varje punkt
är det en konstant tidsskillnad ∆t så att t0 = 0, t1 = ∆t, ..., tn = n∆t. För att
förenkla notationen kommer vi använda beteckningen xn för x(tn ) och v n för
v(tn ) osv. Det den numeriska metoden kommer göra är att hitta en approximativ
lösning till ekvation 3. Grundidén är att om vi vet var vi är vid tiden tn kan
vi approximera var vi kommer vara vid tiden tn+1 och på så sätt stega oss
framåt i tiden. Detta görs genom att approximera derivator, vi kan till exempel
approximera förändringen i x vi tiden tn med
x(tn + ∆t) − x(tn − ∆t)
xn+1 − xn−1
dx
≈
=
dt
2∆t
2∆t
(4)
Steglängden ∆t avgör hur bra approximationen blir och vi ser att om ∆t → 0
blir derivatan exakt. Men om vi tar väldigt små tidssteg kommer simuleringen
att ta lång tid eftersom vi måste ta er steg för att simulera samma mängd tid.
4
Konsten är att välja ett ∆t som är litet nog för att ge en bra approximation
men samtidigt inte för litet så att simuleringen tar orimligt lång tid att köra.
Integrationsmetoden vi ska använda heter Velocity-Verlet, och bygger på två
rekursiva uppdaterings formler som används för att stega framåt i tiden. De två
formlerna en för positionen och en för hastigheten återges här nedan

 xn+1
 v n+1
= xn + v n ∆t + 21 an ∆t2
= v n + 21 an + an+1 ∆t
(5)
En härledning av Velocity-Verlet metoden nns i appendix efter uppgifterna,
läs gärna igenom den, det är bra om man vet vilka antagnaden som ligger
bakom en metod man använder så man inte råkar ut för överraskningar. Det är
trivialt att utöka systemet till er dimensioner, rörelsen i de olika riktningarna
är bara kopplade genom kraften som är positions beroende i övrigt är rörelsen
i de olika riktningarna oberoende av varandra. För två dimensioner får vi fyra
uppdaterings ekvationer x, y , vx och vy



xn+1




 y n+1


vxn+1




 v n+1
y
= xn + vxn ∆t + 21 anx ∆t2
= y n + vyn ∆t + 21 any ∆t2
∆t
= vxn + 12 anx + an+1
x
= vyn + 12 any + an+1
∆t
y
(6)
Mer allmänt kan vi skriva om uppdaterings relationerna på vektorform,

 rn+1
 vn+1
= rn + vn ∆t + 21 a(rn )∆t2
= vn +
1
2
a(rn ) + a(rn+1 ) ∆t
(7)
där r är positionsvektorn, v är hastighets vektorn och a(r) är en känd funktion
som beräknar accelerationsvektorn utifrån en given positionsvektor. Ekvationerna ovan gäller för rörelsen av en enskild partikel/föremål, men vi kan utöka
dem till att gälla era partiklar. Antag att vi har två partiklar med position r1
och r2 och som växelverkar genom accelerationen dvs a(r1 , r2 ) då får vi följande



rn+1

1



 rn+1
2

 v1n+1




 vn+1
2
=
rn1 + v1n ∆t + 21 a1 (rn1 , rn2 )∆t2
=
rn2 + v2n ∆t + 21 a2 (rn1 , rn2 )∆t2
=
v1n +
1
2
=
v2n +
1
2
, rn+1
) ∆t
a1 (rn1 , rn2 ) + a1 (rn+1
1
2
a2 (rn1 , rn2 ) + a2 (rn+1
, rn+1
) ∆t
1
2
(8)
Har man er partiklar är det bara att fortsätta lägga till er ekvationer. Det
ser mycket ut när man skriver ut ekvationerna så här men när man programmerar behöver man bara skriva in ekvationerna en gång, har man era partiklar
löser man det med en loop eller ännu bättre ordnar datat i vektorer så man
kan använda MATLABs inbyggda vektoroperationer istället. Dvs när ni skriver programmet är det bara de fyra uppdateringsformlerna i ekvation 6 som ni
behöver använda, har ni era partiklar så löser ni det genom att vektorisera
ekvationerna.
5
Gravitation
Hittills har vi talat ganska allmänt om hur man simulerar Newtons-rörelseekvationer,
men nu ska vi bli lite mer specika och fokusera på det problem vi ska jobba
med dvs gravitationkrafter och planetbanor. Det vi behöver för att använda
Velocity-Verlet-metoden är ett uttryck för accelerationen och hur den varierar
med positionen. Eftersom a = F/m börjar vi med att specicera vilka krafter
vi har. Antag att vi har två planeter i och j med massa mi och mj . Gravitationskraften Fij mellan planeterna ges av
Fij = G
mi mj
ˆrij ,
2
rij
(9)
där G är gravitationskonstanten, rij är avståndet mellan objekten och ˆrij är en
enhetsvektor som pekar från planet i mot planet j . Eftersom kraft och motkraft
måste balansera så vet vi att Fij = −Fji . Vi ska lösa problemet i två dimensioner
så rij ges av
q
rij = |rj − ri | = (xj − xi )2 + (yj − yi )2
(10)
För att kunna använda lösningsmetoden i ekvation 6 måste vi komposantuppdela kraften. I gur 1 ser vi kraftvektorn ut ritad för planet i. De två
kraftkomposanterna på planet i ges av

 Fx
= −Fij cos θ = −Fij
xi −xj
rij
= −Gmi mj
xi −xj
3
rij
 F
y
= −Fij sin θ = −Fij
yi −yj
rij
= −Gmi mj
yi −yj
3
rij
(11)
från dessa kan vi sedan beräkna accelerations komposanterna

 ax
=
− miji cos θ = − miji
 a
y
=
− miji sin θ = − miji
F
F
xi −xj
rij
= −Gmj
xi −xj
3
rij
F
F
yi −yj
rij
= −Gmj
yi −yj
3
rij
(12)
Vi ser att accelerationen av planet i är oberoende av planetens egen massa mi .
Vi noterar också att vi aldrig behöver räkna ut vinkeln θ. Nu har vi allt vi
y
mj
rij = |rj − ri |
rj
Fx
Fij
θ
Fy
mi
ri
x
Figur 1: Figuren visar de två planeterna i och j , tillsammans med dess positionsvektorer ri och rj , till planet i har även kraftvektorn ritats ut.
6
behöver för att simulera rörelsen av två himlakroppar som påverkar varandra
via gravitationen. Om vi har mer än två planeter så kommer alla planeter att
påverka varandra. För att få den totala kraften som påverkar en planet måste
vi summera kraften från alla andra planeter. Antag att vi har N planeter den
total kraften på planet i ges då av
Fi =
N
X
(13)
Fij
j
Nu när vi har allt vi behöver för att simulera rörelsen för ett godtyckligt antal
himlakroppar ska vi bara snabbt gå igenom olika kvantiteter som vi kan vilja
mäta i vårt system. Lättaste sättet testa att en simulering fungerar som den ska
är att kontrollera att energin bevaras i systemet. Vi har två typer av energi i
systemet, kinetisk energi som ges av
Ek =
N
X
mi v 2
(14)
i
2
i
och potentiellenergi som ges av
Ep = −G
N X
N
X
mi mj
i
rij
j>i
(15)
Den totala energin i systemet E ges av
(16)
Om allt fungera som det ska och systemet inte påverkas av några yttre krafter ska
E vara konstant genom hela simuleringen. En annan kvantitet som också bevaras
förutsatt att inga yttre krafter verkar på systemet är den totala rörelsemängden
p som ges av
E = Ek + Ep
p=
N
X
(17)
mi vi
i
En tredje bevarad kvantitet är rörelsemängdsmomentet
L=
N
X
ri × mi vi
(18)
i
som bevaras om systemet inte påverkas av några vridmoment skapade av yttre
krafter. L är en vektor som är vinkelrätt mot både r och v. För ett tvådimensionellt system, dvs ett systemet som är begränsat till att bara röra sig i
xy -planet, kommer L endast bestå av en z -komponent. Att kontrollera att de
tre kvantiteterna E , p och L bevaras är ett mycket bra test som hjälper en att
hitta eventuella fel i simuleringen. Att de bevaras är dock ingen garanti på att
allt är rätt, men om någon av dem inte bevaras så är det denitivt ett tecken på att något är fel. En annan intressant kvantitet man kan studera är hur
masscentrum för hela systemet rör sig. Vi kan beräkna masscentrums position
genom
1
rCM = PN
i
N
X
mi
i
mi ri
(19)
Om inga yttre krafter verkar på systemet och den totala rörelsemängden p är
noll kommer masscentrum att stå still.
7
Uppgifter
Målet med uppgifterna är att vi ska skriva en simulering där vi kan simulera
rörelsen för solen och de inre planeterna i solsystemet. Vi ska bygga upp simuleringen stegvis och börjar först med några enklare fall som vi sedan kan bygga
vidare på för att få den fulla lösningen.
Uppgift 1: Simulering av satellit
Vi ska börja med att simulera en satellit i omloppsbana kring en planet. Vi antar
att planeten väger mycket mer än satelliten och att planetens rörelse därför inte
kommer påverkas nämnvärt av satelliten. I vår simulering kommer vi därför anta
att planeten står stilla och att det bara är satelliten som rör sig. Detta är förstås
en approximation. I nästa uppgift ska vi även simulera planetens rörelse och se
hur den påverkar resultatet.
y
Fx
m
F
Fy
r
θ
x
M
Figur 2: En liten satellit i omloppsbana kring en planet. Vi ska simulera systemet
under antagandet att m M .
• Skriv en funktion orbit_1body som simulerar banan för en satellit med
massa m i omloppsbana kring en planet med massa M . Vi antar att m M så att planetens rörelse kan försummas. Planeten står still, och är
placerad i origo. Rörelseekvationerna ska integreras med hjälp av VelocityVerlet. Funktionen ska använda följande funktionshuvud:
function [x,y,vx,vy,t]=orbit_1body(G,M,x0,y0,vx0,vy0,dt,tmax)
där G är gravitations konstanten, M är massan hos planeten, x0, y0, vx0 och
vy0 beskriver satellitens position och hastighet vid tiden t = 0. Variabeln
dt är längden på tidsteget och tmax är den totala tiden att simulera dvs vi
ska simulera systemet från t = 0 till t =tmax. Notera att satellitens bana
är helt oberoende av satellitens egen massa m. Funktionen returnerar fem
vektorer x och y som innehåller satellitens koordinater en datapunkt för
8
varje tidsteg i simuleringen, vx och vy som innehåller satellitens hastighet
för varje tidsteg, samt t som innehåller tiden för varje datapunkt. På sista
sidan i instruktionerna nns ett kodexempel som går att använda som
utgångspunkt om ni känner er osäkra på hur ni ska börja.
• Testkör funktionen med följande parametrar och initialvärden: G=1, M=10,
m=0.01, x0=10 y0=0 vx0=0, och vy0=0.75, välj tmax så att satelliten
hinner ca 5-10 varv runt planeten under simuleringen.
Rita upp satellitens bana, och markera planetens position i guren.
Simulera med olika värden på dt och se hur noggrannheten i simuleringen förändras. Vilket värde på dt verkar lämpligt att använda?
Kontrollera att energin bevaras i simuleringen, plotta Ek , Ep och
Ep + Ek och se hur de förändras med tiden.
Kontrollera att rörelsemängdsmomentet är bevarat för systemet.
Kontrollera att rörelsemängden bevaras i systemet. Får vi det förväntade resultatet? Om inte förklara vad som händer i systemet.
Uppgift 2: Omloppstiden
Vi ska nu beräkna omloppstiden för satelliten i föregående uppgift utifrån våra
simulerings resultat. Vi kommer behöva räkna ut er omloppstider i de följande
uppgifterna vi ska därför automatisera den processen genom att skriva en funktion som gör jobbet åt oss.
• Skriv en funktion som beräknar omloppstiden givet koordinatvektorerna
x, y och tidsvektorn t som genererats av funktionen orbit_1body. Detta
går att göra på era olika sätt. Beskriv vilken algoritm ni använder i
rapporten. (Tips börja med att plotta x och/eller y mot tiden för att
avgöra vilken egenskap hos kurvorna som kan användas för att beräkna
omloppstiden.)
• Testkör funktionen, använd data med samma initialvillkor som i föregående
uppgift. Vilken omloppstid har satelliten?
• Hur påverkas omloppstiden om satellitens initialhastighet ökar? Vid vilken
hastighet slutar satelliten att gå i omloppsbana runt planeten? Vad är
den totala energin i systemet när det händer? Simulera med samma initialvillkor som tidigare men öka stegvis värdet på vy0. Kontrollera hur
omloppstiden och den totala energin Ek + Ep varierar.
• Vi ska nu testa vår kod på ett verkligt exempel: Rymdstationen ISS (In-
ternational Space Station) ligger i en nästan cirkulär omloppsbana kring
jorden. Omloppsbanan är en så kallad LEO (Low Earth Orbit) vilket innebär att stationens omloppsbanan ligger strax utanför atmosfären på ca
400km höjd över jordytan. Stationen väger ca 450ton och har en medelhastighet av 7700m/s. Simulera stationens rörelse, anpassa dt så att ni
får en stabil lösning. Vilken omloppstid har stationen?
9
Uppgift 3: Tvåkroppsproblemet
Nu ska vi utöka vår simulering så att vi även simulerar planetens rörelse. Eftersom vi simulerar rörelsen hos båda kropparna behöver vi inte göra några antaganden om hurvida satelliten väger mer eller mindre än planeten. Vårt system
kommer nu se ut som i gur 1.
• Skriv en funktion orbit_2body som simulerar rörelsen hos två himlakrop-
par med hjälp av Velocity-Verlet-metoden, utgå från er tidigare simulering.
Funktionen ska använda följande funktionshuvud:
function [x,y,vx,vy,t]=orbit_2body(G,m,x0,y0,vx0,vy0,dt,tmax)
Skillnaden mot tidigare är att vi nu har två kroppar som rörsig, vilket betyder att vi behöver dubbelt så många initialvillkor. Variablerna m, x0, y0,
vx0 och vy0 kommer därför vara vektorer med två element. Till exempel
den initiala x-positionen för den första kroppen anges i det första elementet x0(1) och motsvarande data för den andra kroppen anges i x0(2)
osv. Vektorerna x, y, vx och vy som funktionen returnerar innehåller data
för rörelsen hos båda kropparna och har dimensionen 2×steps där steps
är antalet tidssteg i simuleringen.
• Testkör funktionen med samma initialvärden som i uppgift 1. För planeten
innebär det m(1)=10, x0(1)=0, y0(1)=0, vx0(1)=0 och vy0(1)=0. Satel-
litens initialvärden är samma som i uppgift 1 och placeras på position 2 i
vektorerna tex m(2)=0.01.
Rita upp satellitens och planetens bana i samma gur och jämför
med resultatet från uppgift 1.
Undersök planetens rörelse, har den en sluten bana? om inte vad
beror det på? Testa att öka massan på satelliten till m(2)=1 för att
få en tydligare eekt. Vad händer med masscentrum?
Vad måste vara uppfyllt för att masscentrum ska stå still? Räkna ut
vilken initialhastighet planeten behöver för att masscentrum ska stå
still. Simulera med de nya initialvillkoren och veriera att det fungerar. Hur förändras planetens bana? Gör en gur som visar planetbanan
före och efter förändringen.
Kontrollera att energin är bevarad.
Kontrolera att rörelsemängdsmomentet är bevarat.
Kontrolera att rörelsemängden är bevarad. Rita en gur med tre
kurvor, rörelsemängden för satelliten, rörelsemängden för planeten
och deras totala rörelsemängd.
10
Uppgift 4: Solsystemet
Nu är vi redo att simulera solsystemet. Vi ska utöka simuleringen så att den
kan hantera rörelsen hos N planeter som alla påverkar varandra genom gravitationen. och sedan testa simuleringen genom att simulera de inre delarna av
solsystemet.
• Skriv en funktion force som beräknar kraften för alla N planeterna. Funk-
tionen ska använda följande funktionshuvud
function [f]=force(G,m,x,y)
där f, m, x och y är vektorer med längd N . Där m innehåller planeternas
massor, x och y är planeternas positioner, G är gravitationskonstanten och
f är den totala kraften som verkar på varje planet som räknas ut genom
att summera kraften från alla andra planeter enligt ekvation 13.
• Skriv en funktion orbit_Nbody som simulerar rörelsen hos N st him-
lakroppar, med hjälp av Velocity-Verlet-metoden, utgå från er tidigare
simulering och använd funktionen force för att beräkna kraften mellan
planeterna. Funktionen ska använda följande funktionshuvud:
function [x,y,vx,vy,t]=orbit_Nbody(G,m,x0,y0,vx0,vy0,dt,tmax)
In parametrarna är samma som tidigare med skillnaden att massan och
initialvärdena nu är vektorer med längd N eftersom vi nu har N kroppar
som behöver initieras. Av samma anledning har nu returvärdena dimension N ×steps för att kunna hålla banorna för de N himlakropparna vi
simulerar.
• Testkör simuleringen med två kroppar och samma initialvillkor som i
uppgift 3. Kontrollera att det blir samma resultat som i uppgift 3.
• Simulera solsystemet! Simulera solen och de inre planeterna Merkurius,
Venus, Jorden och Mars. Dvs vi har N = 5, använd Physics Handbook
eller annan källa för att hitta lämpliga initialvärden för planeterna. För
att förenkla valet av initialvärden kan vi anta att alla planeterna ligger
på en rad vid tiden t = 0. Solens hastighet bör väljas så att den totala
rörelsemängden i systemet blir noll. För den här uppgiften räcker det om
ni använder planeternas medelhastighet och medlavstånd från solen som
initialvärden, men då kommer banorna att bli cirkulära istället för elliptiska. (För att få korrekta elliptiska banor behöver vi veta avståndet till
solen och planetens hastigheten i en specik punkt på banan i stället för
medelvärdet beräknat över hela banan)
Välj tidsteg dt så att alla planeterna får en stabil bana, hur kort
tidsteget måste vara bestäms av den snabbaste rörelsen i systemet, i
det här fallet Merkurius.
Gör en gur som visar planeternas banor.
Kontrollera att energi, rörelsemängd och rörelesmängdsmoment bevaras. Gör tre gurer en för vardera kvantitet, gurerna ska innehålla
en kurva för varje himlakropp och en kurva med det totala värdet.
Beräkna omloppstiderna och jämför med tabellvärden.
11
Appendix
Verlet-integration
Vi ska nu härleda en metod för att integrera accelerationen och beräkna banan
som ett föremål kommer att följa, vi antar att vi vet föremålets position och
hastighet vid t = 0. Vi utgår ifrån Newtons andra lag
F (x(t))
d2 x
= a(x(t)) =
dt2
m
(20)
där F (x(t)) är en konservativ kraft som endast beror på positionen x. Vi diskretiserar tiden och approximerar andraderivatan enligt följande
a(x) =
d2 x
≈
dt2
xn+1 −xn
∆t
−x
∆t
n
−xn−1
∆t
=
xn+1 − 2xn + xn−1
= an
∆t
(21)
Om vi löser ut xn+1 får vi
xn+1 = 2xn − xn−1 + an ∆t2
(22)
vilket är en fullt fungerande integrations metod som kan användas för att räkna
ut framtida positioner xn+1 givet de två föregående positionerna xn och xn−1
och accelerationen. Metoden har utvecklats era gånger av olika forskare genom
historien men kallas oftast Verlets metod efter den senaste upptäckaren som
gjorde metoden känd men man kan även se andra namn som tex Störmers metod.
Verlets metod använder bara positionen och accelerationen för att räkna ut nästa
position, hastigheten används inte. Detta har både för och nackdelar beroende
på vilken typ av problem vi vill lösa. Om vi inte behöver känna till hastigheten
är metoden mycket eektiv eftersom det går att räkna ut positionen direkt utan
att gå omvägen via hastigheten. Men om vi behöver hastigheten blir det lite
omständigt. Hastigheten går att räkna ut från skillnaden mellan två positioner:
vn =
xn+1 − xn−1
2∆t
(23)
men vi ser att för att räkna ut den nuvarande hastigheten v n måste vi redan
känna till den nya positionen xn+1 , vilket inte alltid är praktiskt. En annan sak
att tänka på är att oftast känner man bara till en position vid t = 0 så för
att komma igång måste xn−1 först beräknas med hjälp av någon annan metod.
För att metoden ska fungera måste ∆t vara konstant under hela simuleringen.
Ändrar man tidsteget under simuleringen kommer partikelns rörelse att bli felaktig om man inte samtidigt skalar om de andra termerna för att kompensera
för förändringen.
Velocity-Verlet
Det går att skriva om Verlets metod så att man får ut hastigheten direkt, metoden kallas då Velocity-Verlet. Det är en populär metod som ofta används för att
simulera rörelsen hos tex molekyler. Vi utgår ifrån den vanliga Verlet metoden
och börjar med med att lösa ut xn−1 från ekvation 23, vi får då
xn−1 = xn+1 − 2v n ∆t
12
(24)
Vi sätter in uttrycket för xn−1 i ekvation 22 vilket ger
1
xn+1 = xn + v n ∆t + an ∆t2
2
(25)
som är en ny uppdaterings metod för positionen som även beror på hastigheten.
Men för att få en fungerande metod behöver vi även veta hur vi ska uppdatera
hastigheten. Utifrån ekvation 23 kan vi sätta upp ett uttryck för v n+1 :
v n+1 =
xn+2 − xn
2∆t
(26)
Vi kan ta fram uttryck för xn+2 och xn från ekvation 25 och sätta in dem i
ekvation 26 vilket ger
v
n+1
xn+1 + v n+1 ∆t + 21 an+1 ∆t2 − xn+1 − v n ∆t − 21 an ∆t2
=
2∆t
(27)
Vilket efter förenkling blir
v n+1 = v n +
1 n
a + an+1 ∆t
2
(28)
Vi har nu två uppdateringsekvationer:

 xn+1
 v n+1
= xn + v n ∆t + 21 an ∆t2
= v n + 21 an + an+1 ∆t
(29)
som tillsammans utgör metoden vi kallar Velocity-Verlet. Metoden är själv startande, dvs vi behöver inte använda en separat metod för att beräkna xn−1
eftersom det värdet aldrig används till skillnad från i den vanliga Verlet metoden. Även Velocity-Verlet kräver ett konstant ∆t för att fungera korrekt. Man
bör också vara medveten om att metoden bygger på antagandet att kraften
är konservativ, dvs att kraften går att skriva som gradienten av en potential,
F = −∇V (r), där r är positions vektorn, och V är en funktion som beskriver
den potentiella energin i systemet. De krafter (gravitation) vi ska simulera är
konservativa så det antagandet är ingen begränsning för oss i det här fallet.
Men problem kan uppstå om man har icke konservativa krafter som till exempel
friktion eller luftmotstånd.
Slutligen ska vi bara nämna att det nns en annan populär metod som kallas
Leapfrog som också går att härleda genom att göra en omskrivning av Verletmetoden. Velocity-Verlet och Leapfrog är väldigt lika och har i princip samma
egenskaper.
13
Exempelkod
Här nedan återges ett exempel på strukturen i ett program som använder stegningsmetoder liknande Velocity-Verlet. Koden är tänkt att simulera en planet
med rörelse i två dimensioner x och y . Syftet med koden är att ge lite tips om
hur ni kan börja med uppgiften, men för att inte göra uppgiften för lätt är era
rader kapade och avslutas istället med ...
%starting point for a simple simulation program
%preallocate memory (increases performance)
steps=... %select number of time steps
x=zeros(1,steps);
y=zeros(1,steps);
vx=zeros(1,steps);
vy=zeros(1,steps);
%set initial conditions
x(1)=...
y(1)=...
vx(1)=...
vy(1)=...
%define functions for calculating acceleration based on position
ax=@(...)...
ay=@(...)...
%simulate orbit
for i=1:steps
x(i+1)=...
y(i+1)=...
vx(i+1)=...
vy(i+1)=...
end
%
%update position and velocity
%using Velocity-Verlet
%
%plot and analyse results
14