Jämförelse av implicita numeriska metoder för en styv Van der Pol-ekvation

(1)

2

INOM

EXAMENSARBETE DATATEKNIK, GRUNDNIVÅ, 15 HP

STOCKHOLM SVERIGE 2017,

Jämförelse av implicita numeriska metoder för en styv Van der Pol- ekvation

ANDREAS BRYNOLFSSON BORG

KTH

SKOLAN FÖR DATAVETENSKAP OCH KOMMUNIKATION

(2)

Jämförelse av implicita numeriska

metoder för en styv Van der Pol-ekvation

ANDREAS BRYNOLFSSON BORG

Civilingenjör Datateknik Datum: 5 juni 2017

Handledare: Michael Schliephake Examinator: Örjan Ekeberg

Engelsk titel: Comparison of Implicit Methods for a Stiff Van der Pol Oscillator

Skolan för Datavetenskap och Kommunikation

(3)

ii

Sammanfattning

Att på ett så effektivt sätt som möjligt kunna lösa styva ordinära differentialekvationer med godtycklig noggrannhet är ett fortsatt levande område inom matematiken. Denna studie hade som mål att jämföra tre implicita numeriska metoder lämpade att lösa en styv Van der Pol-ekvation, samt undersöka till vilken grad vektorisering i MATLAB kan användas för att förbättra prestanda för sådana algoritmer. Resultaten visar att MATLAB:s inbyggda metod ODE15s är snabbast och ger god noggrannhet då finkänsliga svar krävs, medan den implicita mittpunktsmetoden genererar mycket goda lösningar till problemet även i fall där storleken på felet inte är centralt. Dock är denna metod mer resurskrävande än de övriga. Dessutom visar sig metoden RADAU5 vara ett bra metodval vid lösning av styva ODE i fall är toleransen är större än 10⁻³, även om ODE:n har stor kurvatur. Effekten av vektorisering vid lösning av ett tvådimensionellt ODE-system uppskattas ge möjlighet till en halvering av kör- tiden, jämfört med en ovektoriserad likvärdig funktion.

(4)

iii

Abstract

Finding as efficient numerical algorithms for solving stiff ordinary differential equations as possible is still an ongoing pursuit in mathematics and computer science. The purpose of this study was to compare three different numerical methods particularly well-suited for solving a stiff Van der Pol oscillator. In addition to this, the effects on performance by vectorizing such methods in MATLAB were studied. The results of the study show that MATLAB’s ODE- solver ODE15s is the fastest method and provides well-approximated solutions when the size of the error is important, whereas the implicit midpoint method generates very exact solutions regardless of the tolerance used. That method is substantially more computationally heavy however. In cases where greater tolerances can be used (greater than 10⁻³), RADAU5 can be used, even if the ODE has significant curvature. The effects of vectorization when solving a 2- dimensional system of ODEs can provide an approximate reduction in runtime by half compared to a non-vectorized version of the same function.

(5)

(6)

Innehåll

Innehåll v

1 Introduktion 1

1.1 Syfte . . . 2

1.2 Avgränsningar . . . 2

1.3 Frågeställning . . . 2

1.4 Rapportens utformning . . . 2

2 Bakgrund 3 2.1 Kort om differentialekvationer . . . 3

2.2 Van der Pol-ekvationen . . . 4

2.3 Numeriska metoder för ODE:s . . . 5

2.4 Styvhet . . . 7

2.5 Stabilitet . . . 8

2.6 Vektorisering i MATLAB . . . 8

2.7 Studerade numeriska metoder . . . 9

2.7.1 RADAU5 . . . 9

2.7.2 ODE15s . . . 10

2.7.3 Implicit mittpunktsmetod . . . 10

3 Metod 12 3.1 Noggrannhet . . . 12

3.1.1 Exakt instans . . . 12

3.1.2 Styv Van der Pol . . . 13

3.2 Tidsåtgång . . . 13

4 Resultat 15 4.1 Noggrannhet . . . 15

4.2 Tidsåtgång . . . 19

v

(7)

vi INNEHÅLL

5 Diskussion 22

5.1 Noggrannhet . . . 22

5.2 Tidsåtgång och vektorisering . . . 24

6 Slutsats 25

Litteraturförteckning 26

(8)

Kapitel 1

Introduktion

Sedan länge har människan haft nytta av att kunna förutspå och förutse di- verse naturliga fenomen: Från att tyda solens rörelser för att så vid rätt tid i hopp om att maximera skördar, till att förstå komplexa samband inom biolo- gi, fysik och ekonomi. För att på ett strukturerat sätt kunna modellera och lösa problem har matematiken ofta använts. I takt med matematikens utveckling har fler problem kunnat lösas på allt mer sofistikerade vis. De senaste sekler- na har differentialekvationer visat sig vara praktiska för att kunna göra goda uppskattningar om verkligheten. Idealt vore att kunna hitta exakta analytiska lösningar till samtliga differentialekvationer, men detta har visat sig vara matematiskt omöjligt (en konsekvens av att vissa funktioner saknar elementära primitiva funktioner). Istället får man förlita sig på approximationer till lös- ningar av differentialekvationer med hjälp av så kallade numeriska metoder. Att beräkna approximerade lösningar för hand är både opraktiskt och tidsödande, varför man gärna tar datorer till hjälp. Att finna algoritmer som är både snabba och ger goda numeriska approximationer till differentialekvationer är ett ämne som har studerats sedan 1970-talet och är än idag ett viktigt forsknings- område [1].

En viss typ av differentialekvationer karaktäriseras av störningar i de approximerade lösningarna som numeriska metoder genererar. Detta gäller i synnerhet när explicita numeriska metoder används vid ekvationslösningen. Differen- tialekvationer som tillhör denna klass kallas för styva, och löses därför hellre med implicita numeriska metoder. Problemet med implicita numeriska metoder ligger i faktumet att de allt som oftast är mer resurskrävande än explicita.

Därför kan avvägningar mellan effektivitet och noggrannhet behöva göras.

En välkänd styv differentialekvation är Van der Pol-ekvationen. Den har både historisk signifikans och nutida relevans, och används för att modellera saker som elektriska kretsar till ekonomiska cykler [2].

1

(9)

2 KAPITEL 1. INTRODUKTION

1.1 Syfte

Eftersom styva differentialekvationer har visat sig dyka upp i många viktiga områden, såväl teoretiska som praktiska, är det av vikt att finna effektiva och tillförlitliga algoritmer som kan lösa dem [3]. Studiens syfte är således att un- dersöka och jämföra tre numeriska metoder särskilt lämpade för att lösa en styv instans av Van der Pol-ekvationen. Centrala mätvärden att titta på i sammanhanget är tidsåtgång och numerisk noggrannhet. I utvecklingsmiljön MAT- LAB kan man dessutom skriva numeriska algoritmer på ett vektoriserat sätt.

