Simulering av solsystemet

Datorlab med MATLAB Daniel Vågberg Institutionen för fysik Umeå Universitet 17 April 2013 Editerat: Pontus Svenmarker 24 mars 2014

1

Innehåll

1 Newtons rörelsekvationer med MATLAB 1.1 Gravitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.2 Newtons rörelseekvationer med Velocity-Verlet . . . . . . . . . .

2 Uppgifter 2.1 Kraft och acceleration mellan kroppar . . . . . . . . . . . . . . .

2.2 En satellit i omloppsbana . . . . . . . . . . . . . . . . . . . . . .

2.3 Jorden och månen . . . . . . . . . . . . . . . . . . . . . . . . . .

7 7 7 9 2.4 Solsystemet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3 3 4 A Verlet-integration B Exempelkod 12 14 2

1 Newtons rörelsekvationer med MATLAB

Fysikaliska beräkningar med datorhjälp har fått en allt mer central roll i civilin genjörsutbildningen i och med att datorkraften ökat samtidigt som användarvän ligheten förbättrats. Idag har vanliga persondatorer tillräckligt med datorkraft för att genomföra relevanta fysikaliska simuleringar. Traditionellt används da torns CPU (Central Processing Unit) för att göra beräkningar, men i och med att utvecklingen av grakkort idag går snabbare än utvecklingen av CPU har det även blivit intressant att använda datorns grakkort för beräkningar. Användar vänligheten utvecklas konstant vilket minskar insteget till att göra beräkningar med hjälp av datorer.

Er uppgift är att simulera vårt solsystem i MATLAB. Valet av MATLAB är medvetet då det erbjuder exibilitet, en hög nivå (behöver inte skriva så många rader kod) och bra hjälp och support. Syftet med datorlaborationen är främst att ni ska träna er färdighet i att använda MATLAB till fysikaliska simuleringar, men även att öka er förståelse för gravitation.

Denna labinstruktion börjar med en genomgång av gravitation och fortsätter med de numeriska metoder som ska användas. Sist kommer fyra uppgifter, en uppgift ska lösas för hand och tre stycken med hjälp av datorn.

1.1 Gravitation

Gravitation är det fenomen som får kroppar att attrahera varandra. I vår vardag möter vi gravition främst som det som ger saker tyngd. Vi ska nu lyfta blicken mot vårt solsystem och förstå hur vi kan förklara planeternas rörelse runt solen med den teori som Isac Newton har presenterat om gravitation.

Tänk dig två kroppar (solen och jorden) med massan attraherande gravitationskraften F ij m i respektive som skapas mellan kropparna ges av m j som attraheras av varandra, schematiskt illustrerade i Fig. 1. Storleken på den F ij = G m i m r 2 ij j ˆ ij , (1) där G är gravitationskonstanten, r ij enhetsvektor som pekar från kropp i måste vara balanserade vet vi att F ij = − F ji .

Om vi har er än två kroppar i vårt system, kommer alla kroppar att påverka varandra. För att få den totala kraften som påverkar en kropp måste vi summera kraften från alla andra kroppar. Om vi antar att vi har N kroppar, då blir den total kraften på kropp i är avståndet mellan objekten och mot kropp ˆ ij är en j . Eftersom kraft och motkraft N F i = X F ij .

(2) j För ett system bestående av tre kroppar är alla krafter utritade i Fig. 2. Totalt sex stycken krafter verkar i systemet, men enbart två per kropp. För kropp 1 kan vi alltså skriva N F 1 = X F ij = F 12 + F 13 .

(3) j Nu har vi allt vi behöver för att förstå rörelsen för ett godtyckligt antal himlakroppar. Det nns ett antal storhter, utöver krafter, som vi vill kunna 3

y r j m j r ij = | r j − r i | θ F ij F y F x m i r i x Figur 1: Två kroppar med massan tionsvektorer r i och m i respektive m j tillsammans med deras posi r j . Gravitationskraften mellan kropparna F ij är utritad i guren tillsammans med kraftkomposanterna i riktningarna och ˆ .