Detta programmeringsparadigm påstås möjliggöra avsevärda prestandaförbätt- ringar och att åskådliggöra effekterna av det är intressant i sammanhanget.

1.2 Avgränsningar

Rapporten är avgränsad till att enbart studera en styv instans av Van der Pol- ekvationen samt att enbart jämföra tre numeriska metoder för att approximera en lösning till denna. Trots att det kan vara intressant att undersöka prestandan av explicita metoder för styva differentialekvationer görs inte detta i rapporten. Istället studeras de implicita metoderna: RADAU5, MATLAB:s ODE15s och den implicita mittpunktsmetoden. De variabler som mäts och analyseras är tidsåtgång för att lösa problemet, samt den numeriska noggrannheten hos de genererade lösningarna. Programmet MATLAB används för samtliga beräk- ningar eftersom det är utformat att effektivt hantera numerisk kalkyl genom så kallad vektorisering av funktioner.

1.3 Frågeställning

Hur skiljer sig RADAU5, ODE15s och den implicita mittpunktsmetoden, vid lösning av en styv Van der Pol-ekvation med avseende på tid och noggrannhet? Vidare, till vilken grad kan prestandan av den implicita mittpunktsmetoden förändras då den implementeras som en vektoriserad funktion för samma Van der Pol-ekvation?

1.4 Rapportens utformning

Det nästkommande kapitlet är avsett att ge läsaren en förståelse för teorin och koncepten bakom differentialekvationer, fenomenet styvhet, Van der Pol- ekvationen och en överblick av numeriska metoder. Dessutom presenteras tidigare arbete inom området, såväl som de tre algoritmer som studeras i rapporten. Kapitel 3 behandlar metodologin för hur mätdata har samlats in och behandlats. Resultaten från studien presenteras i kapitel 4 och diskuteras senare i kapitel 5. Slutsatser om studien redogörs för i kapitel 6.

(10)

Kapitel 2

Bakgrund

Nedan beskrivs den matematiska teorin som ligger till grund för studien. De använda numeriska metoderna och konceptet vektorisering beskrivs dessutom i det omfång som krävs för att till fullo förstå senare kapitel i rapporten.

2.1 Kort om differentialekvationer

I samband med upptäckten av differentialkalkylen av Newton och Liebniz på 1600-talet uppstod även behovet av att kunna lösa differentialekvationer. Myc- ket inom området har sedan dess utvecklats och differentialekvationer är idag en fullvärdig matematisk gren [2]. Kort sagt är differentialekvationer ekvationer som relaterar en funktion med dess derivator. Studien om dem går ut på att hitta sätt att lösa ut funktionen. Ett känt introduktionsexempel är

y⁰(x) = y(x) (2.1)

som har den exakta lösningen y(x) = Ce^x, för någon konstant C. Ekvation (2.1) hör till en särskild typ av differentialekvationer vid namn ordinära differentialekvationer (ODE:s). En ODE är en differentialekvation där funktionen i fråga, samt dess derivator, endast beror på en oberoende variabel. I fall där det tydligt framgår vad den oberoende variabeln är brukar man ofta underlätta läsandet genom att inte skriva ut funktionsargumenten [2]. Exempelvis skulle ekvation (2.1) ovan förkortas till den ekvivalenta differentialekvationen

y⁰ = y. (2.2)

Ett begynnelsevärdeproblem i studiens sammanhang betecknar en ODE med en mängd givna funktionsvärden utvärderade i någon begynnelsepunkt t0. Nota- tionen för detta är vanligtvis

y(t0) = y0. (2.3)

3

(11)

4 KAPITEL 2. BAKGRUND

Varje ODE tillskrivs dessutom en så kallad ordning. Ordningen av en ODE är ett heltal som helt enkelt är lika med den högsta derivatan av funktionen som finns i differentialekvationen (således är ekvationerna ovan av första ordningen). Många metoder ägnade åt att lösa ODE:s fungerar enbart för ekvationer av första ordningen. En vanlig teknik för att lösa n:te ordningens ODE:s är därför att skriva om den som ett system av n första ordningens ODE:s genom att införa nya variabler och ställa upp ett ekvationssystem [4].

Notationen och teorin för system av första ordningens ODE:s är en utökning av fallet med bara en ODE, där funktionerna dock är vektorvärda istället för skalära. Standardformen för ett ODE-system av första ordningen med begynnel- sevillkor skrivs generellt som

du

dt = f (t,u), u(t0) =u0. (2.4)

2.2 Van der Pol-ekvationen

Den här studien ägnas åt att hitta lösningar till Van der Pol-ekvationen. Ek- vationen formulerades först av Balthasar van der Pol på 1920-talet då han ar- betade med oscillationer i elektronrör. Den har dock visat sig vara användbar i många andra sammanhang där så kallade gränscykler (en slags periodisk at- traktor) uppstår [5]. Ekvationen skrivs ofta på ett av två vis. Antingen presenteras den som en ODE av andra ordningen

y⁰⁰+ µ(y²− 1)y⁰+ y = 0, µ ∈ R > 0, (2.5) eller som ett tvådimensionellt system av första ordningens ODE:s, genom att införa beteckningarna y1 = y⁰ och y2 = y₁⁰. Detta resulterar i det följande systemet

y⁰₁= y2

y⁰₂= µ(1 − y²₁)y2− y₁, µ ∈ R > 0.

(2.6)

Van der Pol-ekvationens lösning är periodisk och det har visat sig att perioden ökar i takt med µ [3]. För att lättare studera lösningen på små intervall är det därför fördelaktigt att skala ner perioden genom att manipulera ODE-systemet.

Hairer [3] har visat att det är möjligt att transformera systemet ovan till ett system med lika stor amplitud men med en kortare period, och skrivs på formen

y₁⁰= y₂

y₂⁰= (1 − y²₁)y₂− y₁/, ∈ R > 0.

(2.7)

(12)

KAPITEL 2. BAKGRUND 5

En annan fördel med att reducera perioden på Van der Pol-ekvationen är det kan gå fortare att approximera lösningar till den utan att förändra den egentliga strukturen av problemet. Den transformation som används i denna studie specificeras i [6] och gör det möjligt att lösa problemet på ett intervall tre stor- leksordningar mindre än det ursprungliga.

Ett Van der Pol-problem av samma format som ekvation (2.7) används i denna studie. Systemets begynnelsevärden är y1(0) = 2, y2(0) = −0.66och = 10⁻⁶. Denna uppsättning parametrar ger upphov till ett styvt ODE-system [6]. På vektorform införs beteckningarna

y =



 y1

y2



, y⁰=



 y⁰₁ y⁰₂



, f (t,y) =





y2

10⁶ (1 − y₁²)y2− y₁



. (2.8)

Standardformen för den styva Van der Pol-ekvationen som undersöks i den här studien blir således

y⁰ = f (t,y), y(0) =



 2

−0.66



. (2.9)

2.3 Numeriska metoder för ODE:s

Den grundläggande idén bakom numerisk lösning av ODE:s är att använda tidigare vetskap om funktionen för att extrapolera om dess framtida beteende.

Ofta försöker man göra detta genom att först dela upp (diskretisera) intervallet man är intresserad av att finna lösningen på. Därefter stegar man sig igenom intervallet med någon ändlig steglängd och uppskattar funktionens lutning i olika punkter med hjälp av gamla funktionsvärden och differenskvoter. Lin- järinterpolation mellan de utvärderade lösningspunkterna används därefter för att få en kontinuerlig funktion [2]. Traditionellt har fix steglängder använts, men med hjälp av datorkraft kan man idag utnyttja adaptiva steglängder som varierar beroende på funktionens kurvatur och lutning. Fördelen med metoder som kan nyttja adaptiva steglängder är att de kan bespara beräkningar på intervall där funktionen har ett förhållandevis långsamt förändringsförlopp. På intervall där storleken av förstaderivatan är liten i beloppet kan man då undvi- ka att utföra mer eller mindre identiska beräkningar upprepade gånger [3].

De algoritmer som har utvecklats över åren för att lösa ODE:s tillhör alltsom oftast en av två klasser: explicita och implicita metoder [4]. Explicita metoder är metoder där differentialekvationens nästa steg kan skrivas som en funktion av tidigare funktionsvärden och steglängder. Det går alltså att analytiskt lö- sa ut nästa funktionsvärde och lösningsprocessen reduceras till att substituera värden för att beräkna nya. Implicita metoder å andra sidan sammanfogar det

(13)

nästkommande funktionsvärdet med tidigare data och kräver på så sätt ekva- tionslösning för att få fram nya värden [7]. Fördelar med implicita metoder är att steglängder kan tas grövre än för explicita metoder, samt att metodernas stabilitetsregioner generellt sett är större [8]. Nackdelen med implicita metoder är att de ofta är beräkningstyngre än explicita i och med den ekvationslösning som krävs för att extrahera framtida funktionsvärden. Detta gäller i synnerhet för system av ickelinjära ekvationer som inte kan lösas med hjälp av snabb Gausseliminering.

Ekvationslösning av ickelinjära ekvationssystem kan utföras på flera olika sätt, med godtycklig noggrannhet. I algoritmerna som används i den här studien används Newton-Raphsons metod. Den utnyttjar den så kallade jakobianmatrisen som består av partialderivatorna till ekvationssystemet. För ett allmänt ekvationssystem med ekvationer f1, . . . , f_m med variabler x1, . . . , x_n betecknas jakobianmatrisen

J =







∂f1

∂x1 . . . _∂x^∂f¹

n

... . .. ...

∂fm

∂x1 . . . ^∂f_∂x^m

n







. (2.10)

I fallet för ekvation (2.8) blir jakobianen

J =





0 1

10⁶(−2y₁y₂− 1) 10⁶(1 − y₁²)



. (2.11)

Newton-Raphsons metod arbetar sedan iterativt fram en lösning med kvadra- tisk konvergens. I korta drag kan metoden specificeras som följande.

NewtRaph(ODE på formen f(x) = 0, jakobian J(x), rotgissning x, tolerans ) while kxk∞>

Beräkna en korrektionsterm xδ = −J(x)⁻¹f(x) Uppdatera vektorn x = x + xδ

end while Returnera x

Algoritm 1: Generaliserad Newton-Raphsons metod för flerdimensionella system.

Den returnerade vektorn är således en kandidat för nästa diskretiserade funk- tionsvärde.

(14)

2.4 Styvhet

Trots allt arbete och den forskning som ägnats åt studien av styva ODE:s finns det ingen vedertagen matematisk definition av termen. Curtiss [9], som först upptäckte styva ODE:s år 1952, observerade dem som “/. . ./ equations where certain implicit methods, and in particular backward-differentiation formulae (BDFs), perform better, usually tremendously better, than explicit ones” när försöka att lösa dem numeriskt görs.

En annan definition av styva ODE:s presenterad av Hairer [6] är: “Stiff problems are characterized by the fact that the numerical solution of slow smooth movements is considerably perturbed by nearby rapid solutions”.

Oavsett ordval handlar det i vilket fall om en svårklassad komplexitet hos dessa ODE:s. I rapportens sammanhang är en tillräcklig definition av styva ODE:s att explicita metoder inte går att använda för att lösa dem inom rimlig tid med hyfsad noggrannhet. Vidare sägs ODE:s som inte är styva vara ostyva.

De två stora problem kopplade till styva ODE:s avspeglas i denna definition.

Det första är att de genererade lösningarna kan vara dåliga approximationer. Ett typexempel med den explicita framåteulermetoden visas nedan i figur 2.1. Precis enligt definitionen ovan fungerar det inte att lösa den styva ODE:n y⁰ = −50(y − cos(t)), y(0) = 0på ett explicit vis. Det är helt enkelt i vissa fall svårt att förutse funktionens beteende och kraftiga oscillationer uppstår. In- stabiliteten hos dessa svängningar tenderar heller inte att avta då t → ∞ [6].

Implikationerna av detta är att lösningar kan vara opålitliga och inte vidare användbara.

Figur 2.1: Exempel på numerisk instabilitet hos explicita ODE-lösare.

(15)

Det andra huvudsakliga problemet med styva ODE:s är att det kan ta oskäligt mycket tid att finna lösningar till dem. Den extra tidsåtgången är i viss mån kopplad till det första problemet, eftersom felet på lösningens noggrannhet kan förbättras genom att dra ner på steglängden, vilket resulterar i fler beräk- ningar på integrationsintervallet [4]. Som ovannämnt är detta starkast förknip- pat med explicita metoder men även implicita metoder, oavsett om de brukar adaptiv steglängd eller inte, kan kräva stor processorkraft i förhållande till problemets storlek.

2.5 Stabilitet