mäta i vårt system. Vi har två typer av energi i systemet; kinetisk energi som ges av N E k = X i m i 2 v i 2 (4) och potentiell energi som ges av E p = N N − G X X m i m j .

r ij i j>i (5) Den totala energin i systemet E ges av E = E k + E p .

(6) Om systemet inte påverkas av några yttre krafter ska den totala energin vara konstant. En annan storhet som också bevaras, förutsatt att inga yttre krafter verkar på systemet, är den totala rörelsemängden p = N X m i v i .

i (7) Ytterligare en intressant storhet man kan studera är hur masscentrum för hela systemet rör sig. Vi kan beräkna masscentrums position genom r CM 1 = P N i m i N X m i r i .

i (8) Om inga yttre krafter verkar på systemet och den totala rörelsemängden noll kommer masscentrum att stå still.

p är

1.2 Newtons rörelseekvationer med Velocity-Verlet

Hittills har vi talat allmänt om gravitation, men nu ska vi bli lite mer spe cika och fokusera på det problem vi ska jobba med gravitationkrafter och 4

m 2 F 21 m 3 F 31 F 12 m 1 Figur 2: Tre kroppar med olika massor där gravitationskraften mellan alla kroppar är utritade.

planetbanor. Vi börjar med att studera rörelse i en dimension och när vi löst det endimensionella fallet, kan vi enkelt 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, givet att 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 , (9) där F enligt är kraften som påverkar föremålet, m är föremålets massa och a är före målets acceleration. Accelerationen denieras som andraderivatan av positionen d 2 x a = dt 2 .

(10) Om vi känner till kroppens position och hastighet vid tiden ofta enda lösningen att använda datorhjälpmedel.