Ett koncept starkt kopplat till styvhet är stabilitet. Både ODE:s och numeriska metoder kan vara stabila, om än i lite olika bemärkelser. Enbart numeriska metoders stabilitet behandlas i denna rapport. I korta drag är ett mått på hur pass stabil en numerisk metod är för att lösa styva problem dess stabilitetsregion. Regionerna beskrivs som delmängder till C, det komplexa talplanet.

Något förenklat kan man säga att en metod är stabil i en region om lösningen till en styv ODE inte divergerar exponentiellt ifall begynnelsevärdena skulle ändras något [3]. Det vill säga, likadana ODE:s med störda begynnelsevärden u(t0± δ) = uσ kommer aldrig ha lösningar som obundet avviker från varandra.

En konsekvens av att använda ODE-lösare med högre noggrannhetsordningar är att stabilitetsregionen minskar [7]. På grund av detta fenomen krävs avväg- ningar vid val av klasser av ODE-lösare med olika noggrannhetsordningar.

Metoder vars stabilitetsregion minst täcker C⁻, mängden av alla komplexa tal med negativ realdel, sägs vara A-stabila [3].

2.6 Vektorisering i MATLAB

Kärnan av MATLAB som programspråk och utvecklingsmiljö utgörs av matriser och operationer på dessa (därav namnet MATrix LABoratory). Stor vikt har lagts åt att optimera de kompilerade funktionerna som medföljer programmet.

Eftersom i stort sett alla dessa funktioner har stöd för vektorvärda indata och utdata finns det stora fördelar med att utnyttja vektorer i MATLAB [10].

Det som ligger till grund för detta beteende är att kostnaden (i tid och processorkraft) för funktionsanrop är relativt sett stort. Däremot ökar kostnaden för att utvärdera en funktion långsamt i förhållande till antalet funktionsar- gument funktionen tar [11]. Detta gör det möjligt att operera på flera värden

“simultant” med samma operator på ett resursbesparande sätt (exempelvis kan loopar, därmed upprepade funktionsanrop, undvikas). Matriser kan även reducera antalet cachemissar och på så sätt förbättra körtider.

(16)

Avseende minneshantering så lagras data i matriser kolumnvis i MATLAB. Det innebär att matrisen

A =







a11 a12 a13

a21 a22 a23

a31 a32 a33







i datorns minne representeras som A = [a11, a₂₁, a₃₁, a₁₂, a₂₂, a₃₂, a₁₃, a₂₃, a₃₃].

Att hämta data som är närliggande i minnet är fördelaktigt för prestanda eftersom cachens snabba accesstider utnyttjas maximalt [10]. Hur data i MATLAB- program hämtas och lagras är således viktigt. Kolumnvis är att föredra, varför man hellre använder den vektoriserade wildcard-syntaxen A(:,i) (samtliga element i kolumn i) istället för A(i,:) (alla element på rad i) om möjligheten finns. Lyckligtvis är ju ekvation (2.8) tvådimensionell så den lämpar sig väl för kolumnbaserad vektorisering, samt för MATLAB överlag. Samtliga studerade vektoriserade algoritmer använder sig därför av kolumnvisa operationer.

Vektorisering realiseras även i hårdvaran. Nyare versioner av MATLAB (från och med R2013b) utför vektorisering med hjälp av AVX-instruktioner som är SIMD (same instruction-multiple data) utan behov av särskild multitrådig kod [12]. Stora delar av MATLAB:s linjär algebra-bibliotek består av just detta, vilket har medfört prestandaförbättringar av MATLAB-kod på senare år [11].

Genom att utnyttja konceptet vektorisering i den mån det är möjligt kan man kraftigt reducera resurser som krävs för att köra ett program; tidigare arbe- ten av bland andra Shampine [13] har visat att exekveringstider kan förbättras med åtminstone en storleksordning.

2.7 Studerade numeriska metoder

2.7.1 RADAU5

Den historiska utveckling som ligger till grund för metoden RADAU5 är lång och invecklad. För en längre och mer detaljerad redogörelse för metoden hän- visas det till Hairer [6], som skapade algoritmen. Metoden är som ovannämnt implicit och bygger på så kallad Radau-kvadratur som härstammar från sent

1800-tal, Runge-Kutta metoder, samt femte ordningens BDF (backwards-differentiation formula). Metoden är alltså enligt Curtiss [9] definition av styvhet väl lämpad

att lösa ekvation (2.8).

Implementationen som används i denna rapport är en starkt vektoriserad MATLAB- översättning av originalet skrivet i Fortran [14] [15]. Metoden använder sig av en adaptiv steglängd av effektivitetsskäl och Newton-Raphsons metod för lös- ning av ickelinjära ekvationssystem. Metoden är A-stabil [6], trots den relativt höga ordningen av BDF som används.

(17)

2.7.2 ODE15s

MATLAB har en uppsättning inbyggda mycket användbara ODE-lösare vars prestanda och användningsområden varierar. Just ODE15s är specialanpassad för att hantera styva ODE:s inom rimlig tid. Förutom variabel steglängd kan algoritmen dessutom köras med en samling olika inställningar. En av dessa in- ställningar är användning av BDF av andra ordningen [16], och är av intresse utifrån Curtiss definition.

En annan taktik ODE15s använder sig av för att hantera styva problem är att algoritmen utnyttjar variabel ordning för att påskynda lösningsprocessen. Att uttala sig om metodens stabilitetsregion är därför svårt eftersom ordningen kan variera från 1 till 13, steg emellan. Även här körs Newton-Raphsons metod för att approximera framtida funktionsvärden. ODE15s skiljer sig från de två övriga algoritmerna i det att den kan beräkna jakobianmatriser numeriskt med differenskvoter istället för att ta den analytiskt utvärderade matrisen som en del av dess input. Metoden gör detta på ett resurssnålt sätt, genom att utnyttja snabba och robusta linjär algebra-kommandon för återanvända tidigare utvärderade jakobianer där möjlighet ges [16].

Eftersom ODE15s är en rutin som medföljer MATLAB är rutinen som körs skriven i mer maskinnära språk (Fortran och C [11]) än ordinära MATLAB- skript [4]. Implikationerna av detta är att algoritmen i teorin borde vara snabbare än övriga studerade metoder som annars behöver förbehandlas och tol- kas. Koden som faktiskt exekveras på hårdvaran är dock inte tillgänglig för gängse MATLAB-användare, så exakta detaljer om den går inte att uttala sig kring.

2.7.3 Implicit mittpunktsmetod

Den implicita mittpunktsmetoden är också besläktad med metoderna framtag- na av Runge och Kutta [17]. Metoden är den simplaste av de tre som under- söks i denna rapport. Rent matematiskt är den implicita mittpunktsmetoden definierad som

y_n+1=y_n+ hf

tn+h

2,1

2(y_n+y_n+1)

.

Geometriskt kan man tolka det som att metoden försöker approximera lut- ningen av funktionens tangent i mittpunkten av intervallet y_n till y_n+1. För att lösa ut vektorn y_n+1 i det ickelinjära ekvationssystemet som uppstår används återigen Newton-Raphsons metod. I pseudokod kan metoden sammanfattas som följande.

(18)

ImpMittpunkt(ODE y’ = f(t, y), jakobian J(f), y₀, [t0, t_slut], tolerans ) U ←lösningsmatris

U[0] ← y₀ ti ← t₀ h ←

while ti ≤ t_slut ti← t_i+^h₂

U[i + 1] ← NewtRaph(f(t, y) = 0, J(f), f(ti+^h₂,U[i]), ) Justera h baserat på kU[i] − U[i + 1]k∞

end while Returnera U

Algoritm 2: Den implicita mittpunktsmetoden.

Metoden är A-stabil och är av noggrannhetsordning O(h²) [17]. Hädanefter är denna metod synonymt kallad mittpunktsmetoden. Skillnaderna mellan den vektoriserade och ovektoriserade versionen sammanfattas nedan.

Vektoriserad mittpunktsmetod: Variablerna y1 och y2 lagras kolumnvis i en tvådimensionell vektor. Detta medför att hämtning av element och lagring av nya data kan effektiviseras med MATLAB:s wildcard-syntax (:). Med hjälp av vektorisering krävs det även mindre arbete för att hitta maxnormen mellan två intilliggande utvärderade punkter för att justera h.

Ovektoriserad mittpunktsmetod: I denna version lagras y1 och y2 i separa- ta endimensionella vektorer. Implikationerna blir att MATLAB under huven behöver utföra fler operationer vid i stort sett alla centrala steg i algoritmen.

En ineffektiv MATLAB-operation är sammanslagning och storleksändring av vektorer, något som behöver göras upprepade gånger i denna version.

(19)

Kapitel 3

Metod

I detta kapitel beskrivs de tillvägagångssätt som använts för att samla in mät- data vid lösning av Van der Pol-ekvationen. Dessutom förklaras det hur denna mätdata statistiskt behandlas för att senare analyseras.

3.1 Noggrannhet

På grund av att Van der Pol-ekvationen över huvudtaget saknar exakta lös- ningar går det inte att mäta algoritmernas faktiska noggrannhetsgrad när de försöker lösa problemet [5]. Dock kan en del ansatser göras för att ändå få be- lägg för vilken eller vilka algoritmer som presterar som bäst rent numeriskt sett.

3.1.1 Exakt instans

Inledningsvis kan man se hur algoritmerna presterar vid lösning av en ostyv ODE som går att lösa analytiskt. Därefter jämförs skillnaderna mellan de approximerade lösningarna och den egentliga lösningen och detta kan ligga till grund för diskussion för hur algoritmerna presterar för godtyckliga problem [3]. Följande ODE går att lösa analytiskt

y⁰= −y cos(t), y(0) = 1 (3.1)

och har lösningen y = e^{− sin(t)}. Samtliga algoritmer fick alltså i uppgift att lösa ekvation (3.1) på intervallet 0 ≤ t ≤ 10. Hur pass fel metodernas lösningar var kan kvantifieras på olika sätt. Det globala felet [18] av en lösning är

globalt fel = |X

(exakta värdet − uppskattade värdet)| (3.2) och har ett närbesläktat noggrannhetsmått kallat medelvärdesfelet som är lika med

medelvärdesfel = globalt fel

antal uppskattade värden. (3.3)

12

(20)

KAPITEL 3. METOD 13

Medelvärdesfelet är är kanske i synnerhet intressant vid jämförelse av metoder som utnyttjar adaptiva steglängder. Eftersom alla tre metoderna inte nödvän- digtvis försöker lösa ekvationen i lika många punkter för identiska indata är det fördelaktigt att använda ett mått som i någon bemärkelse utjämnar dessa skillnader metoderna emellan. Att undersöka respektive metods maxfel (det approximerade funktionsvärde där avvikelsen är som störst från det faktiska värdet) är också centralt.

Hur dessa noggrannhetsmått påverkas av feltoleransen som används för Newton- Raphsons metod är dessutom av intresse. Eftersom de studerade metoderna är deterministiska kommer inte antalet steg varje algoritm tar variera mellan kör- ningar för samma indata. Alltså räcker det att för varje algoritm lösa ekvation (3.1) en gång för vardera tolerans: 10⁻¹, 10⁻², 10⁻⁵ och 10⁻⁷.

För att minimera avrundningsfel vid samtliga beräkningar användes MATLAB- kommandot format long så att flyttal beräknas med upp till 15 decimalers noggrannhet. Därefter extraherades de signifikanta siffror som är korrekta gentemot vardera tolerans. Dessutom användes inställningen BDF on för ODE15s.

3.1.2 Styv Van der Pol

Därefter löstes den styva Van der Pol-ekvationen som beskrevs av ekvation (2.8) med hjälp av RADAU5, ODE15s och mittpunktsmetoden. Eftersom ekvationen över huvudtaget saknar analytiska lösningar behövs noggrannheten sättas i relation till något som anses vara exakt nog. En annan av MATLAB:s ODE-lösare, ODE113, är väl lämpad att lösa ostyva problem med hög noggrannhet (toleransen sattes till 10⁻¹⁰), men den kan inte lösa styva problem effektivt. Lösningen som ODE113 genererade användes därför som en simulerad “korrekt” en, men dess körtid var inte av intresse. Samma mått på fel från avsnittet ovan användes (linjärinterpolation nyttjades i de fall då ODE113 inte hade löst ekvationen i samma punkter som övriga funktioner, för att skatta rimliga mellanliggande värden). I tidsbesparande syfte användes grövre to- leransgränser för denna instans. Dessa var 10⁻¹, 10⁻², 10⁻³ och 10⁻⁴. Antalet korrekta signifikanta siffror avspeglas av detta val. Även här använde ODE15s sig av BDF on.

3.2 Tidsåtgång

Trots att alla tre studerade metoder är deterministiska så är deras körtider i praktiken inte det. På grund av andra processer som körs i bakgrunden, cachemissar, m.m. varierar körtider av datorprogram alltid något. Återigen finns det flera benchmarkingmått värda att använda i detta sammanhang.

(21)

14 KAPITEL 3. METOD