t = 0 samt om dess acceleration (eller den totala kraften som verkar på den) ä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 Ekvation (10) ä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 hjälpvariabel. Problemet som ska lösas blir då ( dx ( t ) dt dv ( t ) dt = v ( t ) = a ( t ) (11) där initialvärden x (0) = x 0 och v (0) = v 0 är kända, och metoden vi ska använda förutsätter att kraften F a ( t ) är en känd funk tion. Både ekvation (10) och (11) beskriver samma fysik. Numeriskt är det ofta lättare att hantera första ordningens dierentialekvationer så vi kommer använ da ekvation (11) som grund för våra simuleringar. Notera att den numeriska är konservativ. En konser vativ kraft är en sådan kraft som utför ett arbete som enbart är beroende av 5

dess start och slutposition. Gravitationskraft och coulumbkraft är exempel på konservativa krafter, medan friktionskraft är ett exempel på en dissipativ kraft.

För att lösa problemet numeriskt måste vi diskretisera det. Vi börjar med att dela in tiden i diskreta punkter t 0 , t 1 , t 2 ..., t n , t n +1 , ...

. Mellan varje punkt råder konstant tidsskillnad ∆ t så att t 0 = 0 , t 1 = ∆ notationen kommer vi använda beteckningen t , ..., x n för (11). Grundidén är att om vi vet var vi är vid tiden t n = x ( t n ) t n n ∆ t . För att förenkla och v n för v ( t n ) osv.

Den numeriska metoden kommer att hitta en approximativ lösning till ekvation kan vi approximera var vi kommer vara vid tiden vi genom att approximera tidsderivatan, alltså en förändring av positionen med avseende på tid i en viss punkt. Vi kan tillexempel approximera förändringen i x vid tiden t n med t n +1 och på så sätt stega oss framåt i tiden. Detta gör dx ≈ dt x ( t n + ∆ t ) − x ( t n − ∆ t ) 2∆ t = x n +1 − 2∆ t x n − 1 (12) Steglängden ∆ t avgör hur bra approximationen blir och vi ser att om 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 en given total tid.

Konsten är att välja ett ∆ t formlerna, en för positionen och en för hastigheten, är ∆ t → men samtidigt inte för litet så att simuleringen tar orimligt lång tid att köra.

0 som är litet nog för att ge en bra approximation, Integrationsmetoden vi ska använda heter Velocity-Verlet, och bygger på två rekursiva uppdateringsformler som används för att stega framåt i tiden. De två x v n n +1 +1 = = x n v n + + v 1 2 n ∆ a n t + 1 2 a n ∆ + a n +1 t 2 ∆ t (13) En härledning av Velocity-Verlet-metoden nns i appendix 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 syste met till er dimensioner; rörelsen i de olika riktningarna är bara kopplade genom kraften som är positionsberoende och i övrigt är rörelsen i de olika riktningarna oberoende av varandra. För två dimensioner får vi fyra uppdateringsekvationer x , y , v x och v y , enligt        x n +1 y n +1 v n x +1 v n y +1 = = = = x n y n v n x v n y + v n x ∆ t + + v n y ∆ t + + + 1 2 1 2 a n x a n y 1 2 1 2 a n x ∆ t 2 a n y ∆ t 2 + a n x +1 ∆ t + a n y +1 ∆ t .

Mer allmänt kan vi skriva om uppdateringsrelationerna på vektorform som (14) r v n n +1 +1 = = r n v n + + v 1 2 n ∆ t + 1 2 a a ( r n ) + ( r a n ( r )∆ n t +1 2 ) ∆ t , (15) 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. Ekvatio nerna ovan gäller för rörelsen av en enskild kropp, men vi kan utöka dem till att gälla era kroppar. Antag att vi har två kroppar med positionerna r 1 och 6

r 2 som växelverkar genom accelerationen, dvs tionssystem: a ( r 1 , r 2 ) . Då får vi följande ekva        r n 1 +1 r n 2 +1 v n 1 +1 v n 2 +1 = = = = r n 1 r n 2 v n 1 v n 2 + + + v 1 2 1 2 n 1 ∆ + v n 2 ∆ t + a 1 ( r n 1 , r n 2 ) + a 1 ( r n +1 a t 2 ( + r n 1 1 2 1 2 , a a 2 ( r n 1 , r n 2 )∆ t 2 r 1 n 2 ( r n 1 ) + , r a n 2 2 ( )∆ r 1 n 1 t 2 +1 , r n 2 +1 , r n 2 +1 ) ) ∆ t ∆ t .

(16) Har man er kroppar ä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 kroppar i sitt system löser man det med en loop, eller ännu hellre, ordnar datat i vektorer så man kan använda MATLABs inbyggda vektoroperationer. När ni skriver programmet är det med andra ord bara de fyra uppdateringsformlerna i ekvation (14) som ni behöver använda. Har ni era kroppar i ert system så löser ni det genom att vektorisera ekvationerna.

2 Uppgifter

Målet med uppgifterna är att simulera rörelsen för solen och de inre planeterna i vårt solsystem. Vi kommer att bygga upp simuleringen i steg och starta med ett enkelt fall, för att sedan bygga vidare på det till den fullständiga lösningen.

2.1 Kraft och acceleration mellan kroppar

Som bekant säger Newtons teori att gravitationskraften mellan två kroppar beskrivs av ekvation (1). Ekvationen är uttryckt i vektorform och beskriver kraften mellan kropparna med massan m i och m j , där kropparna benner sig på respektive position r i och r j . Vi vill enbart studera kraft och acceleration i två dimensioner och kan därmed förenkla uttrycket i ekvation (1). Vidare behöver vi uttrycka acceleration för varje koordinataxel var för sig för att kunna använda Velocity-Verlet-metoden. Viktigt är att ställa upp uttryck utan sinus och cosinus som annars senare komemr ställa till teckenbesvär.

• Ställ upp ett uttryck för avstådet och kartersiska koordinater.

r ij mellan kropp i och j i två dimensioner • Uttryck gravitationskraften mellan kropp med massan m j i med massan m i och kropp i två dimensioner och kartersiska koordinater.

j • Uttryck accelerationen mellankropp massan m j i med massan m i och kropp i två dimensioner och kartersiska koordinater.

j med

2.2 En satellit i omloppsbana

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 nämnvärt kommer påverkas 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.

7

y r F x F m F y M θ x Figur 3: 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 i om loppsbana kring en planet. Rörelseekvationerna ska integreras med hjälp av Velocity-Verlet. Funktionen ska använda följande funktionshuvud: function [x,y,vx,vy,ax,ay,t]=orbit_1body(G,m,x0,y0,vx0,vy0,dt,tmax) där G är gravitationskonstanten, m är massan hos planeten, x0,y0 och vx0,vy0 beskriver satellitens position och hastighet vid tiden t=0. Varia beln dt är längden på tidsteget och tmax är den totala simuleringstiden.

Vi kommer alltså att simulera systemet från t=0 till t=tmax. Notera att satellitens bana är helt oberoende av dess massa. Funktionen returnerar sju vektorer x,y som innehåller satellitens koordinater med en datapunkt för varje tidsteg i simuleringen, vx,vy som innehåller satellitens hastighet för varje tidsteg, ax,ay som innehåller satellitens acceleration vid varje tidssteg, samt t som innehåller tiden för varje datapunkt. På sista sidan i instruktionen nns ett kodexempel. Använd gärna detta som en utgångs punkt.

• 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.

Rita upp satellitens bana, hastighet och acceleration. Markera även planets position i guren. För att få rätt proportioner på guren, se till att använda axis equal eller axis image. Hastigheten och accelerationen kan med fördel plottas med hjälp av quiver.

Simulera med olika värden på dt och se hur noggrannheten i simule ringen förändras. Vilket värde på dt verkar lämpligt att använda?

Kontrollera att energin bevaras i simuleringen genom att plotta den kinetiska energin E k , den potentiella energin Kan ni hitta all energi i systemet, så att E p E + p E samt k = 0 E ?

p + E k och se hur de förändras med tiden. Stämmer det med vad ni förväntar er?

8

Kontrollera att rörelsemängden bevaras i systemet. Är met.

p = 0 ? Får ni det förväntade resultatet? Om inte förklara vad som händer i syste Räkna ut masscentrum för systemet och plotta det i samma gur som ritar upp satellitens bana. Står masscentrum still?

Beräkna omloppstiden

Vi ska nu beräkna omloppstiden för satelliten i föregående uppgift. Vi kommer att behöva räkna ut er omloppstider i de följande uppgifterna och vi ska därför automatisera den här 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. Tips: Börja med att plotta för att beräkna omloppstiden.

x och/eller y mot tiden för att avgöra vilken egenskap hos kurvorna som kan användas Testkör funktionen. Använd data med samma initialvillkor som i fö regående uppgift. Vilken omloppstid har satelliten? Är omloppstiden konstant? Om inte, fundera över om ett korrekt tidssteg används.

Hur påverkas omloppstiden om satellitens initialhastighet ökar? Vid vilken hastighet slutar satelliten att gå i omloppsbana runt planeten?

Simulera med samma initialvillkor som tidigare men öka stegvis vär det på vy0. Vilket är förhållandet mellan den kinetiska energin E k och den potentiella energin E p när satelliten lämnat omloppsbanan?

Rymdstationen ISS

Rymdstationen ISS (International Space Station) ligger i en nästan cirkulär om loppsbana kring jorden. Omloppsbanan är en så kallad LEO (Low Earth Orbit) vilket innebär att stationens omloppsbana 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 och anpassa tidssteget dt så att ni får en stabil lösning. Vilken omloppstid har stationen? Hur väl stämmer simuleringen överens med verkligheten?

2.3 Jorden och månen

Nu ska vi utöka simuleringen till att även inkludera planetens rörelse. Eftersom vi simulerar rörelsen hos båda kropparna behöver vi inte göra några antaganden om kropparnas massa, till skillnad från i förra uppgiften då vi antog att satelliten vägde mycket mindre än planeten. Vårt system kommer nu att 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 den tidigare simule ringen. Funktionen ska använda följande funktionshuvud: function [x,y,vx,vy,ax,ay,t]=orbit_2body(G,m,x0,y0,vx0,vy0,dt,tmax) 9

Skillnaden mot tidigare är att vi nu har två kroppar i rörelse, vilket kräver dubbelt så många initialvillkor. Variablerna m, x0, y0, vx0 och vy0 kommer därför att vara vektorer med två element vardera. Den initiala hos båda kropparna och har dimensionen 2 × x -positionen för den första kroppen anges till exempel i det första elementet x0(1) och motsvarande data för den andra kroppen anges i x0(2) osv. Vektorerna x, y, vx, vy, ax och ay som funktionen returnerar innehåller data för rörelsen 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 desamma som i uppgift 2 och placeras på position 2 i vektorerna, t.ex. m(2)=0.01.

Rita upp satellitens och planetens bana, hastighet och acceleration i samma gur och jämför med resultatet från uppgift 1.

Undersök planetens rörelse. Har den en sluten bana, dvs kommer den tillbaka till samma punkt efter ett varv? Vad händer med masscent rum för systemet? Rita ut masscentrum för systemet i samma gur som satellitens och planetens bana. Testa att öka massan på satelliten till m(2)=1 för att få en tydligare eekt.

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?

Kontrollera att energin är bevarad.

Kontrollera att rörelsemängden är bevarad.

Månens omloppstid runt jorden

Månen är jordens enda naturliga satellit. Det är den femte största månen i vårt solsystem och den rör sig i en elliptisk bana runt jorden. Med en periodtid runt jorden på en månad har den varit en naturlig referens för tidiga kalendrar.

I relation till 'stillastående' stjärnor har månen en omloppstid på 27.3 dagar runt jorden, men eftersom jorden i sig själv roterar runt solen tar det 29.5

dagar från en fullmåne till nästa. Använd din förnade kod för att räkna ut månens omlopptid runt jorden. Hitta själv initialvilkor från lämplig källa. Hur väl stämmer periodtiden överens med verkligheten?

2.4 Solsystemet

Nu är vi redo att simulera solsystemet. Vi ska utöka simuleringen så att den kan hantera rörelsen hos hantera N N planeter som alla påverkar varandra genom gravitation.

Beroende på hur du valde att implementera lösningen till uppgift 3 kan det ha blivit många index att hålla reda på. När vi nu ska utöka simuleringen till att planeter behöver vi generalisera koden.

• Skriv en funktion acceleration som beräknar accelerationen för alla N planeterna. Om du följt kodexemplet på baksidan av dessa instruktioner 10

har du förmodligen redan implementerat en funktion som räknar ut acce lerationen, fast då som en anonym funktion. Nu behöver du utöka denna funktion till mer än en rad. Funktionen ska använda följande funktions huvud function [ax,ay]=acceleration(G,m,x,y) där ax, ay, m, x och y är vektorer med längd planeternas massor, x och y är planeternas positioner, G är gravitations konstanten och ax, ax är den totala accelerationen som verkar på varje planet som räknas ut genom att summera accelerationen från alla andra planeter.

N . Parametern m innehåller • Skriv en funktion orbit_Nbody som simulerar rörelsen hos N st himla kroppar med hjälp av Velocity-Verlet-metoden. Utgå från er tidigare simu lering och använd funktionen acceleration för att beräkna accelerationen mellan planeterna. Funktionen ska använda följande funktionshuvud: function [x,y,vx,vy,ax,ay,t]=orbit_Nbody(G,m,x0,y0,vx0,vy0,dt,tmax) Inparametrarna är samma som tidigare med skillnaden att massan och initialvärdena nu är vektorer med längd sionen N × N eftersom vi nu har N kroppar som behöver initieras. Av samma anledning har nu returvärdena dimen steps för att kunna inrymma informationen om 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.

• Simulera solsystemet! Simulera solen och de inre planeterna Merkurius, Venus, Jorden och Mars. Detta blir alltså en simulering med N = 5 .

Använd lämplig källa för att hitta initialvärden för planeterna. För att förenkla valet av initialvärden kan vi anta att alla planeterna ligger längs en rät linje 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 ellip tiska. (För att få korrekta elliptiska banor behöver vi veta avståndet till solen och varje planets hastighet i en specik punkt längs sin bana 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 rörelse.

Gör en gur som visar planeternas banor.

Kontrollera att energi och rörelsemängd bevaras.

Beräkna omloppstiderna och jämför med tabellvärden.

11

A 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, med antagandet att föremålets position och hastighet vid t = 0 är kända. Vi utgår från Newtons andra lag, d 2 x dt 2 = a ( x ( t )) = F ( x ( t )) , m (A.1) där F ( x ( t )) är en konservativ kraft som endast beror på positionen kretiserar tiden och approximerar andraderivatan enligt följande: x . Vi dis a ( x ) = d 2 x dt 2 ≈ x n +1 − x n ∆ t − x n − x n − 1 ∆ t ∆ t = x n +1 − 2 x n + x n − 1 ∆ t 2 = a n .

Om vi löser ut x n +1 får vi (A.2) x n +1 = 2 x n − x n − 1 + a n ∆ t 2 (A.3) vilket är en fullt fungerande integrationsmetod som kan användas för att räkna ut framtida positioner en given position x n x n omgärdar den i tid, enligt +1 givet de två föregående positionerna x n och x n − 1 samt accelerationen. Metoden har utvecklats era gånger av olika forskare ge nom historien. Den kallas oftast Verlets metod efter den senaste upphovsmanen 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 vid en given tidpunkt för att räkna ut nästa position och använder inte den momentana has tigheten. 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 hastig heten. Men om vi behöver hastigheten blir det lite omständligt. Hastigheten i går att räkna ut från skillnaden mellan två positioner som v n = x n +1 − x n − 1 , 2∆ t (A.4) men vi ser att för att räkna ut denna hastigheten nästkomande positionen t = 0 x n +1 , vilket inte alltid är praktiskt. En annan sak att tänka på är att man oftast bara känner till position , så för att komma igång måste x n − 1 v n måste vi redan känna till x först beräknas med hjälp av någon annan metod. För att Verlet-metoden ska fungera måste 0 , alltså positionen vid ∆ 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. Meto den 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 att lösa ut x n − 1 från ekvation (A.4). Vi får då x n − 1 = x n +1 − 2 v n ∆ t.

(A.5) 12

Vi sätter in uttrycket för x n − 1 i ekvation (A.3) vilket ger x n +1 = x n + v n ∆ t + 1 2 a n ∆ t 2 , (A.6) som är en ny uppdateringsmetod 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 (A.4) kan vi sätta upp ett uttryck för v n +1 : v n +1 = x n +2 − x n .

2∆ t (A.7) Vi kan ta fram uttryck för ekvation (A.7) vilket ger x n +2 och x n från ekvation (A.6) och sätta in dem i v n +1 = x n +1 + v n +1 ∆ t + 1 2 a n +1 ∆ t 2 2∆ t − vilket efter förenkling blir x n +1 − v n ∆ t − 1 2 a n ∆ t 2 , v n +1 = v n + 1 2 Vi har nu två uppdateringsekvationer: a n + a n +1 ∆ t.

(A.8) (A.9) x v n n +1 +1 = = x n v n + v n ∆ t + 1 2 a n + 1 2 a n ∆ + a n +1 t 2 ∆ t (A.10) som tillsammans utgör metoden vi kallar Velocity-Verlet. Metoden är självstar tande, dvs vi behöver inte använda en separat metod för att beräkna Även Velocity-Verlet kräver ett konstant F = −∇ V egenskaper.

( r luftmotstånd.

) , där r ∆ t är positionsvektorn, och V x n − 1 ef tersom det värdet aldrig används till skillnad från i den vanliga Verlet-metoden.

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, är en funktion som beskri ver den potentiella energin i systemet. Den kraft vi ska simulera är konservativ 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 Slutligen kan nämnas att det nns en annan populär metod som kallas Le apfrog som också går att härleda genom att göra en omskrivning av Verlet metoden. Velocity-Verlet och Leapfrog är väldigt lika och har i princip samma 13

Exempelkod

Nedan återges ett exempel på strukturen i ett program som använder stegnings metoder liknande Velocity-Verlet. Koden är tänkt att simulera en planet med rörelse i två dimensioner x och y kapade och avslutas istället med . 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 ...

.

%starting point for a simple simulation program %preallocate memory (increases performance) steps=... %select number of time steps x=nan(steps,1); y=nan(steps,1); vx=nan(steps,1); vy=nan(steps,1); ax=nan(steps,1); ay=nan(steps,1); %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