Det första är att uppskatta den ungefärliga medelvärdestiden algoritmen tar för att lösa problemet. Detta beräknades genom att köra varje algoritm upprepade gånger, summera den totala körtiden och sedan ta fram det aritmetiska medelvärdet [11]. Även tillhörande standardavvikelser beräknades.

Man kan även titta på ytterligheter i körtiderna, det vill säga de uppmätta bästa- respektive värstafallstiderna. I denna rapport kommer bästafallstider- na vara viktigast eftersom de främst representerar de egentliga körtiderna av algoritmerna under optimala förutsättningar.

Samtliga metoder löste den styva Van der Pol-ekvationen (2.8) n = 20 gång- er för vardera toleransgränserna: 10⁻¹, 10⁻², 10⁻³ och 10⁻⁴. Till skillnad från avsnittet ovan har MATLAB inte stöd för tidtagning med 15 decimalers noggrannhet. Istället fick körtider mätas med upp till sex värdesiffrors noggrannhet med hjälp av tidtagningskommandona tic och toc. Eftersom MATLAB använder sig av just in time-kompilering (JIT) och behöver interpretera skript som körs tar filer överlag avsevärt längre tid vid en första körning jämfört med senare gånger [19]. Under ytan utförs dessutom optimeringar av data- strukturer och funktioner som inte är synliga för användaren [11]. För att mot- verka denna effekt och simulera att metoderna kördes för första gången varje gång föregicks varje körning av kommandot clear all.

(22)

Kapitel 4

Resultat

Nedan presenteras resultaten som erhölls från metoderna specificerade i kapitel 3.

4.1 Noggrannhet

Eftersom mittpunktsmetoden är deterministisk så genererar de två implementationerna av den identisk utdata. Därför presenteras noggrannheten av mittpunktsmetoden samlad under en paraplyterm.

Datan i tabellerna 4.1–4.3 nedan komplementeras visuellt av figurer 4.1 och 4.2. Det framgår tydligt att majoriteten av respektive globalfel uppstår i samband med funktionens kurvatur. Detta är i synnerhet iögonfallande för grövre toleranser. Intressant är att mittpunktsmetoden går från att vara den absolut noggrannaste metoden till att bli relativt felande i förhållande till övriga metoder. Notera att skalan på globalfels-axlarna varierar för att kunna belysa skillnaderna mellan kurvorna på en mer detaljerad nivå.

Tabell 4.1: Felskattningar, RADAU5

Tolerans Globalt fel Medelvärdesfel Maxfel

10⁻¹ 3.7 0.3 0.9

10⁻² 0.78 0.05 0.14

10⁻⁵ 0.00050 0.00001 0.00003

10⁻⁷ 0.0000099 0.0000001 0.0000003

15

(23)

16 KAPITEL 4. RESULTAT

Tabell 4.2: Felskattningar, ODE15s

10⁻¹ 2.8 0.1 0.4

10⁻² 1.63 0.05 0.14

10⁻⁵ 0.00609 0.00006 0.00023

10⁻⁷ 0.0004556 0.0000022 0.0000066 Tabell 4.3: Felskattningar, implicit mittpunktsmetod

10⁻¹ 0.6 0.0 0.2

10⁻² 0.26 0.00 0.01

10⁻⁵ 0.01645 0.00001 0.00004

10⁻⁷ 0.0025022 0.0000001 0.0000005

Figur 4.1: Numeriska lösningar och den exakta kurvan för olika toleranser.

(24)

KAPITEL 4. RESULTAT 17

Figur 4.2: Kumulativa globalfel som funktioner av t för olika toleranser.

Samma presentationsordning som användes för den exakta instansen nytt- jas även här. Datan i tabellerna 4.4–4.6 hör ihop med av figurer 4.3–4.4. Mitt- punktsmetoden är utritad med streckade kurvor för att samtliga approximerade kurvor ska vara synbara. Figurerna har inte heller här enhetliga globalfels- axlar, av samma anledning som ovan. Återigen visar funktionens skarpa kur- vaturförändringar vara kopplade till var approximationerna blir som sämst.

Den framträdande storleken på globalfelen för ODE15s vid toleranserna 10⁻¹ och 10⁻² härstammar från den förlängning respektive förkortning av den approximerade lösningskurvans period. Storleken på medelvärdesfelen för dessa två lösningar beror på de relativt små mängder punkter ODE15s har försökt löst ekvationen i.

Tabell 4.4: Felskattningar, RADAU5

10⁻¹ 0.4 0.0 0.1

10⁻² 0.38 0.02 0.07

10⁻³ 0.328 0.016 0.076

10⁻⁴ 0.1767 0.0084 0.0287

(25)

Tabell 4.5: Felskattningar, ODE15s

10⁻¹ 13.6 0.6 3.3

10⁻² 6.81 0.32 3.09

10⁻³ 0.067 0.003 0.015

10⁻⁴ 0.0170 0.0008 0.0050

Tabell 4.6: Felskattningar, mittpunktsmetoden Tolerans Globalt fel Medelvärdesfel Maxfel

10⁻¹ 0.0 0.0 0.0

10⁻² 0.01 0.00 0.01

10⁻³ 0.013 0.001 0.010

10⁻⁴ 0.0125 0.0006 0.0091

Figur 4.3: Numeriska lösningar och simulerad exakt kurva för olika toleranser.

(26)

Figur 4.4: Kumulativa globalfel som funktioner av t för olika toleranser.

4.2 Tidsåtgång

Här skiljs den vektoriserade mittpunktsmetoden från den ovektoriserade implementationen. Samtliga standardavvikelser beräknades med 15 decimalers noggrannhet men har därefter avrundats till 6 signifikanta siffror i alla nedan- stående tabeller (4.7–4.10). De utritade standardavvikelserna i figur 4.5 använ- de sig dock av de oavrundade värdena. Enheterna som används i samtliga tabeller och figurer är sekunder.

(27)

Tabell 4.7: Körtider, RADAU5

Tolerans Snittid Standardavvikelse Värstafall Bästafall 10⁻¹ 0.333807 0.0147395 0.370819 0.312999 10⁻² 0.337528 0.0149815 0.379134 0.320095 10⁻³ 0.382760 0.0334471 0.461156 0.354039 10⁻⁴ 0.448018 0.0552061 0.645812 0.406034

Tabell 4.8: Körtider, ODE15s

Tabell 4.9: Körtider, vektoriserad mittpunktsmetod

Tabell 4.10: Körtider, ovektoriserad mittpunktsmetod

(28)

Figur 4.5: Statistik över genomsnittliga körtider med standardavvikelser.

(29)

Kapitel 5

Diskussion

5.1 Noggrannhet

Som tidigare nämnt framstår det som att globalfelen för samtliga algoritmer växer som mest då ODE:n har som kraftigast svängningar. Detta är som tydli- gast vid grövre toleranser (10⁻¹ och 10⁻²) och mest iögonfallade för RADAU5 i figur 4.1. För toleransen 10⁻¹ vid den exakta instansen ser man tydligt att metoden har tagit en för stor steglängd vid flera lösningspunkter, vilket resulterar till en relativt lågkvalitativ numerisk lösning. Som följd av detta blir även medelvärdesfelet större än för övriga metoder, se tabell 4.1.

Även ODE15s genererar en märkbart avvikande kurva från den exakta lös- ningen vid toleransen 10⁻¹. Metoden visar sig även prestera sämst av de tre metoderna avseende globalfelet vid toleransen 10⁻². Detta är något förvånande utifrån graferna i övre högra hörnet i figur 4.1 (där RADAU5 intuitivt ser ut att ha störst globalfel), men beror på att ODE15s märkbart har felskattat lös- ningen från och med t ≈ 5, se figur 4.2.

Vid finare toleranser än 10⁻² ger RADAU5 i stort lösningar som är oskiljba- ra från den exakta lösningen. Även ODE15s genererar kurvor av noggranna- re karaktär än mittpunktsmetoden som i övriga fall presterat bäst. MATLAB- dokumentation rekommenderar (och har som standard) att samtliga inbyggda ODE-lösare bör köras med toleranser ≤ 10⁻³ [16]. Detta kan vara en förkla- ring till varför felen för ODE15s minskar oproportionerligt för de senare två fallen. Samma resonemang kan appliceras på RADAU5:s goda resultat vid små toleranser, enligt dokumentation från Hairer [14] har metoden en defaultin- ställning för toleransen på 10⁻³. De mer sofistikerade steglängdsväljarna hos RADAU5 och ODE15s lär även vara en bidragande faktor till varför de är nog- grannare än mittpunktsmetoden vid mindre toleranser än så.

22

(30)

KAPITEL 5. DISKUSSION 23

Som det framgår i figurer 4.3-4.4 och avsnittet ovan så skiljer sig metodernas noggrannhet probleminstanserna emellan. Alltså är det inte fullt möjligt att extrapolera om en metods noggrannhet på andra ODE:s utifrån lösningen till en annan. I vilket fall indikerar lösningarna till ekvation (2.8) på att den föreslag- na övre gränsen för toleransen till ODE15s är rimlig. För större toleranser så tas åtminstone ett steg för stort på det ungefärliga intervallet 0.8 ≤ t ≤ 1.7 vilket resulterar i lösningskurvor med deformerade perioder. Att det är just i detta intervall (där punkterna där förstaderivatan till funktionen är som störst i beloppet ligger) som lösningskvaliteten påverkas negativt stöds av styvehets- definitionen enligt Hairer [6]. De trappstegsliknande globalfelskurvorna i figur 4.4 belyser tydligt exakt var dessa punkter ligger på t-intervallet.

Steglängden tas stor för ODE15s vid dessa toleranser och detta påverkar även medelvärdesfelet, vilket är mångdubbelt större än övriga metoders. Faktumet att maxfelet är större än den approximativa amplituden av Van der Pol-kurvan talar ytterligare för att metoden är olämplig att användas vid dessa toleranser.

Någonstans mellan toleranserna 10⁻² och 10⁻³ blir RADAU5 avsevärt mindre noggrann än ODE15s, som i övrigt har haft jämförbar noggrannhet med mittpunktsmetoden. Till skillnad från ovan så växer inte RADAU5-kurvorna i figur 4.4 (för toleranserna 10⁻³ och 10⁻⁴) märkbart vid de kritiska punkterna där derivatan är som störst. Istället ökar de på ett mer linjärt sätt. Visserligen är RADAU5-lösningarna i dessa fall mindre noggranna, men faktumet att de inte verkar påverkas märkbart av ekvationens svängningar kan vara indikativt på att det är en metod som lär fungera bra för andra styva ODE:s.

För samtliga toleranser ger mittpunktsmetoden de absolut minsta globalfelen.

Då metoden tar betydligt fler steg än övriga metoder (upp till två storleksord- ningar fler, beroende på tolerans) blir även medelvärdesfelen mycket små. För- modligen bottnar detta i den simpla steglängdsväljaren metoden använder, jämfört med de som används i RADAU5 och ODE15s. En förfining av denna borde leda till att färre steg krävs, utan att lösningskvaliteten påverkas negativt (eftersom ODE15:s bevisligen presterar nästintill lika bra vid toleranser

≤ 10⁻³).

I ett lite större perspektiv är det även viktigt att värdera hur pass noggranna lösningar man verkligen behöver. Generellt sett varierar detta beroende på användningsområde och situation. Val av ODE-lösare får göras därefter. Kort sagt kan man säga att ODE15s inte bör användas för toleranser större än 10⁻³, medan RADAU5 ser ut att fungera bra för styva ODE:s för just grövre toleranser. Krävs mycket god noggrannhet oavsett tolerans rekommenderas mittpunktsmetoden.

(31)

24 KAPITEL 5. DISKUSSION

5.2 Tidsåtgång och vektorisering

Ur figur 4.5 och tabeller 4.7-4.8 framgår det tydligt att RADAU5 och ODE15s i stort sett är likvärdiga vad gäller körtider, oavsett tolerans. Dock har ODE15s snabbare medelvärdes-, värstafalls- och bästafallstider. Tidsåtgången ökar dessutom sublinjärt i förhållande till toleransen som används vilket tyder på att de kan köras med godtyckligt små toleranser och ändå ha god prestanda, förut- satt att MATLAB kan hantera tillräckligt små tal för Newton-Raphsons metod.

Implementationerna av mittpunktsmetoden verkar dock inte ha denna önsk- värda egenskap. Båda versionerna tar markant längre tid på sig ju mindre toleransen blir. I synnerhet gäller detta för den ovektoriserade implementationen. Som nämndes i diskussionen om noggrannhet lär det finnas möjlighet att förbättra båda versionernas prestanda om ytterligare arbete läggs på steg- ländgsväljaren i metoden.

De uppmätta standardavvikelserna för mittpunktmetoderna är även de värda att titta på. I de allra flesta fall ligger dessa på några tiondelars sekund. Stan- dardavvikelsen för den ovektoriserade mittpunktsmetodens körtid för toleransen 10⁻⁴ är dock nästan tredubbelt större än för övriga toleranser, se tabell 4.10 och figur 4.5.

I vilket fall åskådliggörs effekten av vektorisering i figur 4.5. Körtiderna för den ovektoriserade implementationen ser ut att närma sig den dubbla av den vektoriserade versionen. Detta är inte helt oväntat då den ovektoriserade versionen gör nästan dubbelt så många funktionsanrop som den vektoriserade (ett för varje endimensionell yi-vektor). Därtill behöver allt längre vektorer slås samman i takt med att toleransen minskar och antalet lösningspunkter ökar.

Detta talar för att skillnaden i tidsåtgång implementationerna emellan ökar som en funktion av toleransen.

Återigen behöver man ställa sig frågan hur pass noggranna lösningar man be- höver och sätta detta i proportion till tiden som krävs för att generera dessa.

Utifrån diskussionen i avsnitt 5.1.2 så lämpar sig ODE15s då både snabb- och noggrannhet efterfrågas, RADAU5/ODE15s då enbart snabbhet är av intresse och (den vektoriserade) mittpunktsmetoden om noggrannhet är den viktigaste aspekten att ta hänsyn till.

(32)

Kapitel 6

Slutsats

Sammantaget handlar det hela om en avvägning mellan noggrannhet och tidsåt- gång vid val av ODE-lösare. Detta har tidigare beskrivits av bland andra Pul- liam [8], och belyses ytterligare av den här rapportens resultat. Den implicita mittpunktsmetoden ger i samtliga observerade testfall klart exaktast lösningar, men tar avsevärt längre tid för finare toleranser, förmodligen på grund av en mindre sofistikerad steglängdsväljare. Föga förvånande är den inbyggda och finjusterade metoden ODE15s bäst om både noggrannhet och effektivitet är ef- terfrågat. Lösningarna den genererar har visat sig vara av kraftigt varierande noggrannhet vid grövre toleranser. Detta kan vara indikativt på faktumet att metoden i vanliga fall har en förinställd standardtolerans på 10⁻³. Fram till denna tolerans är nådd är RADAU5 att föredra över ODE15s, då skillnaden i körtid mellan metoderna är nästintill försumbar.

Fördelarna med att använda konceptet vektorisering vid design av MATLAB- funktioner är tydliga. Resultaten visar att den ovektoriserade implementationen tar uppemot två gånger så lång tid vid lösningen av den styva Van der Pol-ekvationen, särskilt då toleransen för Newton-Raphsons metod blir allt mindre. Graden till vilken prestandan kan förbättras av vektorisering uppskattas därför till runt ×2.

Genom att förfina steglängdsväljaren för mittpunktsmetoden går det troligtvis att förbättra prestandan märkbart. Ifall ytterligare funktionsvärden kan viktas och används för att utöka metoden finns det dessutom utrymme för att utnyttja vektorisering till större grad vilket kan leda till förbättrad noggrannhet i för- hållande till sublinjärt många extra funktionsanrop. Förbättrad parallellisering med MATLAB är även möjlig men kommer att kräva särskilda utvecklings- verktyg. Att testa metoderna på andra styva ODE:s är även något som kan undersökas i framtiden, för att kunna föra mindre anekdotiska resonemang.

25

(33)

Litteraturförteckning

[1] B. Bjurel, G. & Lindberg. Methods for Stiff Differential Equations. ACM SIGNUM Newsletter, 8:16–17, 1973.

[2] W. Brannan, J. & Boyce. Differential Equations: An Introduction to Modern Methods and Applications. Wiley, 2015.

[3] G. Hairer, E. & Wanner. Solving Ordinary Differential Equations II: Stiff and Differential-algebraic Problems. Springer, 1996.

[4] G. Bäckström. Praktisk matematik med MATLAB. Studentlitteratur, 1997.

[5] J.-M. Ginoux and C. Letellier. Van der Pol and the History of Relaxation Oscillations: Toward the Emergence of a Concept. Chaos, 22(2):023120, June 2012. doi: 10.1063/1.3670008.

[6] G. Hairer, E. & Wanner. Stiff Differential Equations Solved by Radau Met- hods. Journal of Computational and Applied Mathematics, 111:93–111, 1998.

[7] T. Bui. Explicit and Implicit Methods in Solving Differential Equations.

University of Connecticut Libraries, Honors Scholar Theses, (119), 2010.

[8] T. H. Pulliam. Time Accuracy and the Use of Implicit Methods. NASA Ames Research Center, 1993.

[9] J. O. Curtiss, C. F. & Hirschfelder. Integration of Stiff Equations. Pro- ceedings of the National Academy of Sciences of the United States of America, 38:

235–243, 1952.

[10] Mathworks. MATLAB Data Storage.

http://mathworks.com/help/matlab

/matlab_external/matlab-data.htmlf22019, 2014.

[11] Y. M. Altman. Accelerating MATLAB Performance. Chapman & Hall/CRC, 1st edition, 2014. ISBN 1482211297, 9781482211290.

[12] Mathworks. MATLAB R2013b Release Notes.

https://mathworks.com/help/releases/R2013b/matlab /release-notes.html, 2013.

26

(34)

LITTERATURFÖRTECKNING 27

[13] L. F. Shampine. Vectorized Solution of ODEs in MATLAB with Control of Residual and Error. pages 259–271, 2011.

[14] G. Hairer, E. & Wanner. Fortran Implemen- tation of the Radau IIA Method of Order 5.

http://www.unige.ch/ hairer/prog/stiff/radau5.f, 1996.

[15] C. Engstler. MATLAB Implementation of the Radau IIA Method of Order 5 by Ch. Engstler After the Fortran Code RADAU5 of Hairer/Wanner.

http://na.uni-tuebingen.de/pub/codes/radau5.tar.gz, 1999.

[16] Lawrence F. Shampine and Mark W. Reichelt. The MATLAB ODE Suite. SIAM J. Sci. Comput., 18(1):1–22, January 1997.

ISSN 1064-8275. doi: 10.1137/S1064827594276424. URL http://dx.doi.org/10.1137/S1064827594276424.

[17] J. Flaherty. One-step Methods. Rensselaer Polytechnic Institute Press, 1989.

[18] L. F. Shampine. Error Estimation and Control for ODEs. Journal of Scientific Computing, 25(1):3–16, 2005. ISSN 1573-7691. doi: 10.1007/BF02728979.

URL http://dx.doi.org/10.1007/BF02728979.

[19] Mathworks. MATLAB Execution Engine.

https://mathworks.com/products/matlab /matlab-execution-engine.html, 2015.

(35)

(36)

www.kth.se