Handleiding MATLAB W en BMT trimester 1.1., 1.2.
F.J.L. Martens J.C. van der Meer Faculteit Wiskunde en Informatica 2000/2001
2
Inhoudsopgave 1 Deze handleiding
7
2 Basiselementen van MATLAB
9
1
Het pakket MATLAB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
2
Het starten en stoppen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
3
Opdrachten in MATLAB . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10
4
Het afbreken van berekeningen of van uitvoer . . . . . . . . . . . . . . . . . .
12
5
Onvolledige opdrachten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
6
Variabelen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
7
Het opslaan van databestanden . . . . . . . . . . . . . . . . . . . . . . . . . .
14
8
Het opslaan van de tekst van een MATLABsessie . . . . . . . . . . . . . . . .
14
9
Het beheer van files en het zoekpad . . . . . . . . . . . . . . . . . . . . . . . .
15
10
Wiskundige expressies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
16
11
De arraystructuur
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
18
12
Arrays met ´e´en rij of ´e´en kolom . . . . . . . . . . . . . . . . . . . . . . . . . .
19
13
Array-operaties in MATLAB . . . . . . . . . . . . . . . . . . . . . . . . . . .
20
13.1
Algebra¨ısche array-operaties . . . . . . . . . . . . . . . . . . . . . . . .
20
13.2
Relationele array-operaties . . . . . . . . . . . . . . . . . . . . . . . .
23
13.3
Andere array-operaties . . . . . . . . . . . . . . . . . . . . . . . . . . .
24
13.4
Transponeren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
26
14
Grafieken . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
27
15
Script files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
29
16
Eigen functies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
30
3
4
INHOUDSOPGAVE
17
De helpfaciliteiten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
31
18
Opgaven . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
33
3 De numerieke gereedschapskist
35
1
Inleiding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
35
2
De grafiek . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
35
3
De nulpunten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
36
4
De extrema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
36
5
Numerieke integratie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
38
6
Opgaven . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
39
4 De symbolische gereedschapskist
41
1
Inleiding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
41
2
Expressies met variabelen . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
42
3
Substitutie
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
43
4
Differenti¨eren en integreren . . . . . . . . . . . . . . . . . . . . . . . . . . . .
44
5
Numerieke waarden . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
45
6
Het manipuleren van expressies . . . . . . . . . . . . . . . . . . . . . . . . . .
46
7
Symbolen, strings en getallen . . . . . . . . . . . . . . . . . . . . . . . . . . .
47
8
Opgaven . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
49
5 Lineaire Algebra
53
1
Inleiding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
53
2
Matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
53
2.1
Wat is een matrix? . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
53
2.2
Het invoeren van matrices . . . . . . . . . . . . . . . . . . . . . . . . .
54
2.3
Matrices met symbolische elementen . . . . . . . . . . . . . . . . . . .
55
2.4
Het invoeren van vectoren . . . . . . . . . . . . . . . . . . . . . . . . .
55
2.5
Bijzondere matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . .
57
2.6
Indices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
58
2.7
Matrixbewerkingen in MATLAB . . . . . . . . . . . . . . . . . . . . .
59
INHOUDSOPGAVE
5
2.8
Matrixfuncties in MATLAB . . . . . . . . . . . . . . . . . . . . . . . .
62
Oplossen van stelsels lineaire vergelijkingen . . . . . . . . . . . . . . . . . . .
63
3.1
Het commando A \ b . . . . . . . . . . . . . . . . . . . . . . . . . . . .
63
3.2
De rijgereduceerde trapvorm . . . . . . . . . . . . . . . . . . . . . . .
64
3.3
Het oplossen van stelsels met de Symbolic Toolbox . . . . . . . . . . .
65
Numerieke aspecten bij gebruik van MATLAB . . . . . . . . . . . . . . . . .
66
4.1
Getalrepresentatie en afrondfouten, een voorbeeld . . . . . . . . . . .
66
4.2
Datafouten, een voorbeeld . . . . . . . . . . . . . . . . . . . . . . . . .
67
4.3
Effici¨entie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
68
3
4
5
MATLAB en Lineaire Algebra 5.1
. . . . . . . . . . . . . . . . . . . . . . . . . .
69
Opdrachten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
69
6 Differentiaalvergelijkingen 1
2
3
4
71
M-files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
71
1.1
M-files maken . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
71
1.2
Script- en functionfiles . . . . . . . . . . . . . . . . . . . . . . . . . . .
71
Differentiaalvergelijkingen . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
73
2.1
Het numeriek oplossen van differentiaalvergelijkingen . . . . . . . . . .
73
2.2
Stelsels differentiaalvergelijkingen . . . . . . . . . . . . . . . . . . . . .
75
2.3
Het richtingsveld . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
76
2.4
Plotten van integraalkrommen . . . . . . . . . . . . . . . . . . . . . .
77
Het symbolisch oplossen van differentiaalvergelijkingen . . . . . . . . . . . . .
77
3.1
Eerste-orde differentiaalvergelijkingen . . . . . . . . . . . . . . . . . .
77
3.2
Stelsels eerste-orde differentiaalvergelijkingen . . . . . . . . . . . . . .
78
3.3
Hogere-orde differentiaalvergelijkingen . . . . . . . . . . . . . . . . . .
79
MATLAB en differentiaalvergelijkingen - Opdrachten . . . . . . . . . . . . . .
79
4.1
Opdrachten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
79
Index
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
82
Referenties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
82
6
INHOUDSOPGAVE
Hoofdstuk 1
Deze handleiding Deze handleiding is bedoeld om u snel met enkele aspecten van MATLAB vertrouwd te maken. Het is de bedoeling dat u, verdeeld over twee trimesters in een viertal dagdelen met dit materiaal oefent. Doe dit samen met iemand anders, dat werkt effici¨enter. Als u weinig met computers gewerkt heeft, dan zult u misschien niet alles in deze twee dagdelen af kunnen werken. Haal dit snel op een geschikt moment in! Het programma is als volgt: Trimester 1.1. dagdeel 1 De hoofdstukken “De handleiding”, “Basiselementen van MATLAB” en “De numerieke gereedschapskist”. dagdeel 2 De hoofdstukken “De symbolische gereedschapskist”. Trimester 1.2. dagdeel 1 Het hoofdstuk “Lineaire Algebra”. dagdeel 2 Het hoofdstuk “Differentiaalvergelijkingen”. De trainingen aan het begin van trimester 1.1 sluiten aan bij het college Calculus (2Y130) en de basisvaardigheden uit deze trainingen zullen in dit college, met name bij de oefeningen, verder worden ontwikkeld. Het eerste dagdeel in het tweede trimester sluit aan op het eerste college Lineaire Algebra (2Y650) en is de aanzet tot het gebruik van MATLAB bij lineaire algebra. Het tweede dagdeel in het tweede trimester sluit aan op het zevende college Lineaire Algebra (2Y650) en de colleges Dynamische systemen (4Q220) en Dynamica (4A230) en is een eerste aanzet voor het gebruik van MATLAB bij colleges in het derde trimester. MATLAB speelt een zeer belangrijke rol bij het OGO, het Onwerp Geori¨enteerd Onderwijs. Hierbij zult u MATLAB moeten gebruiken. Eventueel moet u zichzelf met behulp van de 7
8
HOOFDSTUK 1. DEZE HANDLEIDING
MATLABdocumentatie commando’s eigen maken om bepaalde berekeningen uit te kunnen voeren. Daarom is een goede beheersing van de elementaire MATLABcommando’s essentieel. Iedere sectie is voorzien van oefeningen. Voer deze oefeningen nauwgezet uit en schrijf uw uitwerkingen op. Vermeld bij uw uitwerkingen de gebruikte MATLABopdrachten. Als u geen aantekeningen maakt, zijn de oefeningen volstrekt zinloos! Het is verstandig om uw aantekeningen te bewaren zodat u in de toekomst nog eens kunt naslaan hoe iets ook alweer moet. Aan het eind van een trimester moet u zo bedreven zijn in het gebruik van MATLAB, dat u het deel van de handleiding, dat bij de genoten trainingen hoort, niet meer hoeft in te kijken. De hulpfaciliteiten van MATLAB bieden ook alle informatie. Een belangrijk onderdeel van deze trainingen is het leren gebruiken van deze faciliteiten. We gaan ervan uit dat u met Windows95/98 vertrouwd bent. U moet in staat zijn om ‘files’ te beheren. U moet ze kunnen aanmaken, vinden, veranderen en weggooien. Op uw machine moet MATLAB inclusief de ‘Extended Symbolic Math Toolbox’ of de ‘Symbolic Math Toolbox’ en een Internetbrowser staan. De versies van deze programma’s doen minder ter zake. Een commando kan dus net iets anders werken dan in deze handleiding staat, maar met wat experimenteren komt u er vast wel uit. Er zijn enkele redenen om niet al te gedetailleerd te omschrijven hoe u MATLAB op moet starten, welke menukeuzes u moet maken, wat u moet invoeren en hoe de resultaten eruit zien. De eerste reden is dat op uw machine MATLAB op een andere manier ge¨ınstalleerd kan zijn en dat uw versie net even anders is dan de versie die uitgangspunt is geweest bij deze handleiding. De tweede reden is dat bij iedere nieuwe versie op detailpunten de programma’s veranderen. De helpfaciliteiten helpen u wel op de goede weg.
Hoofdstuk 2
Basiselementen van MATLAB 1
Het pakket MATLAB
MATLAB is een pakket voor technische berekeningen. Met het pakket kunnen berekeningen, visualisaties en en zelfgeschreven programma’s gecombineerd worden. MATLAB kan worden uitgebreid met zogenaamde ‘toolboxes’. In het algemeen zijn dit collecties MATLABfuncties die geschikt zijn voor bijzondere klassen van problemen. De “(Extended) Symbolic Math Toolbox” wijkt af. Deze ‘toolbox’ bestaat voor een groot deel uit ingredi¨enten van het pakket Maple dat sterk is in formulemanipulatie. MATLAB voert berekeningen uit met behulp van matrices. Matrices zijn rechthoekige schema’s van getallen en worden ook arrays genoemd. De “Symbolic Math Toolbox” wijkt nogal af van MATLAB zelf. Symbolische berekeningen maken nauwelijks gebruik van arrays. MATLABfuncties maken in tegenstelling van functies uit de “Symbolic Math Toolbox” gebruik van berekeningen met matrices. Om deze reden is het van belang om te weten of een gebruikte functie uit MATLAB of uit de “Symbolic Math Toolbox” is. Op voorhand zal dit niet duidelijk zijn, maar na geregeld gebruik van MATLAB zal dit duidelijk worden.
2
Het starten en stoppen
Starten kan op verschillende manieren: • Klik op ‘Start’ en ‘Run’ en type in de box de naam van de directory, waarin “matlab.exe” is opgeslagen, gevolgd door \matlab . Bij versie 5.3 van MATLAB is dit matlabr11\bin\matlab . • Klik het MATLABicoon aan. In beide gevallen verschijnt het ‘Command Window’ waarin opdrachten kunnen worden ingetypt na de MATLABprompt 9
10
HOOFDSTUK 2. BASISELEMENTEN VAN MATLAB
>> Opdrachten worden afgesloten met de ‘return’-toets. De opdracht wordt dan door MATLAB uitgevoerd. Als MATLAB klaar is, dan verschijnt na eventuele uitvoer een nieuwe prompt. Er kan dan weer een opdracht gegeven worden. Het verlaten van MATLAB kan ook op verschillende manieren: • Door het commando >> quit • Via ‘File\Exit MATLAB’ Oefening 2.1 Start MATLAB, bekijk de mededelingen op het scherm en verlaat het programma.
3
Opdrachten in MATLAB
Na de MATLABprompt >> met de knipperende verticale streep, de cursor |, kan een opdracht worden ingetikt. Na indrukken van de ‘Return’-toets wordt de opdracht uitgevoerd en verschijnt het resultaat eventueel op het scherm. Een opdracht en zijn resultaat ziet er als volgt uit:
>> a = 1 + 2 + 3 a = 6
Het resultaat van 1 + 2 + 3 wordt aan de variabele a toegekend. Met deze variabele kan men verder rekenen. >> b = 2*a + a/3 b = 14
In MATLAB worden variabelen ge¨ıntroduceerd door er een waarde aan toe te kennen. Men kan dan met de variabele verder rekenen en er later eventueel weer een nieuwe waarde aan toekennen. Het is niet mogelijk om berekeningen uit te voeren met een variabele waaraan
3. OPDRACHTEN IN MATLAB
11
geen waarde is toegekend. Een opdracht hoeft niet met een toekenning van de vorm variabele = te beginnen. In zo’n geval wordt het resultaat automatisch aan de variabele ans toegekend. Men kan dan met ans verder rekenen. >> 4*a + 1 ans = 25 >> ans*ans ans = 625
Wilt u nu weten wat de waarde van a is, dan volstaat de volgende opdracht: >> a a = 6
MATLAB zet normaal de uitvoer onder de opdracht. De uitvoer wordt onderdrukt door achter de opdracht een puntkomma te plaatsen. Na uitvoering van de opdracht s = 1 + 2; ziet het scherm er als volgt uit >> s = 1 + 2; >> De variabele s heeft de waarde 3 gekregen. Als u de waarde van deze variabele wil weten, dan kunt u s opnieuw opvragen. >> s s = 3 Als een opdracht meer dan een regel beslaat, dan kunt u voor het eind van de regel drie punten neerzetten en de ‘Return’-toets indrukken. U kunt dan op de volgende regel verder gaan. Het scherm ziet er na het voltooien van de opdracht als volgt uit:
12
HOOFDSTUK 2. BASISELEMENTEN VAN MATLAB
>> s = 1 + 2 + ... 3 + 4 s = 10
Oefening 2.2 Het is de bedoeling dat u onderstaande opdrachten in de aangegeven volgorde uitvoert. Bedenk, alvorens een opdracht uit te voeren, wat er op het scherm verschijnt en wat het resultaat van de opdracht zal zijn.
a.
>> a = 1 + 2 + 3 + 4 + 5 + 6
b.
>> b = 1 + 2 + 3 + 4 + 5 + 6;
c.
>> b + 3
d.
>> 1 + 2 + 3 + 4 + 5 + 6
e.
>> c = 10*ans
4
Het afbreken van berekeningen of van uitvoer
Een langdurige berekening van MATLAB kan met de toetscombinatie ‘Control-C’ worden afgebroken. Na deze combinatie komt er een nieuwe prompt. Nieuwe opdrachten kunnen weer worden uitgevoerd.
Oefening 2.3 Evaluatievan de opdracht pause(30) heeft als effect dat MATLAB 30 seconden pauze houdt. Voer deze opdracht uit en pas snel de toetscombinatie ‘Control-C’ toe. Wat is het effect?
5
Onvolledige opdrachten
Als een opdracht onvolledig is en u drukt op de ‘Return’-toets, dan zijn er twee mogelijkheden: • Er verschijnen foutmeldingen en een nieuwe MATLABprompt. U kunt weer opnieuw aan de slag. • De knipperende cursor staat helemaal links aan het begin van de regel onder de opdracht. Er zijn dan twee mogelijkheden. – De opdracht kan worden afgemaakt en vervolgens worden ge¨evalueerd. – De toetscombinatie ‘Control-C’ zorgt ervoor dat de invoering van de opdracht wordt afgebroken. Er verschijnt ook een nieuwe MATLABprompt.
6. VARIABELEN
13
Oefening 2.4 Door de opdracht a = [2,3,4] wordt aan de variabele a het array (2, 3, 4) toegekend. (a) Voer a = [2,3,4 in en druk op de ‘Return’-toets. Maak de opdracht af. (b) Voer a = [2,3,4 in en druk op de ‘Return’-toets. Onderbreek de invoering van de opdracht.
6
Variabelen
De naam van een variabele moet met een letter beginnen. Daarna mag de naam uit een willekeurig aantal cijfers en/of letters bestaan. MATLAB maakt onderscheid tussen grote en kleine letters. MATLAB gebruikt enkele speciale variabelen: ans eps
i j pi Inf NaN
Deze variabele bevat de uitkomst van de laatste berekening die niet aan een andere variabele is toegekend. De waarde van deze variabele is ongeveer 2.2204 · 10−16 en wordt intern door MATLAB gebruikt. Verwar eps niet met de e-macht. Want e = exp(1) ≈ 2.7183 . Het complexe getal i met de eigenschap i2 = −1. Het complexe getal i. De notatie j voor het getal i komt veel voor. 3.14159 . . . De waarde is oneindig. Delen van 1 door 0, 1/0, levert Inf op. Deze variabele is een representatie van “Not a Number”. Een berekening met NaN resulteert altijd in NaN. De opdracht 0/0 levert NaN op.
Het commando workspace resulteert in een ‘window’ met de naam “Workspace Browser” waarin alle door u zelf gebruikte variabelen vermeld staan. Er staat ook aangegeven tot welke klasse de variabelen behoren. Eventueel kunt u via dit ‘window’ enkele variabelen verwijderen. Oefening 2.5 Gebruik in enkele opdrachten de variabelen u en v correct en verwijder deze vervolgens. Enkele andere belangrijke commando’s voor het variabelenbeheer zijn: who whos clear clear x y
geeft lijst met gebruikte variabelen. geeft lijst met gebruikte variabelen en extra informatie. verwijdert alle variabelen. verwijdert de variabelen x en y.
Waarschuwing: Het is mogelijk om aan de MATLABvariabelen zelf een andere waarde toe te kennen. Berekeningen kunnen hierdoor be¨ınvloed worden. Geef nooit een andere waarde aan een MATLABvariabele. Als u per ongeluk aan een MATLABvariabele een andere waarde geeft, dan kunt u deze waarde verwijderen met behulp van clear of de de “Workspace Browser”.
14
7
HOOFDSTUK 2. BASISELEMENTEN VAN MATLAB
Het opslaan van databestanden
Resultaten van berekeningen van MATLAB worden vaak elders weer opnieuw gebruikt. De waarden van variabelen kunnen bewaard worden in een ‘.mat file’. Dit gebeurt met de opdracht >> save filenaam waardoor alle gebruikte variabelen worden bewaard in de file met naam ‘‘filenaam.mat’’. De extensie ‘.mat’ hoeft niet opgegeven te worden. Hoeven er maar een aantal variabelen te worden opgeslagen, dan dienen die gespecificeerd te worden na de filenaam. De opgeslagen file kan in MATLAB worden ingelezen met de opdracht >> load filenaam waarbij de extensie ‘.mat’ mag worden weggelaten. De variabelen met hun waarde kunnen nu weer gebruikt worden. Oefening 2.6 Voer de twee opdrachten >> f = 2 en >> g = 3 + f uit. Sla de variabelen f en g op in de file “probeer.mat”. Verlaat MATLAB, start MATLAB opnieuw op en laad de file “probeer.mat”. Kijk of de variabelen f en g bekend zijn. Waar is de file “probeer.mat” terecht gekomen? Verwijder de file “probeer.mat”.
8
Het opslaan van de tekst van een MATLABsessie
Het is ook mogelijk om alle in- en uitvoer van een MATLABsessie in een tekstfile op te slaan. U kunt deze tekst niet voor een nieuwe sessie gebruiken, maar u kunt wel zien wat u allemaal gedaan heeft. De opdracht >> diary filenaam zorgt ervoor dat alle in- en uitvoer naar de file “filenaam” wordt weggeschreven. Met de opdracht >> diary off wordt het opslaan be¨eindigd. Na het be¨eindigen kunt u de file “filenaam” bekijken.
9. HET BEHEER VAN FILES EN HET ZOEKPAD
15
Oefening 2.7 Open een tekstfile “afschrift” om in- en uitvoer op te slaan. Voer de opdrachten >> f = 2*3 en >> g = f+2 uit en be¨eindig het opslaan van in- en uitvoer. Bekijk de inhoud van de file “afschrift”. Waar is de file terecht gekomen? Verwijder deze file.
9
Het beheer van files en het zoekpad
Als MATLAB draait, dan is er een directory actief. Aangemaakte files in een MATLABsessie komen hier terecht tenzij er expliciet directories worden vermeld. In MATLAB kunt u opvragen welke directory actief is. >> pwd ans = C:\MATLABR11\work De opdracht is een afkorting van “path working directory”. Het antwoord is natuurlijk afhankelijk van de installatie van MATLAB op uw machine en van de manier van opstarten. Er kan een andere directory, bijvoorbeeld “C:\tmp”, gekozen worden waarin de aangemaakte files terechtkomen. >> cd C:\tmp Bij het zoeken naar files doorloopt MATLAB een zoekpad. De directories waarin zelf gemaakte files staan, moeten in het zoekpad staan, want anders kan MATLAB de zelf gemaakte files niet terugvinden. Maak een directory “D:\matlabtr” op uw D-schijf. De opdracht >> matlabpath laat een lijst van alle directories zien waar MATLAB naar files zoekt. Het toevoegen gaat als volgt: >> addpath D:\matlabtr Het opnieuw opvragen van het zoekpad laat zien of het toevoegen gelukt is. Het opstarten van MATLAB in de directory “D:\matlabtr” en het toevoegen van deze directory aan het zoekpad kan worden geautomatiseerd door een file “startup.m” in die directory
16
HOOFDSTUK 2. BASISELEMENTEN VAN MATLAB
te plaatsten, waar ook al de file “matlabrc.m” staat. Deze file wordt dan bij het opstarten van MATLAB automatisch gelezen en de opdrachten erin worden uitgevoerd. >> which matlabrc C:\MATLABR11\toolbox\local\matlabrc.m In dit geval moet “startup.m” in “C:\MATLABR11\toolbox\local” worden gezet. Maak een tekstfile “startup.m” met de tekst % Voeg mijn directory met eigen M-files toe aan zoekpad. addpath D:\matlabtr % Ga naar mijn eigen directory. cd D:\matlabtr % Zet formaat voor papier op A4. set(0,’DefaultFigurePaperType’,’a4’) % Gebruik centimeters in plaats van inches. set(0,’DefaultFigurePaperUnits’,’centimeters’) en zet deze file op de juiste plaats. Een regel achter een %-teken is commentaar. Het gebruik van deze file zorgt ervoor dat ook meteen enkele papierinstellingen worden veranderd. Verlaat MATLAB en start MATLAB opnieuw op. Controleer uw instellingen met onderstaande opdrachten. >> pwd ans = D:\matlabtr >> get(0,’DefaultFigurePaperType’) ans = A4 Vanaf nu worden tijdens een MATLABsessie zelf gemaakte files zonder plaatsaanduiding automatisch naar “D:\matlabtr” geschreven en ook daar gezocht!
10
Wiskundige expressies
Het invoeren van wiskundige expressies gaat niet letterlijk. De volgende tabellen maken duidelijk welke vertaalslag u moet uitvoeren. De variabelen in deze tabellen stellen getallen voor.
10. WISKUNDIGE EXPRESSIES
17
Bewerkingen standaard MATLAB a+b a + b a−b a - b ab a * b a a / b b ab a ^ b Opmerking De volgorde van de bewerkingen in MATLAB is de gewone volgorde: Eerst machtsverheffen, dan vermenigvuldigen en delen en vervolgens optellen en aftrekken. Constanten Standaard MATLAB e exp(1) π pi i i of j ∞ inf of Inf
Standaard sin(x) cos(x) tan(x) arcsin(x) arccos(x) arctan(x)
Functies MATLAB Standaard √ sin(x) x cos(x) ex tan(x) ln(x) 10 log(x) asin(x) acos(x) |x| atan(x) sign(x)
MATLAB sqrt(x) exp(x) log(x) log10(x) abs(x) sign(x)
Opmerking De ronde haken ‘(’ en ‘)’ worden gebruikt om de volgorde van de berekeningen vast te leggen. Oefening 2.8 Bereken in MATLAB de volgende expressies:
(a)
3 + 22 16 + 52
(b) 42/3 (c) log(e) π (d) sin( ) 4
Oefening 2.9 Geef de variabele x de waarde 2. Voer de expressies aan de linkerkant in MATLAB in en bereken deze. Het resultaat van de berekeningen staat aan de rechterkant.
18
HOOFDSTUK 2. BASISELEMENTEN VAN MATLAB
Expressie x3 6 x √ 1 + x2 1
23 e1+x 3
Resultaat 1.3333 0.8944 1.2599
2
148.4132 2
x sin(x )
-6.0544
arctan(x) 1 + x2
0.2214
11
De arraystructuur
MATLAB gebruikt als basiseenheid voor het rekenen matrices ofwel arrays. Een array is een rechthoekige schema van getallen, die elementen worden genoemd, gerangschikt in m rijen en n kolommen: a11 a12 . . . a1n a21 a22 . . . a2n .. .. .. .. . . . . . am1 am2 . . . amn Een array kan op verschillende manieren worden ingevoerd.
1.
Tussen “[” en “]” waarbij de elementen gescheiden worden door spaties of komma’s en de rijen gescheiden worden door een puntkomma. De opdrachten zijn >> a = [1 2 3;4 5 6;7 8 9] of >> a = [1,2,3;4,5,6;7,8,9]
2.
Tussen “[” en “]” waarbij de elementen gescheiden worden door spaties of komma’s en de rijen steeds op een nieuwe regel met behulp van de ‘Return’-toets worden ingevoerd. De opdracht luidt nu >> a=[1 2 3 4 5 6 7 8 9]
Het resultaat is dan a =
´ EN ´ RIJ OF E ´ EN ´ KOLOM 12. ARRAYS MET E
1 4 7
2 5 8
19
3 6 9
Met arrays kan men als volgt manipuleren:
12
MATLAB
Betekenis van opdracht
a(1,2)
Geef het element uit de eerste rij en de tweede kolom van a
a(1,:)
Geef de eerste rij van a.
a(:,3)
Geef de derde kolom van a.
a(:,2)=[10;10;10]
Verander de tweede kolom van a in een kolom met alleen maar tienen.
a(6)
Geef het 6e element van de rij, ontstaan uit het aan elkaar plakken van de kolommen van a.
Arrays met ´ e´ en rij of ´ e´ en kolom
Een array (matrix) met ´e´en rij wordt ook wel een rij of een rijvector genoemd, en een array (matrix) met ´e´en kolom een kolom of een kolomvector. Als het verschil tussen de opeenvolgende elementen van de rij steeds hetzelfde is, kan de rij ook gevormd worden door ‘variabele’ = ‘beginwaarde’ : ‘stapgrootte’ : ‘eindwaarde’ Hierbij is ‘beginwaarde’ het eerste element van de rij, ‘eindwaarde’ het laatste element van de rij en ‘stapgrootte’ het verschil tussen de opeenvolgende elementen. Bijvoorbeeld >> x = 1:-.2:0 geeft de rij x = 1.0000
0.8000
0.6000
0.4000
0.2000
0
Als het verschil tussen de opeenvolgende elementen van de rij ´e´en is kan de ‘stapgrootte’ worden weggelaten. Bijvoorbeeld >> x=1:5 geeft als resultaat de rij
20
HOOFDSTUK 2. BASISELEMENTEN VAN MATLAB
x = 1
2
3
4
5
Men kan arrays (matrices, vectoren) samenvoegen tot nieuwe arrays. Rangschikt men de arrays in een rij, dan moet het aantal rijen in de arrays overeenkomen. Rangschikt men de arrays in een kolom, dan moet het aantal kolommen in de arrays overeenkomen. Het aaneenplakken van de arrays a = [2 , 3 , 4] en b = [6 , 5] gaat als volgt: >> [a,b] ans = 2
3
4
6
5
Elementen van rijen en kolommen kunnen worden aangegeven door te verwijzen naar de index. In dit geval kan men met ´e´en index volstaan. MATLABopdracht
Resultaat
b(1)
6, het eerste element van b.
a([3,1])
[ 4 2 ], een array bestaande uit het 3e en 1e element van a.
Oefening 2.10 Maak de rijvector met beginwaarde -1, eindwaarde 9 en stapgrootte 0.5 . Oefening 2.11 Maak de rijvector met beginwaarde 9, eindwaarde 0 en stapgrootte −1 .
13
Array-operaties in MATLAB
Een array is een verzameling getallen gerangschikt in rijen en kolommen. Een “arrayoperatie” is een operatie (zoals de optelling of de vermeingvuldiging) die men op overeenkomstige elementen van twee arrays van dezelfde vorm toepast. Wanneer men een bewerking op alle elementen van ´e´en array toepast, dan spreekt men ook wel van een “array-operatie”.
13.1
Algebra¨ısche array-operaties
De bewerkingen optellen en aftrekken worden automatische als arraybewerkingen opgevat. Na de toekenningsopdrachten >> a = [ 1 , 2, 3]; b = [ 6 , 5, 4]; levert het optellen van de arrays
13. ARRAY-OPERATIES IN MATLAB
21
>> a+b ans = 7
7
7
en het aftrekken van de arrays >> a-b ans = -5
-3
-1
De arraybewerkingen vermenigvuldigen (*), delen (/) en machtsverheffen (ˆ) geeft men aan door voor het bewerkingsteken een punt ( . ) te plaatsen. Dus >> a.*b ans = 6
10
12
>> a./b ans = 0.1667
0.4000
0.7500
>> a.^b ans = 1
32
81
Bij al deze array-operaties gelden de volgende conventies: (1) Als bij een array-operatie een van de arrays gelijk is aan een getal, dan wordt dit getal gezien als een array waarvan ieder element gelijk is aan dat getal en waarvan de afmetingen gelijk zijn aan de afmetingen van het andere array. (2) Als ieder element van een array met hetzelfde getal vermenigvuldigd moet worden dan hoeft men geen punt voor het vermenigvuldigingsteken te zetten. Voorbeelden: >> 2.^b
22
HOOFDSTUK 2. BASISELEMENTEN VAN MATLAB
ans = 64
32
16
>> [2 2 2].^b ans = 64
32
16
>> a./2 ans = 0.5000
1.0000
1.5000
3*a ans = 3
6
9
Informatie over deze bewerkingen kan worden opgevraagd met het commando help arith. Het maken van een array van functiewaarden x3 . Beschouw de functie f (x) = 1 + x2 We willen aan de variabele y de rijvector (f (0), f (0.2), . . . , f (1.8), f (2)) toekennen. Eerst maken we een array met de argumenten (0, 0.2, . . . , 1.8, 1). >> x = 0:0.2:2 x = Columns 1 through 7 0
0.2000
0.4000
0.6000
1.8000
2.0000
0.8000
1.0000
Columns 8 through 11 1.4000
1.6000
Het maken van de rij van functiewaarden gaat met array-operaties! >> y = x.^3 ./ (1 + x.^2)
1.2000
13. ARRAY-OPERATIES IN MATLAB
23
y = Columns 1 through 7 0
0.0077
0.0552
0.1588
1.3755
1.6000
0.3122
0.5000
0.7082
Columns 8 through 11 0.9270
1.1506
In feite wordt hier een rij van tellers componentsgewijs gedeeld door een rij van noemers. Oefening 2.12 (a) Voer de opdracht t = x.^3 / (1 + x.^2) uit. U krijgt geen foutmelding omdat u iets heel anders heeft uitgerekend! Opmerking: Ook al krijgt u uitvoer, wees altijd kritisch! (b) Voer de opdracht t = x^3 ./ (1 + x.^2) uit. Wat betekent de foutmelding?
Oefening 2.13 Maak voor de volgende functies f het array (f (0), f (0.1), . . . , f (2)). (a) f (x) = x2 + 2x + 1. (b) f (x) =
3x . 1 + 3x
(c) f (x) =
x √ . 1+ x
13.2
Relationele array-operaties
Relationele array-operaties zijn operaties waarbij arrays elementsgewijs worden vergeleken. Bijvoorbeeld voor de matrices 3 5 20 4 3 11 A= en B = 11 9 1 8 2 4 is de opdracht om te kijken welke elementen van A groter dan die van B zijn: >> A > B ans = 0 1
1 1
1 0
24
HOOFDSTUK 2. BASISELEMENTEN VAN MATLAB
Een 1 geeft aan welke elementen van A groter zijn dan de overeenkomstige van B. Het vergelijken van elementen kan alleen als de matrices A en B dezelfde afmetingen hebben of ´e´en van de matrices een getal is. Zo geeft >> A > 2 ans = 1 1
1 1
1 0
De relationele operaties staan in de volgende tabel. Operatie > >= < <= == ~=
Betekenis groter dan groter dan of gelijk aan kleiner dan kleiner dan of gelijk aan gelijk aan niet gelijk aan
Relationele operaties worden vaak gebruikt in combinatie met de functie find. Bekijk het resultaat van help find. Oefening 2.14 Voer de opdracht a = rand(5,5) uit. Deze opdracht genereert een array met 5 rijen en 5 kolommen waarbij ieder getal willekeurig gekozen wordt tussen 0 en 1. Maak een array met al die elementen van a, die kleiner dan 0.3 zijn. Hint: Gebruik bij deze opdracht find.
13.3
Andere array-operaties
In MATLAB werken de wiskundige standaardfuncties sin, cos, sqrt, atan, asin etc. automatisch op array’s. >> a = [1 2 3]; >> sin(a) ans = 0.8415
0.9093
0.1411
1.4142
1.7321
>> sqrt(a) ans = 1.0000
13. ARRAY-OPERATIES IN MATLAB
25
Beschouw de functie f (x) = x sin(x). Een rij van functiewaarden wordt met de volgende opdracht gegenereerd >> x = 0 : 0.25 : 4*pi; >> y = x .* sin(x) Het resultaat is een array met 51 functiewaarden. t √ . 1+ t Maak m.b.v. array-operaties en de functie sqrt de rij (f (0), f (0.1), f (0.2), . . . , f (1)).
Oefening 2.15 Beschouw de functie f (t) =
Enkele andere operaties zijn >> size(A) Het resultaat geeft de afmetingen van A aan, dat wil zeggen het aantal rijen en het aantal kolommen. >> sum(A) Deze functie berekent de kolomsommen van een rechthoekig array A en geeft deze sommen weer in een rijvector. (Een kolomsom is de som van de elementen van een kolom.) Als A een rij is, dan wordt de som van de elementen van A berekent. De functie >> prod(A) doet hetzelfde voor het product.
Oefening 2.16 Geef een ´e´enregelige opdracht die bij een array a het aantal elementen van a bepaalt die groter dan 5 zijn. Test uw opdracht uit met onderstaande arrays. (a) a = [1, 2, 3; 4, 5, 6; 7, 8, 9] (b) a = 1:10 (b) a = 10*rand(6,6) Dezelfde opdracht moet in alle drie de gevallen werken!
Oefening 2.17 Geef een ´e´enregelige opdracht die bij een array a de gemiddelde waarde van de elementen van a berekent. Test uw opdracht uit met onderstaande arrays.
26
HOOFDSTUK 2. BASISELEMENTEN VAN MATLAB
(a) a = [1, 2, 3; 4, 5, 6; 7, 8, 9] (b) a = 1:10 (b) a = rand(6,6) Dezelfde opdracht moet in alle drie de gevallen werken!
13.4
Transponeren
In een array (matrix) kunnen de rijen en kolommen verwisseld worden. Dit heet transponeren. Het resultaat wordt de getransponeerde array genoemd. Het transponeren gebeurt in MATLAB met de functie transpose. Na de opdrachten >> a = [1 2 3] a = 1
2
3
>> A = [1 2;3 4;5 6] A = 1 3 5
2 4 6
gaat het transponeren als volgt: >> transpose(a) ans = 1 2 3 en >> transpose(A) ans = 1 2
3 4
5 6
14. GRAFIEKEN
27
Het transponeren kan ook met behulp van de combinatie van een punt en een accent ( .’ ). De opdrachten transpose(A) en A.’ leveren hetzelfde resultaat op.
14
Grafieken
Grafieken worden getekend met de ‘plot’-opdracht. Het commando >> plot(v) waarbij v een rijvector is die de getallen v1 t/m vN bevat, genereert een plaatje waarin de punten (1, v1 ), (2, v2 ) t/m (N, vN ) achtereenvolgens met lijnstukken verbonden zijn. Het commando >> plot(v,w) met v en w twee rijvectoren met evenveel elementen, tekent de punten (vi , wi ) en verbindt deze achtereenvolgens met lijnstukken. Het commando >> plot(v,w,x,y) tekent in dezelfde figuur de punten (vi , wi ) en (xi , yi ). sin(x) + 2x willen tekenen. Het tekenen 1 + x2 van de grafiek van f (x) met MATLAB gebeurt door een groot aantal punten (xi , f (xi )) achtereenvolgens met rechte lijnstukken te verbinden. De vloeiendheid van de grafiek wordt natuurlijk door het aantal punten bepaald. Het plotten van de grafiek kan met de volgende opdrachten:
Veronderstel dat we de grafiek van de functie f (x) =
>> x=0:pi/100:2*pi; >> y=(sin(x)+2*x)./(1+x.^2); >> plot(x,y) Bij het tekenen van grafieken worden de assen automatisch gekozen. Men kan zelf de lengte van de assen kiezen met >> axis([XMIN XMAX YMIN YMAX]) waarbij de grafiek getekend wordt in de rechthoek XMIN ≤ x ≤ XMAX en YMIN ≤ y ≤ YMAX. Na iedere nieuwe ‘plot’-opdracht verdwijnt de informatie van de oude ‘plot’. Wil men verschillende grafieken met verschillende ‘plot’-opdrachten in dezelfde figuur laten tekenen, dan moet men na de eerste ‘plot’-opdracht het commando
28
HOOFDSTUK 2. BASISELEMENTEN VAN MATLAB
>> hold on geven. MATLAB onthoudt dan de grafieken. Het commando >> hold off herstelt de toestand waarin iedere volgende ‘plot’-opdracht de voorgaande plot doet verdwijnen. Voor het omgaan met het grafische scherm zijn ook de volgende commando’s van belang: Het tonen van een grafiek gaat door activeren of met de opdracht >> shg Het grafische scherm blijft bewaard tot een nieuw ‘plot’-commando wordt gegeven of tot het wordt leeggemaakt. De opdracht voor het leegmaken is >> clf Voor het plaatsen van tekst en dergelijke in of bij een tekening zijn de volgende commando’s van belang: >> grid voorziet de grafiek van een raster. >> title(’titel’) zet een titel boven een ‘plot’, en wel de string aangegeven tussen de aanhalingstekens, in bovenstaand voorbeeld het woord “titel”. >> text(x,y,’string’) zet in het punt (x, y) van de laatste ‘plot’ de tekst “string”. >> xlabel(’tekst’) schrijft informatie langs de x-as. >> ylabel(’tekst’) doet hetzelfde voor de y-as. Het opslaan van een figuur in een file gebeurt met het commando
15. SCRIPT FILES
29
>> print -depsc filename MATLAB maakt een file aan volgens het opgegeven ‘format’, hier aangegeven met de optie -depsc ( Level 1 color Encapsulated PostScript (EPS) ). De file wordt weggeschreven naar de file “filename.eps”. Figuren in EPS-formaat kunnen in documenten van allerlei soort worden toegevoegd. Oefening 2.18 Teken de functies sin(t) en cos(t) op het interval 0 ≤ t ≤ 2π in een figuur.
Oefening 2.19 Teken de functie namen van de assen.
√
t sin(2t) op het interval [0, π]. Zet bij deze grafiek de 2
Oefening 2.20 Teken de functie e−t op het interval [−2, 3].
15
Script files
Als het maken van een plaatje uit meerdere opdrachten bestaat, dan kan men beter de opdrachten bij elkaar in een file zetten en deze door MATLAB in een keer laten uitvoeren. De file moet een naam met achtervoegsel ‘.m’ hebben. Zo’n file waarin opdrachten staan die moeten worden uitgevoerd heet een ‘script file’. Het tekenen van de functie f (x) = |x sin(2πx)| op het interval [0, 2] met een ‘script file’: Via het menu “File\New\M-file” wordt de “MATLAB Editor/Debugger” geopend. Tik de volgende regels in. x = 0:0.1:2; y = abs(x .* sin(2*pi*x)); plot(x,y) title(’f(x) = |x sin(2 pi x)| op interval [0,2]’) axis([0 2 0 2]) shg Sla de file op via het menu “File\Save” onder de naam “plaat.m”. Als u in MATLAB de opdracht plaat geeft, dan zal de grafiek van f op het interval [0, 2] verschijnen. Het is beter om de opdrachten, die een figuur genereren, te bewaren dan de figuur zelf. Oefening 2.21 Waarom kan men beter de opdrachten voor een figuur dan de figuur zelf bewaren?
Oefening 2.22 Beschrijf de betekenis van de regels in de file “plaat.m”.
30
16
HOOFDSTUK 2. BASISELEMENTEN VAN MATLAB
Eigen functies
In MATLAB kan men ook eigen functies defini¨eren. MATLAB gaat er standaard van uit, dat alle functies op arrays werken. Men moet bij het maken van een eigen functie rekening met array-operaties houden. Men kan dan de eigen functies met MATLABfuncties combineren. Een functiedefinitie moet worden opgeslagen in een ‘fuction file’, een file met de extensie ‘.m’. De naam van de file moet gelijk zijn aan de naam van de functie met de extensie ‘.m’. Beschouw de functie f (x) = x2 + ex . In MATLAB maken we een functie met dezelfde naam. Deze definitie moet worden opgeslagen in de file “f.m”. Via het menu ‘File\New\M-file” wordt de “MATLAB Editor/Debugger” geopend. Als men de volgende twee regels intikt function y = f(x) y = x.^2 + exp(x); en de file onder de naam “f.m” opslaat, dan heeft men binnen MATLAB de beschikking over de functie f . We maken enkele opmerkingen over bovenstaande file. • De eerste regel van de file moet het woord “function” bevatten. De gebruikte variabelen zijn locaal en MATLAB zelf hoeft ze niet te kennen. • Als x een array is, dan wordt y een array van functiewaarden. • De puntkomma op het eind verhindert dat bij iedere functie-evaluatie er onnodige uitvoer op het scherm komt. Test altijd de functie f uit. Oefening 2.23 Voer de definitie van f zelf uit. Waar komt de file “f.m” te staan?
U kunt de definitie van f veranderen m.b.v. het commando edit f. Oefening 2.24 Laat in de definitie van de functie f de puntkomma weg. Bereken enkele waarden. Wat is het effect?
x , x≥0 0 , x<0 Maak een eigen functie g in MATLAB. Gebruik hierbij de MATLABfunctie max. Controleer uw eigen functie met Oefening 2.25 Beschouw de functie g(x) =
>> x = -1:0.2:1; >> g(x)
17. DE HELPFACILITEITEN
31
ans = Columns 1 through 7 0
0
0
0
0.8000
1.0000
0
0
0.2000
Columns 8 through 11 0.4000
0.6000
Waarschuwing: De namen van functies en variabelen mogen niet hetzelfde zijn. Oefening 2.26 Maak een eigen functie met de naam h bij de functie h(x) = x2 die bestand is tegen arrray-operaties. Voer de volgende opdrachten in de juiste volgorde uit. (a)
>> h([1,2,3,4])
(b)
>> h = 2.5
(c)
>> h([1,2,3,4])
(d)
>> h(1)
(e)
>> h(1.4)
(f)
>> clear h
Wordt bij de onderdelen (c) t/m (e) h als een functie of als een variabele gezien? Waarom volgt er bij onderdeel (d) geen foutmelding? Kent MATLAB de functie h nog?
17
De helpfaciliteiten
Dit is een van de belangrijkste paragrafen. In principe is alle informatie over MATLAB ook te vinden in de helpfaciliteiten van MATLAB. Deze handleiding is tot nu toe zo uitgebreid geweest om u een snelle start in MATLAB te geven. Omdat de commando’s en functies in MATLAB vaak in meer algemene situaties gebruikt kunnen worden, is de informatie algemener dan voor u nodig is. U zult de voor u essenti¨ele bestanddelen uit de aangeboden informatie moeten halen. We behandelen drie manieren om informatie te verkrijgen.
Directe informatie Met het commando
32
HOOFDSTUK 2. BASISELEMENTEN VAN MATLAB
>> help of >> help ‘‘trefwoord’’ Zonder specificatie krijgt u een algemene lijst van commando’s, functies en ‘toolboxes’ die binnen MATLAB ter beschikking staan. Met specificatie krijgt u gerichte informatie en worden ook vaak suggesties gedaan om bij andere trefwoorden informatie op te vragen. Oefening 2.27 Vraag informatie op over de sinus en de functie load.
Informatie via het “Help Window” Via aanklikken van ‘Help\Help Window’ of door de opdracht >> helpwin Dan verschijnt het gevraagde ‘window’ en kan men via klikken informatie opvragen. Oefening 2.28 Open het “Help Window” en dubbelklik op de regel met “matlab\elfun”. Er verschijnt een lijst met alle elementaire wiskundige functies in MATLAB. Dubbelklik op de regel met “atan” en bekijk de informatie. Probeer ook de mogelijkheden van het “See also”-menu uit. De mogelijkheid “More sin help(HTML)” hangt samen met de derde manier van informatie opvragen. Informatie via de ‘hypertext’-documentatie Via aanklikken van ‘Help\Help Desk(HTML)’ of door de opdracht >> helpdesk Dan wordt via een ‘Web browser’ de “Help Desk” opgestart. Via klikken en menu’s kan men nu gericht naar informatie zoeken. Het zoeken van de informatie kan alfabetisch of per onderwerp. Oefening 2.29 Zoek via de “Help Desk” informatie over de sinus op.
Oefening 2.30 Klik in “Help Desk” achtereenvolgens ‘Getting Started’ en ‘Command Line Editing’ aan. Probeer zelf de pijltjestoetsen en de ‘Escape’-toets uit. Via de “Help Desk” is het ook mogelijk de offici¨ele MATLABdocumentatie te bekijken. Oefening 2.31 Probeer de handleiding “Getting Started with MATLAB” te bekijken. Klik hiertoe ‘Getting Started’ aan in de “Help Desk” en vervolgens het PDF-icoon van de genoemde handleiding aan. Als alles goed gaat, verschijnt er een ‘Window’ met rechts een titelpagina en links een strook met driehoekjes en rechthoeken. Klik er een paar aan en zie wat er gebeurt. Als het bovenstaande niet lukt, sla deze oefening dan verder over. In de loop van de tijd zult u de helpfaciliteiten steeds gemakkelijker gaan gebruiken. De “Help Desk” geeft het gemakkelijkst toegang tot allerlei informatie.
18. OPGAVEN
18
33
Opgaven
Oefening 2.32 Maak in MATLAB een rijvector met de vijfde machten van de eerste 40 natuurlijke getallen. Oefening 2.33 Bereken in MATLAB het exacte product 1 · 2 · . . . · 24 · 25.
Oefening 2.34 Bereken in MATLAB de exacte som 1 + 23 + 33 + 43 + . . . + 1993 + 2003 . Oefening 2.35 Teken de grafiek van de functie f (x) = x/(1 + x) op het interval [0, 4].
Oefening 2.36 Teken de cirkel met straal 1 en middelpunt (0, 0) in het vierkant met −2 ≤ x ≤ 2 en −2 ≤ y ≤ 2. Zorg ervoor dat de cirkel er rond uitziet! Hints: Maak eerst een lijst van x-co¨ordinaten en een lijst van y-co¨ordinaten. Teken de cirkel. Pas de vorm van de figuur aan door het commando axis in verschillende combinaties te gebruiken. Oefening 2.37 De periodieke functie f met periode 2 is gegeven door 2 , 0≤t≤2, t f (t − 2) , 2 ≤ t , f (t) = f (t + 2) , t < 0 . Definieer deze functie in MATLAB. Maak gebruik van mod. Voor een getal r geeft mod(r,2) een getal s terug, zodanig dat 0 ≤ s < 2 ´en s en r een veelvoud van 2 van elkaar verschillen. Dus f (r) = f (s) = s2 . Zorg ervoor dat de eigen functie ook op arrays werkt.
34
HOOFDSTUK 2. BASISELEMENTEN VAN MATLAB
Hoofdstuk 3
De numerieke gereedschapskist 1
Inleiding
Het pakket MATLAB z`elf is de numerieke ‘toolbox’. We maken enkele algemene opmerkingen en behandelen enkele belangrijke commando’s. Bekijk van alle genoemde commando’s de helpinformatie. Het aantal cijfers van getallen op het scherm kan veranderd worden met behulp van de commando’s format long of format short. De berekeningen zelf veranderen niet. MATLAB rekent intern met 16 cijfers nauwkeurig. We gaan de functie f (x) = x3 − x2 − 3 arctan(x) + 1 op het interval [−2, 3] onderzoeken. Maak bij een te onderzoeken functie altijd een ‘M-file’. Oefening 3.1 Maak een eigen functie f met f (x) = x3 − x2 − 3 arctan(x) + 1.
2
De grafiek
Het tekenen van de grafiek van f kan op de eerder beschreven manier met behulp van een array van functiewaarden. Een andere manier is gebaseerd op het gebruik van ‘strings’.
>> ezplot(’f(x)’,[-2,3])
De expressie f(x) moet tussen enkele aanhalingstekens (’) staan. Het interval wordt als een rijvector bestaande uit het linkereindpunt en rechtereindpunt gegeven. Het voordeel van deze opdracht is dat het interval gemakkelijk is te veranderen om de grafiek van f gedetailleerder te onderzoeken. 35
36
3
HOOFDSTUK 3. DE NUMERIEKE GEREEDSCHAPSKIST
De nulpunten
Nulpunten van verreweg de meeste functies kunnen niet exact bepaald worden. Deze kunnen dan wel numeriek benaderd worden. Men spreekt in het algemeen toch van het numeriek bepalen van nulpunten. De functie f heeft drie nulpunten. Het eerste nulpunt ligt in het interval [−2, −1] en de functiewaarden f (−2) en f (−1) zijn verschillend van teken. >> fzero(’f’,[-2,-1]) ans = -1.2780 >> fzero(’f(x)’,[-2,-1]) ans = -1.2780 Het is dus mogelijk om ’f’ of ’f(x)’ te gebruiken. In het laatste geval moet per se de variabele x gebruikt worden. >> f(ans) ans = -4.4409e-016 Aan de laatste uitvoer kunt u zien dat MATLAB een benadering voor het nulpunt geeft. Oefening 3.2 Laat zien dat de benadering hooguit 10−4 van het nulpunt afligt.
Oefening 3.3 Bepaal de overige nulpunten van de functie f numeriek.
4
De extrema
In het algemeen kunnen ook de extrema niet exact bepaald worden. In principe kunnen (de plaatsen van) de extrema numeriek worden benaderd. Men spreekt ook wel van het numeriek bepalen van de extrema. De plaatsen waar de extrema van f ongeveer liggen, kunnen uit de grafiek worden afgelezen. Er is een minimum in het interval [0, 2]. >> x1 = fmin(’f’,0,2)
4. DE EXTREMA
37
x1 = 1.0878 >> x1 = fmin(’f(x)’,0,2) x1 = 1.0878 >> f(x1) ans = -1.3784 De functie f heeft een minimum in x = 1.0878 met waarde −1.3784 . Het is wederom mogelijk om ’f’ of ’f(x)’ te gebruiken. In het laatste geval moet per se de variabele x gebruikt worden. Het interval, waarin het minimum ligt, moet als een rijvector worden ingevoerd. De functie f heeft een maximum op het interval [−1, 0]. Helaas heeft MATLAB geen commando fmax. Oefening 3.4 Teken de grafiek van de functie x 7→ −f (x) op het interval [−2, 3]. Verander de file “f.m” niet! De volgende opdrachten zullen nu duidelijk zijn. >> x2 = fmin(’-f(x)’,-1,0) x2 = -0.5902 >> f(x2) ans = 2.0456 De functie f heeft in x = −0.5902 een maximum ter waarde 2.0456 . Door fmin wordt de plaats van ´e´en minimum in het opgegeven interval benaderd. De waarde van het minimum moet u zelf uitrekenen. De nauwkeurigheden waarmee de plaatsen van de extrema worden uitgerekend, kunnen worden aangepast.
38
HOOFDSTUK 3. DE NUMERIEKE GEREEDSCHAPSKIST
Oefening 3.5 Zorg ervoor dat de getallen met 15 cijfers op het scherm verschijnen. Voer de opdracht >> fmin(’f(x)’,0,2,[0,10^n]) vier keer uit, waarbij u achtereenvolgens de exponent n vervangt door −4, −6, −8 en −10. Bekijk de gelijkblijvende cijfers van de verschillende uitkomsten. Oefening 3.6 De functie f (x) = (x3 − 2)2 heeft precies ´e´en nulpunt. Bepaal het nulpunt numeriek en vergelijk de gevonden waarde met de werkelijke waarde. Hint: Teken de grafiek van f over het interval [0, 2]. Oefening 3.7 Van de functie f (x) = e2x −4ex +2 moeten alle nulpunten en extrema bepaald worden. Kies een geschikt interval om de functie met MATLAB te tekenen. Het interval is geschikt als alle nulpunten en extrema erin liggen en duidelijk te zien is waar deze liggen. Leg uit waarom het een geschikt interval is. Bepaal de nulpunten en extrema (aard en waarde).
5
Numerieke integratie Z
b
Integralen van de vorm f (x)dx kunnen meestal niet exact worden bepaald. a Z 3 De integraal (x3 − x2 − 3 arctan(x) + 1)dx is wel exact te berekenen. De uitkomst is 115 12
−2
− 9 arctan(3) + 23 log(2) + 6 arctan(2) ≈ 6.0245344593535 .
Er zijn verschillende commando’s om de integraal numeriek te bepalen (d.w.z. te benaderen). Zie ook de helpfaciliteiten. We gebruiken het commando quad en laten de getallen in 15 cijfers op het scherm verschijnen. >> quad(’f’,-2,3) ans = 6.02456548407172 Helaas is het niet mogelijk om ’f(x)’ als eerste argument te gebruiken. De benadering wijkt al in het vijfde cijfer af. ‘In principe’ probeert quad vier correcte cijfers te geven. De nauwkeurigheid van de benadering kan worden opgeschroefd. >> quad(’f’,-2,3,[10^-7 10^-7]) ans = 6.02453446058438
6. OPGAVEN
39
Z
3
1 )dx numeriek. x2 1 Experimenteer met nauwkeurigheden. Vergelijk de resultaten met 88 3 , de exacte waarde van de integraal.
Oefening 3.8 Bepaal met behulp van quad de integraal
(x3 + x2 +
Zorg ervoor dat de getallen weer met vijf `a zes cijfers op het scherm verschijnen.
6
Opgaven
Oefening 3.9 Beschouw de functie f (x) = x2 exp(−x). Beantwoord de onderdelen (a) t/m (c) zonder gebruik te maken van de afgeleide f 0 . (a) Waarom heeft de functie f een minimum in x = 0? (b) Waarom heeft f een maximum op het interval (0, ∞)?. (c) Bepaal met fmin een maximum op (0, ∞). (d) Hoever ligt de in onderdeel (c) gevonden benadering af van de werkelijke waarde?
Oefening 3.10 Beschouw de functie f (x) = x sin2 (x).
(a) Teken de functie f op het interval [0, 2π]. (b) Bepaal numeriek de plaats en de waarde van de maxima van f op het interval [0, 2π]. Z 2π (c) Bepaal numeriek de integraal f (x) dx. 0
40
HOOFDSTUK 3. DE NUMERIEKE GEREEDSCHAPSKIST
Hoofdstuk 4
De symbolische gereedschapskist 1
Inleiding
Dank zij de “(Extended) Symbolic Math Toolbox” is het mogelijk om binnen MATLAB aan formulemanipulatie te doen. In toekomstige versies van MATLAB zal deze ‘toolbox’ verder in MATLAB worden ingebouwd en er zullen nog de nodige veranderingen optreden. Voer de volgende opdracht uit: >> syms x Als er geen foutmeldingen verschijnen, dan heeft u de beschikking over de “(Extended) Symbolic Math Toolbox”. De betekenis van deze opdracht wordt later uitgelegd. Met deze ‘toolbox’ wordt het aantal functies uitgebreid. Als u in de “Help Desk” aan de rechterkant ‘Symbolic Math Toolbox Ref’ aanklikt, dan verschijnt er een lijst van functies die ´of nieuw zijn ´of aangepast zijn vanwege deze ‘toolbox’. Voor directe informatie is het help commando geschikt. >> help atan ATAN Inverse tangent. ATAN(X) is the arctangent of the elements of X. See also ATAN2. Overloaded methods help sym/atan.m >> help sym/atan ATAN
Symbolic inverse tangent. 41
42
HOOFDSTUK 4. DE SYMBOLISCHE GEREEDSCHAPSKIST
Oefening 4.1 Verwijder alle variabelen uit de “Work Space”.
2
Expressies met variabelen
In MATLAB kan men met variabelen, die een waarde hebben, werken. Als men in een expressie een variabele gebruikt die geen waarde heeft, dan krijgt men een foutmelding. >> clear x >> exp(-x)+sin(x) ??? Undefined function or variable ’x’. In MATLAB kunnen ‘strings’ worden ingevoerd, zoals ’exp(-x) + sin(x)’ of ’x^2 + 3 x’. Het zijn expressies, aan beide kanten voorzien van het ’-teken. In feite hebben we dit gedaan bij ezplot, fzero en fmin. Beschouw het resultaat van >> ezplot(’exp(-x) + sin(x)’,[0,2*pi]) De functie f heeft in het interval [3, 4] een nulpunt en in het interval [4, 5] een minimum. >> fzero(’exp(-x)+sin(x)’,[3,4]) Zero found in the interval: [3, 4]. ans = 3.1831 >> fmin(’exp(-x)+sin(x)’,4,5) ans = 4.7213 Het is zonder meer mogelijk om aan een variabele een string toe te kennen. >> clear x >> y = ’exp(-x) + sin(x)’ y = exp(-x) + sin(x) De variabele y is nu ge¨ıntroduceerd, omdat deze een ‘string’ als waarde gekregen heeft, maar de variabele x is in MATLAB onbekend. Het gebruik van x aan de rechterkant van het is-teken geeft in dit geval g´e´en foutmeldingen. Merk op dat aan de uitvoer niet te zien is of
3. SUBSTITUTIE
43
het een expressie met symbolen of een ‘string’ betreft! Men kan MATLAB met symbolen laten werken als men de beschikking heeft over de “(Extended) Symbolic Math Toolbox” en variabelen tot symbool verklaart. Dit gaat als volgt. >> x=sym(’x’) x = x of >> syms x Het declareren van meerdere variabelen of symbolen tegelijkertijd: >> syms x y z Men kan nu met deze variabelen naar hartelust manipuleren. >> y = x^3 + x^2 + 1 y = x^3+x^2+1 >> v=sin(x) v = sin(x) >> y*v ans = (x^3+x^2+1)*sin(x)
3
Substitutie
>> syms a b >> y = 3*sin(a) + cos(b)
44
HOOFDSTUK 4. DE SYMBOLISCHE GEREEDSCHAPSKIST
y = 3*sin(a)+cos(b) >> y =subs(y,a,2) y = 3*sin(2)+cos(b) >> subs(y,b,5) ans = 3.0116
4
Differenti¨ eren en integreren
We illustreren de opdrachten aan de hand van het voorbeeld: >> syms x y >> y = atan(x) Een keer en drie keer differenti¨eren: y = atan(x) >> diff(y,x) ans = 1/(1+x^2) >> diff(y,x,3) ans = 8/(1+x^2)^3*x^2-2/(1+x^2)^2 Z De onbepaalde integraal >> int(y,x)
arctan(x) dx:
5. NUMERIEKE WAARDEN
45
ans = x*atan(x)-1/2*log(x^2+1) Z De bepaalde integraal
7
arctan(x) dx: 2
>> int(y,x,2,7) ans = 7*atan(7)-1/2*log(2)-1/2*log(5)-2*atan(2) Dit is de exacte waarde van de integraal. De numerieke waarde is dan >> double(ans) ans = 6.6367 De opdracht double maakt van een exact getal een decimaal getal. Verreweg de meeste integralen kunnen niet worden uitgerekend. Beschouw het voorbeeld Z 4 p arctan(x) 1 + x3 dx. 1
>> int(exp(-x) * sqrt(1 + x^3), x, 1, 4) ans = int(exp(-x)*(1+x^3)^(1/2),x = 1 .. 4) >> double(ans) ans = 1.0017 De MATLABfunctie double zorgt voor een numerieke benadering van de integraal.
5
Numerieke waarden
De functie sym zet numerieke getallen om naar exacte getallen en double doet het omgekeerde.
46
HOOFDSTUK 4. DE SYMBOLISCHE GEREEDSCHAPSKIST
MATLAB sym(sqrt(2)) sym(2.3654) double(sin(2)+log(3))
Resultaat sqrt(2) 11827/5000 2.0079
Z
3
(cos(x − 1) +
Oefening 4.2 Bereken in MATLAB de integraal 1
Wat is de numerieke waarde van deze integraal?
1 ) dx. x
De manier om een idee te krijgen van de grootte van een exacte uitdrukking is het toepassen van het commando double.
6
Het manipuleren van expressies
Voer eerst de opdracht >> syms a b c d x y uit en experimenteer met de onderstaande opdrachten. MATLAB expand((a+b)^3) factor(a^3+3*a^2*b+3*a*b^2+b^3) [t n] = numden(a/b+c/d) simplify((x^2+2*x+1)/(x+1))
Betekenis van opdracht Werk de haakjes weg. Ontbind in factoren. Breng onder een noemer. Dan wordt t de teller en n de noemer. Vereenvoudig de uitdrukking.
Het schrijven van de uitdrukking x/(1 + x/(1 + x)) als teller gedeeld door noemer gaat als volgt: >> syms x >> y = x/(1+x/(1+x)); >> [t,n]=numden(y) t = x*(1+x)
n = 1+2*x Teller en noemer zijn niet verder te vereenvoudigen.
7. SYMBOLEN, STRINGS EN GETALLEN
7
47
Symbolen, strings en getallen
De functie f (x) = 2xe−x − sin(x) wordt met MATLAB onderzocht op het interval [0, π]. Voer de volgende commando’s uit. >> syms x >> y = 2*x*exp(-x)-sin(x) Nu zijn de variabelen x en y in MATLAB als symbolen gedefinieerd. Het tekenen van de functie f (x) in MATLAB gaat met het commando >> ezplot(y,[0,pi]) Aan het plaatje is meteen te zien dat f een maximum, een minimum en twee nulpunten heeft in het interval [0, π]. De plaatsen van de extrema en de nulpunten zijn niet exact te bepalen. Deze punten moeten numeriek benaderd worden. De MATLABfuncties fmin en fzero werken echter met strings, dat wil zeggen uitdrukkingen die tussen enkele aanhalingstekens (’) staan. De strings mogen expressies in de variabele x zijn. Van een symbolische uitdrukking kan eenvoudig een string gemaakt worden met de functie char. >> fzero(char(y),[0.5,1]) Zero found in the interval: [0.5, 1]. ans = 0.8030 >> fmin(char(y),1.5,2) ans = 1.8409 >> fzero(char(y),[2.5,3]) Zero found in the interval: [2.5, 3]. ans = 2.7923 De plaats van het minimum en de nulpunten zijn gevonden. Ook de plaats van het maximum is gemakkelijk te bepalen. >> fmin(char(-y),0,0.5) ans =
48
HOOFDSTUK 4. DE SYMBOLISCHE GEREEDSCHAPSKIST
0.3384 Oefening 4.3 Bepaal de extreme waarden van de functie f in het interval (0, π). De randpunten van het interval doen niet mee. Het zal vaak voor komen dat u van een symbolische uitdrukking een string moet maken of omgekeerd. U kunt met class altijd achterhalen met wat voor expressie u bezig bent. >> class(y) ans = sym >> z = char(y); >> class(z) ans = char Veronderstel dat u een array van functiewaarden van f wil maken zonder eerst een eigen functie in te voeren, bijvoorbeeld f (−1), f (−0.6), f (−0.2), . . . , f (1). Dat gaat als volgt. >> x = -1:0.4:1; >> z = vectorize(y) z = 2.*x.*exp(-x)-sin(x) >> eval(z) ans = -4.5951
-1.6219
Een controle is >> syms x >> subs(y,x,0.6) ans = 0.0939
-0.2899
0.1288
0.0939
-0.1057
8. OPGAVEN
49
De be¨ınvloeding van symbolen, strings en getallen onderling maakt het gebruik van MATLAB moeilijk. Door veel oefening krijgt u er op den duur toch gevoel voor. We geven ´e´en voor voorbeeld. >> syms x >> 2.3456 + x ans = 1466/625+x De omschrijving van het numerieke waarde 2.3456 naar de exacte breuk 1466/625 doet vermoeden dat het laatste resultaat in zijn geheel een symbolische waarde geworden is, want in MATLAB zelf wordt iedere exacte breuk meteen omgezet naar een numerieke waarde. Een controle met class leert dat dit inderdaad zo is.
8
Opgaven
Een wezenlijk onderdeel van het maken van de volgende opgaven is het op papier zetten van de uitwerkingen en de resultaten. Schrijf de essenti¨ele stappen op en niet iedere toetsaanslag. Symboliseer eventueel de benodigde variabelen met syms. Oefening 4.4 Teken de grafieken van de volgende functies op het interval [0, 2]. (a)
f (x) =
arctan(x) . 1 + x2
(b)
f (x) =
1 + sin(x) . (1 + cos(x))2
Oefening 4.5 Voer de volgende opdrachten letterlijk in en verklaar het verschil. (a)
x ∧ 1/3 .
(b)
x ∧ (1/3) .
Oefening 4.6 Substitueer met MATLAB in de uitdrukking x2 y 3 + sin(πxy) voor x en y respectievelijk 2 en 3. Z Oefening 4.7 Bepaal met MATLAB de onbepaalde integraal Laat met MATLAB zien dat dit antwoord correct is. Hint: Gebruik de MATLABfunctie numden. Oefening 4.8 Bepaal onderstaande integralen met MATLAB. (a)
Rπ
(b)
R∞
0
0
x/2 sin(x) dx. e−2x cos(3x) dx.
1 dx. 1 + x4
50
HOOFDSTUK 4. DE SYMBOLISCHE GEREEDSCHAPSKIST
(c)
R1
(d)
R4
0
log10 x dx.
|(x2 − 2)3 | dx. Vraag de numerieke waarde van het resultaat op. Bereken deze integraal door het integratie-interval in twee stukken te splitsen en de absolute waarden weg te werken. Zijn de resultaten met elkaar in overeenstemming? 0
Oefening 4.9 Bepaal met MATLAB de volgende twee integralen en leg uit wat het resultaat betekent: √ R2 R −t2 2 3 a. b. e dt. 0 sin(t ) 1 + t dt. Oefening 4.10 Teken met MATLAB de grafiek van f (x) = 3+sin(x) op het interval [−1, 2]. Zorg ervoor dat de oorsprong (0, 0) zichtbaar is. Oefening 4.11 Bekijk het resultaat van de MATLABopdrachten >> syms x >> ezplot(asin(sqrt(x)),[-5,5]) Waarom is de grafiek alleen voor 0 ≤ x ≤ 1 getekend? Maak met MATLAB een betere tekening van de grafiek.
Oefening 4.12 Bereken na uitvoering van de opdracht >> syms x n Z met MATLAB de integraal
xn dx .
(a) Voor welke n is dit antwoord incorrect? (b) Geef aan n de waarde uit onderdeel (a) en bereken de integraal opnieuw. Is het antwoord nu correct? Oefening 4.13 Verklaar het resultaat van de opdracht sqrt(x^2).
Oefening 4.14 Beschouw de functie f (x) = (log(1 + x) + x2 /2)/(1 + x)2 . Voer de volgende MATLABopdrachten uit >> syms x y >> y = (log(1+x)+x^2/2)/(1+x)^2 Beantwoord de volgende vragen:
8. OPGAVEN
51
(a) Bereken f 0 (x). (b) Bereken f 0 (1). (c) Teken de grafieken van f (x) en f 0 (x) in een plaatje.
Oefening 4.15 Gegeven is de functie f (x) = e2x − 5ex + 6. (a) Teken met MATLAB de grafiek van f op een geschikt interval. (b) Bepaal numeriek met fmin de plaats van het minimum van de functie f . (c) Bereken de functiewaarde en de waarde van de afgeleide in het door u gevonden punt. (d) Hoeveel nulpunten heeft f ?
p Oefening 4.16 Beschouw de functie f (x) = x − 2 + 8 3 (x − 1)2 . (a) Teken de grafiek (eerst even uit de hand) van f . (b) Bepaal de nulpunten en de extrema van f . (c) Laat zien dat u alle nulpunten en extrema gevonden heeft.
Oefening 4.17 Beschouw de functie f (x) = x sin(x). (a) Bepaal de extrema van de functie f , die binnen het interval (0, 2π) liggen. (b) Teken de grafieken van f (x) en sin(x) voor 0 ≤ x ≤ 2π in een plaatje. (c) Verklaar de ligging van de extrema van f ten opzichte van die van de sinus.
Oefening 4.18 Beschouw de vergelijking √
x sin
1 1 = , x>0. x 4
(a) Schets de functie
√
x sin( x1 ) en de constante functie
1 4
in een plaatje.
(b) Waarom heeft de vergelijking voor x dicht bij 0 geen oplossingen? (c) Waarom heeft de vergelijking voor grote x geen oplossingen? (d) Hoeveel oplossingen zijn er? (e) Benader op z’n minst een oplossing.
Oefening 4.19 De vergelijking ex = ax(1 − x) in de onbekende x heeft voor ´e´en positieve a precies ´e´en oplossing tussen 0 en 1. Geef met behulp van plaatjes een schatting voor deze a. Opmerking: Het getal a is exact te berekenen.
52
HOOFDSTUK 4. DE SYMBOLISCHE GEREEDSCHAPSKIST
√ 2 Oefening 4.20 De functie f (x) = e−x (x+ x) heeft ´e´en maximum in de buurt van x = 0.6. (a) Bepaal (een benadering van) de plaats x0 van het maximum. (b) Bereken f 0 (x0 ).
Oefening 4.21 Bepaal met MATLAB f (10) (x), x 6= 1, van de functie f (x) = sin( Het antwoord past op ´e´en regel! Verklaar het resultaat.
x2 + 2x − 3 ). x−1
Hoofdstuk 5
Lineaire Algebra 1
Inleiding
MATLAB is een programmapakket voor het uitvoeren van matrixberekeningen. Een matrix is in dit verband niets anders dan een rechthoekig schema van getallen dat we ook kunnen opvatten als een array. De naam MATLAB is een afkorting van matrix laboratory. Er bestaan verschillende versies van MATLAB. Bij deze handleiding zal worden uitgegaan van MATLAB voor Windows versie zoals ge¨ınstalleerd op uw notebook, inclusief de ’Extended Symbolic Math Toolbox’. Zoals uit de naam MATLAB blijkt, is de basisrekengrootheid een matrix. MATLAB kent voor het rekenen met matrices twee soorten bewerkingen; de matrixbewerkingen, zoals deze een rol spelen in de lineaire algebra, en de array-bewerkingen, waarmee u al heeft kennis gemaakt in de training MATLAB in het eerste trimester. De array-bewerkingen onderscheiden zich van de in de lineaire algebra gebruikelijke matrixbewerkingen doordat er een punt voor het bewerkingsteken wordt geplaatst. MATLAB kent ook matrixfuncties. Matrixfuncties zijn veelal opdrachten die, voor een gegeven matrix, binnen de lineaire algebra gebruikelijke berekeningen uitvoeren. Een aantal van deze functies zal binnen de cursus Lineaire Algebra worden behandeld.
2 2.1
Matrices Wat is een matrix?
Zij m, n ∈ IN . Een m × n matrix is een rechthoekig schema van getallen, gerangschikt in m rijen en n kolommen: a11 a12 . . . a1n a21 a22 . . . a2n .. .. .. .. . . . . am1 am2 . . .
amn
De getallen aij (1 ≤ i ≤ m, 1 ≤ j ≤ n) heten de elementen van de matrix. In deze notatie heet i de rij-index, j de kolom-index. We schrijven i.h.a. bovenstaande matrix kort op als de 53
54
HOOFDSTUK 5. LINEAIRE ALGEBRA
matrix A, of ook wel (aij ). MATLAB kent in wezen maar ´e´en rekengrootheid, namelijk een m × n-matrix (een matrix met m rijen en n kolommen met m ≥ 1 en n ≥ 1). Zo wordt een scalar in MATLAB opgevat als een 1 × 1-matrix.
2.2
Het invoeren van matrices
Een matrix kan op verschillende manieren ingevoerd worden, namelijk: 1. tussen [ en ], waarbij de rijen gescheiden worden door ; en de elementen gescheiden worden door spaties of komma’s. Bijvoorbeeld >>
A = [1 2 3;4 5 6;7 8 9]
of >>
A = [1,2,3;4,5,6;7,8,9]
2. tussen [ en ], waarbij elke rij op een nieuwe regel begint, b.v.: >>
A = [1 4 7
2 5 8
3 6 9]
3. via externe files, waarin de waarden van de variabelen zijn opgeslagen of waarin deze berekend worden. Een matrix element mag uit alle MATLAB uitdrukkingen bestaan. Zo is b.v.: >>
a = [-1.3
sqrt(3)
(1+2+3) * 4/5]
een toegelaten opdracht. We kunnen ook matrices samenvoegen tot een nieuwe matrix. Zo geeft C=[A;a] de matrix A met als vierde rij daaraan toegevoegd de matrix a. De formaten van de matrices moeten hiervoor wel bij elkaar aansluiten. Oefening 5.1 Voer de volgende 4 bij 4 matrix in. 1 2 −1 3 7 2 0 1 A = 3 −2 −1 −1 0 1 4 8
2. MATRICES
55
Oefening 5.2 Voer de volgende 2 bij 3 matrix in. 1.5 2 − 17 √ B = 3 2 0.625 Oefening 5.3 Voer de volgende matrices in
A=
1 2 3 4
, B=
5 6 7 8
1 1 , C= 1 1
2 2 2 2
Maak nu een matrix D met in de eerste twee kolommen de matrices A en B boven elkaar, en als 3e en 4e kolommen de matrx C.
2.3
Matrices met symbolische elementen
Het commando ‘sym’ construeert symbolische getallen, variabelen en objecten. Lees wat ”help sym”je over dit commando te vertellen heeft. Met het commando ‘syms’ kunnen meerdere symbolische objecten gelijkijdig geconstrueerd worden. De namen voor de objecten moeten in dit geval met een letter beginnen. Symbolische objecten kunnen als matrixelementen optreden. >> syms a b >> A=[ 1 2 a 3 b 4 1 2 3] A = [ 1, 2, a] [ 3, b, 4] [ 1, 2, 3] Oefening 5.4 Voer de volgende 2 bij 3 matrix in. a11 a12 a13 A = a21 a22 a23
2.4
Het invoeren van vectoren
Een rijvector is in feite een matrix met ´e´en rij, en een kolomvector een matrix met ´e´en kolom. Vectoren worden dus net zo ingevoerd als matrices. Als het verschil tussen de opeenvolgende elementen van de rij steeds hetzelfde is, kan de rij ook gevormd worden door ‘variabele’=‘beginwaarde’ : ‘stapgrootte’ : ‘eindwaarde’
56
HOOFDSTUK 5. LINEAIRE ALGEBRA
Hierbij is ‘beginwaarde’ het eerste element van de rij, ‘eindwaarde’ het laatse element van de rij en ‘stapgrootte’ het verschil tussen de opeenvolgende elementen. Bijvoorbeeld >> x=1:-.2:0 geeft x = 1.0000
0.8000
0.6000
0.4000
0.2000
0
Als het verschli tussen de opeenvolgende elementen van de rij ´e´en is kan de ‘stapgrootte’ worden weggelaten. Bijvoorbeeld >> x=1:5 geeft als resultaat x = 1
2
3
4
5
Voor een kolomvector is het vaak het makkelijkst om eerst een rijvector te maken en deze vervolgens te transponeren (rijen als kolommen schrijven, zie §2.7). ’ geeft in MATLAB de getransponeerde van een matrix aan. Bijvoorbeeld >> x=[1:5]’ geeft als resultaat de getransponeerde van de rij [1 2 3 4 5]. D.w.z. x = 1 2 3 4 5 Oefening 5.5 Maak de kolomvector met beginwaarde -1, eindwaarde 9 en stapgrootte 0.5. Oefening 5.6 Maak de kolomvector met beginwaarde 9, eindwaarde 0 en stapgrootte 1.
2. MATRICES
57
Oefening 5.7 Maak de kolomvector met beginwaarde a, eindwaarde a+20 en stapgrootte 2. Merk op dat bij het transponeren het programma er van uit gaat dat de symbolen complexe waarden kunnen aannemen. (zie de opmerkingen bij het transponeren van matrices blz.60) Men kan dit oplossen door te transponeren met het commando x.0 (arraybewerking). (Een andere mogelijkheid is a als ’real’ te declareren. Zoek op hoe dit moet met de help-functie (‘help sym’ of ‘help syms’)).
2.5
Bijzondere matrices
Behalve door direkte invoer, kunnen enkele bijzonder matrices ook rechtstreeks gevormd worden met MATLAB-opdrachten. Deze matrices zijn: – Een m × n nulmatrix, d.w.z. een matrix waarin alle elementen de waarde nul hebben, wordt gevormd met de opdracht ‘zeros(m, n)’. – Een m × n matrix, d.w.z. een matrix waarin alle elementen de waarde 1 hebben, wordt gevormd met de opdracht ‘ones(m, n)’. – Een m×n matrix met willekeurige getallen wordt gevormd met de opdracht ‘rand(m, n)’. – Een n × n diagonaalmatrix, d.w.z. een matrix met alleen op de diagonaal (de elementen met gelijke rij- en kolomindex) elementen ongelijk aan nul, wordt gevormd met de opdracht ‘diag(v)’, waarbij v een vector is met de n diagonaalelementen. – Een m×n eenheidsmatrix, d.w.z. een diagonaalmatrix met enen op de diagonaal, wordt gevormd met de opdracht ‘eye(m, n)’ Een vierkante m × n-matrix wordt gevormd als alleen m als index gebruikt wordt. Bijvoorbeeld een 3 × 3-eenheidsmatrix wordt gevormd met de opdracht: >>
eye(3)
Wil men een matrix met dezelfde dimensies als een gegeven matrix A dan kan dit door gebruik te maken van het commando ‘size’. Het commando >>
zeros(size(A))
vormt een nulmatrix met dezelfde dimensies als A. Van deze commando’s accepteert alleen het commando ‘diag’ een matrix als argument. Wanneer bij het commando ‘diag’ een matrix A wordt gebruikt i.p.v. de vector v krijgt men als antwoord een kolomvector met de diagonaalelementen van A. Oefening 5.8 Voer met behulp van bovenstaande commando’s de volgende matrices in: a. De 6 bij 6 eenheidsmatrix. b. De 5 bij 10 nulmatrix.
58
HOOFDSTUK 5. LINEAIRE ALGEBRA
c. De 5 bij 15 matrix waarin alle elementen de waarde 1 hebben. d. Een willekeurige 5 bij 5 matrix. e. Een diagonaalmatrix met op de diagonaal de getallen 1 t/m 5.
2.6
Indices
Vaak is het gewenst om met onderdelen van matrices of met afzonderlijke elementen te werken. Ook dit is eenvoudig te realiseren in MATLAB. Elk element van een matrix is te bereiken met behulp van indices. De opdracht hiervoor heeft de algemene vorm A(m, n). Hierbij is A de matrix, waarvan het element met rij-index m en kolomindex n wordt gespecificeerd. Als A een vector is, is ´e´en n index voldoende. Bij de matrix A uit oefening 5.1 levert A(3, 2) de waarde van dit element, namelijk − 2. Hiermee is het o.a. mogelijk om ´e´en afzonderlijk element van de matrix te veranderen. De opdracht: >>
A(3,2)=3
verandert bij de matrix A de oorspronkelijke waarde − 2 in de nieuwe waarde 3. Alle andere elementen van A blijven hetzelfde. Men kan in plaats ´e´en element van A ook meerdere elementen van A tegelijk aangeven. Dit doet men door in plaats van een losse index een vector van indices mee te geven. >>
A(3,[2 4])
Geeft de elementen a32 en a34 uit A gerangschikt als in A, d.w.z. een 1 × 2 matrix. Het commando >>
A([1 2 3],[2 4])
geeft dus de 3 × 2 matrix gevormd door de elementen van A in de eerste drie rijen en in de kolommen twee en vier. Gezien de opmerkingen uit §2.4 kunnen we dit laatse ook bewerkstelligen met het commando >>
A(1:3,[2 4])
Wanneer alleen ‘:’ wordt gebruikt, worden alle rijen of kolommen bedoeld. Dit kan in de volgende vormen voorkomen: – A(:, j)
is de j-de kolom van A.
– A(i, :)
is de i-de rij van A.
– A(:, :)
is hetzelfde als A.
2. MATRICES
59
Met name wanneer je met hele rijen of kolommen operaties uit wil voeren is deze mogelijkheid goed te gebruiken. Zo levert, bij de oorspronkelijk matrix A het commando: >>
A(2,:)-7*A(1,:)
als resultaat de rijvektor op, die ontstaat door zeven maal de eerste rij van de tweede rij af te halen. Met het commando: >>
A(2,:)=A(2,:)-7*A(1,:)
wordt hetzelfde resultaat berekend maar dit wordt nu toegekend aan de tweede rij van de matrix A. De matrix A is dus door dit commando veranderd. Oefening 5.9 Voer bovenstaande commando’s uit met de matrix A gegeven in oefening 5.1. Oefening 5.10 Ga na wat het effect is van het commando >>
A([1 3],[1 3])=10*ones(2)
op de matrix A gegeven in oefening 5.1.
2.7
Matrixbewerkingen in MATLAB
Binnen MATLAB is het mogelijk om rekenkundige bewerkingen en opdrachten uit te voeren met matrices als geheel. Hiertoe kent MATLAB een groot aantal commando’s. Een aantal van deze commando’s zullen in de loop van de cursus Lineaire Algebra (zie [2]) nog besproken en gebruikt worden. Meer over de wiskundige context vindt men in [3]. Voor dit moment zijn de volgende operaties binnen MATLAB van belang: vermenigvuldigen: Het vermenigvuldigen van twee matrices gebeurt op speciale wijze. Voor een m×n matrix A = (aij ) en een n × p matrix B = (bjk ) defini¨ Pn eren we het product AB van A met B als de m × p matrix (cik ) waar cik = j=1 aij bjk . Dus: ... b .. .. .. .. 1k . . . . . . . . . . b2k . . . . . . cik . . . ai1 ai2 . . . ain = . .. .. .. .. .. . . . . . . . bnk . . . Het product van de matrices A en B kan alleen gevormd worden als de formaten (dimensies) van A en B passend zijn: het aantal kolommen van A moet gelijk zijn aan het aantal rijen van B. Als AB bestaat hoeft BA niet te bestaan. Als zowel AB als BA bestaan geldt in het algemeen niet AB = BA. Het matrixproduct van A en B is binnen MATLAB altijd gedefinieerd als A of B een getal is. Het vermenigvuldigen van een matrix met een getal (scalaire vermenigvuldiging) komt overeen met het vermenigvuldigen van alle matrixelementen met dat getal. Voor matrixvermenigvuldiging gebruikt men in MATLAB het symbool ‘*’ dus
60
HOOFDSTUK 5. LINEAIRE ALGEBRA
>> A*B optellen en aftrekken: Het optellen >> A+B en aftrekken >> A-B van twee matrices A en B komt overeen met het optellen en aftrekken van de afzonderlijke elementen van de twee matrices. Deze bewerkingen zijn gedefinieerd als de dimensies van de matrices gelijk zijn of als een van de twee matrices een getal is. In dat laatste geval wordt het getal bij elk element van de matrix opgeteld of afgetrokken. machtsverheffen: Het commando >> A^p verheft de matrix A tot de p-de macht. Als p een positief geheel getal is, dan wordt A ∧ p berekent door herhaald vermenigvuldigen van A met zichzelf. transponeren: Zij A = (aij ) een m × n matrix. dan is de getransponeerde AT van A de n × m matrix (bij ) (1 ≤ i ≤ n, 1 ≤ j ≤ m) met bij = aji . Dus: T a11 . . . a1n a11 . . . am1 .. = .. .. AT = ... . . . am1 . . .
amn
a1n . . .
amn
In MATLAB is A0 de getransponeerde van de matrix A. Bijvoorbeeld: >>
[1 2 3]’
vormt de kolommatrix: >>
1 2 3
Als: >>
x=[1 2 3]
wordt dezelfde matrix verkregen met de opdracht: >>
x’
Merk op dat indien een matrix A complexe, niet re¨ele, elementen bevat het commando A0 de matrix niet alleen spiegelt maar ook van ieder element de complex geconjugeerde neemt. Wil men alleen de gespiegelde dan kan dit met het commando A.0 Ter illustratie een tweetal voorbeelden.
2. MATRICES
61
>> A=[2+3*i 4+5*i 2 3] A = 2.0000 + 3.0000i 2.0000
4.0000 + 5.0000i 3.0000
>> A’ ans = 2.0000 - 3.0000i 4.0000 - 5.0000i
2.0000 3.0000
>> A.’ ans = 2.0000 + 3.0000i 4.0000 + 5.0000i
2.0000 3.0000
Hiermee dient men bij het gebruik van symbolische matrixelementen rekening te houden. >> A=[a b 1 2] A = [ a, b] [ 1, 2] >> A’ ans = [ conj(a), [ conj(b),
1] 2]
>> A.’ ans = [ a, 1] [ b, 2] De comando’s ‘\’ en ‘/’: De ‘deling’ is bij matrices niet gedefinieerd. In MATLAB wordt echter bij matrices
62
HOOFDSTUK 5. LINEAIRE ALGEBRA
wel de deelstreepnotatie gebruikt. We zullen spreken van een symbolische deling (binnen MATLAB ook aangeduid met ‘mldivide’ en ‘mrdivide’ ). De betekenis is alsvolgt >> X=A\B geeft de oplossing X van de matrixvergelijking A ∗ X = B, en >> X=B/A geeft de oplossing X van de matrixvergelijking X ∗ A = B. Omdat AX in het algemeen ongelijk is aan XA moet er onderscheid gemaakt worden tussen een symbolische deling aan de ‘rechterkant’ en een symbolische deling aan de ‘linkerkant’. Wanneer de op te lossen vergelijking meerdere oplossingen heeft geeft MATLAB er slechts ´e´en. Wanneer de vergelijking geen oplossingen heeft geeft MATLAB de ‘kleinste kwadraten oplossing’. Het verdient dus aanbeveling het antwoord van MATLAB te controleren. We komen hierop terug in §3 Meerinformatie kunt u vinden in bijv. de “helpdesk”door onder “index”te kijken bij “arithmetic operators”. Oefening 5.11 Bepaal van de eerder ingevoerde matrix A (Oefening 5.1) de getransponeerde matrix AT d.m.v. het commando: B=A’. Bereken AAT via: C = A ∗ B. Ga na dat C symmetrisch is. Bereken AT A via: D = B ∗ A. Vergelijk het resultaat met C = A ∗ B. Bereken 3A + 5B 3 via: E = 3 ∗ A + 5 ∗ B ∧ 3. Oefening 5.12 a. Voer de vector v = (1 2 3 4) in. Bereken vervolgens v ∗ A en A ∗ v 0 , met A de matrix uit de voorgaande oefening. Wat gebeurt er door deze vermenigvuldigingen met de rijen en kolommen van A? b. Maak een diagonaalmatrix D met op de diagonaal de elementen van de vector v. Bereken D ∗ A en A ∗ D. Wat gebeurt er door deze vermenigvuldigingen met de rijen en kolommen van A? Oefening 5.13 Bereken met de vector v uit de vorige oefening v ∗ v 0 en v 0 ∗ v. Ga na dat er sprake is van ’gewone’ matrixvermenigvuldiging. Oefening 5.14 Herhaal de opdrachten uit oefening 5.11 voor de matrix 1 2 a A = b −3 2 . 3 1 c
2.8
Matrixfuncties in MATLAB
In deze paragraaf noemen we een aantal commando’s die een matrix als argument hebben (Denk b.v. ook aan ‘zeros(A)’, blz. 57). Bij de training hoeft u de commando’s in deze
3. OPLOSSEN VAN STELSELS LINEAIRE VERGELIJKINGEN
63
paragraaf niet uit te voeren. De meeste van deze commando’s zullen worden gebruikt in het college Lineaire Algebra (zie [2]) waar ook duidelijk zal worden wat de genoemde begrippen precies betekenen. Op dit moment gaan we hier niet verder op in. Het commando >> inv(A) berekent de inverse van een vierkante matrix A als A inverteerbaar is, d.w.z. berekent een matrix B zodat A ∗ B = B ∗ A = I. Als er geen inverse bestaat (A singulier) of als het numerieke resultaat niet betrouwbaar is zal MATLAB dit melden. ([2, week 2]). Het commando >> det(A) berekent de determinant van A. A moet een vierkante matrix zijn. ([2, week 2]). Het commando >> rank(A) geeft de rang van de matrix A. ([2, week 3]). Het commando >> eig(A) geeft een vector met daarin de eigenwaarden van de A. De matrix A moet vierkant zijn. Dit wordt vaak gebruikt in de vorm >> [S,D]=eig(A) Dit geeft als resultaat een matrix S en een matrix D, zodat A ∗ S = S ∗ D. Iedere kolom van S een eigenvector en D een diagonaalmatrix met op de diagonaal de eigenwaarden van A behorend bij de eigenvectoren in de overeeenkomstige kolommen van S. ([2, week 6]). Opmerking: Bovenstaande commando’s werken ook met symbolische matrices. I.h.a. moet men er op bedacht zijn dat matrixfuncties mogelijkerwijs niet werken voor symbolische matrices. Omgekeerd zijn er ook commando’s die alleen werken voor symbolische matrices.
3 3.1
Oplossen van stelsels lineaire vergelijkingen Het commando A \ b
Binnen MATLAB is het mogelijk om rechtstreeks de oplossing van een stelsel vergelijkingen te bepalen. Het commando om Ax = b rechtstreeks op de lossen luidt:
64
>>
HOOFDSTUK 5. LINEAIRE ALGEBRA
x=A\b
(Binnen MATLAB kan men meer informatie vinden via het trefwoord ”mldivide”.) Soms geeft MATLAB, na de opdracht x = A\b, de volgende waarschuwing. >>
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = .... .
We gaan op de precieze betekenis van deze waarschuwing hier niet in. De waarde van x die op het scherm verschijnt, moeten we niet vertrouwen. Een andere mogelijke waarschuwing is >>
Warning: Matrix is singular to working precision.
Dus: Verschijnt er een van bovenstaande waarschuwingen, dan is de waarde van x onbetrouwbaar.
Indien deze meldingen niet worden gedaan, wil het nog niet zeggen dat A∗x = b. I.h.a. zal dit commando het stelsel oplossen met een Gauss-algoritme. Echter, dit commando geeft soms ook een antwoord x wanneer een stelsel vergelijkingen geen (exacte) oplossing heeft. In dat geval wordt een ander algoritme gebruikt. De betekenis van deze ’kleinste-kwadratenoplossing’ x komt in de cursus Lineaire algebra nog aan bod ([2, week 4]). Het is daarom belangrijk een antwoord dat je met dit commando bepaald hebt expliciet te controleren door A ∗ x uit te rekenen en te vergelijken met de gewenste uitkomst b. Geldt niet A ∗ x = b, dan is het stelsel overbepaald. Geldt wel A ∗ x = b dan is het stelsel regulier of onderbepaald. MATLAB levert namelijk hoogstens ´e´en oplossing. Bij onderbepaalde stelsels, d.w.z. stelsels met meerdere oplossingen, moet de algemene oplossing dan anders bepaald worden, bijv. door vegen van de matrix. Bovenstaande methode werkt ook indien A en b symbolische matrices zijn. Men krijgt een waarschuwing als de oplossing niet bestaat.
3.2
De rijgereduceerde trapvorm
MATLAB kent een commando om een matrix te vegen met het Gauss-Jordan algoritme, d.w.z. om de rijgereduceerde trapvorm van een matrix te bepalen. Binnen MATLAB kunnen we de rijgereduceerde trapvorm (reduced row echelon form) van een (symbolische) matrix A bepalen met het commando >> rref(A) Passen we dit commando toe op de uitgebreide co¨effici¨entenmatrix van een stelsel vergelijkingen dan kunnen we uit het resultaat de oplossing bepalen.
3. OPLOSSEN VAN STELSELS LINEAIRE VERGELIJKINGEN
65
Oefening 5.15 Beschouw de matrix 1 2 3 6 2 en de vector b = 14 A = 2 −3 3 1 −1 −2 en het stelsel Ax = b. Los dit stelsel op met het commando A\b en met het rref commando. Oefening 5.16 Beschouw de matrix −4 1 2 −3 A= en de vector b = 2 1 −3 4 en het stelsel Ax = b. Los dit stelsel op met commando A\b. Is de oplossing uniek? Bepaal alle oplossingen door gebruik te maken van rref([A,b]). Maak symbolische matrices AA = sym(A) en bb = sym(b) en los het stelsel ook symbolisch op met het commando AA\bb. Oefening 5.17 1 1 A= 1
Beschouw de matrix 10 0 2 en de vector b = 16 17 4
en het stelsel Ax = b. Ga met rref([A,b]) na dat dit stelsel strijdig is. Wat geeft het commando A\b? Is dit een oplossing? Maak symbolische matrices AA = sym(A) en bb = sym(b) en kijk wat het resultaat is van het commando AA\bb.
3.3
Het oplossen van stelsels met de Symbolic Toolbox
De symbolic toolbox van MATLAB maakt het mogelijk ook vergelijkingen op te lossen waarin de co¨effici¨enten symbolische uitdrukkingen zijn. Het commando waarmee symbolische vergelijkingen (niet noodzakelijk lineair) kunnen worden opgelost is het commando ’solve’. Dit commando wordt alsvolgt gebruikt solve(’expressie1’,’expressie2’,’onbekende1’,’onbekende2’) Hierin zijn expressie1 en expressie2 de vergelijkingen, en onbekende1 en onbekende2 de variabelen waarnaar men wil oplossen. Voorbeeld: >> syms a b c d p q >> [x1,x2]=solve(’a*x1 +b*x2=p’,’c*x1 +d*x2=q’,’x1’,’x2’) x1 =
66
HOOFDSTUK 5. LINEAIRE ALGEBRA
(-b*q+p*d)/(-b*c+d*a)
x2 = (-c*p+q*a)/(-b*c+d*a) Meer detailinformatie over dit commando is te verkrijgen via de online-help van het programma. Oefening 5.18 Beschouw het stelsel vergelijkingen uit oefening 5.16 en los dit stelsel op met behulp van het commando ‘solve’. Hoe ziet de algemene oplossing eruit in vectornotatie? Oefening 5.19 Maak BK opgave 1.5.11 door de oplossing via symbolische matrices te bepalen in termen van de parameter a. Is de oplossing voor alle a gedefinieerd? Beschouw de waarden van a waarvoor deze oplossing niet is gedefinieerd apart, bijv. door gebruik te maken van het commando voor substitutie ‘subs’. Oefening 5.20 Herhaal de opdracht uit oefening 5.19 voor BK opgave 1.5.13. Onderzoek √ of er oplossingen zijn voor a = 6 via ‘rref’ en via ‘A\b’.
4
Numerieke aspecten bij gebruik van MATLAB
Wanneer we van een pakket als MATLAB gebruik maken dan moeten we er rekening mee houden dat het een numeriek pakket is. Bij numerieke algoritmen spelen de aspecten van effici¨entie, afrondfouten en datafouten een rol. We zullen ons beperken tot een korte illustratie. Deze onderwerpen komen in latere colleges in meer detail aan bod. Opmerking. Wanneer men binnen de symbolic toolbox numeriek rekent maakt men geen gebruik van de numerieke algoritmen van MATLAB maar van de algoritmen van Maple.
4.1
Getalrepresentatie en afrondfouten, een voorbeeld
Alle waarden waarmee gerekend wordt worden gerepresenteerd met een beperkt aantal cijfers. De computer rond alle re¨ele getallen af tot een getal dat past in de gebruikte representatie. Voorbeeld (standaard outputrepresentatie in MATLAB): 1 3 = 0.3333 √ 2 = 1.4142 −10 e = 4.5400e-005 e10 = 2.2026e+004 1 10 = 0.1000 1 = 1.0000e-006 106
4. NUMERIEKE ASPECTEN BIJ GEBRUIK VAN MATLAB
67
Intern rekent MATLAB met een grotere precisie dan op het scherm wordt weergegeven. Dat de outputrepresentatie tot onverwachte uitkomsten aanleiding kan geven wordt duidelijk in de volgende oefening. Oefening 5.21 Beschouw de matrix 9 0 A= 0 1 Bereken met de hand A2 en A3 . Kunt u een formule geven voor An ? Bereken nu met MATLAB A10 . Is dit antwoord exact? Wanneer we een antwoord met meer cijfers willen representeren kan dit m.b.v. het commando >> format long Geef het commando ’format long’ en daarna opnieuw de opdracht A10 . vergelijk de uitkomst met het eerder verkregen antwoord. Geef vervolgens het commando >> format om weer terug te keren naar een representatie met 5 cijfers. We zien dus dat ’afrondfouten’ gevolgen kunnen hebben voor het resultaat. Gebruik ’help format’ voor meer informatie over de mogelijke representaties van matlab output. De rol die afrondfouten in berekeningen spelen kan vaak beperkt worden door het slim inrichten van de algoritmen.
4.2
Datafouten, een voorbeeld
Andere fouten waar je mee te maken kunt krijgen zijn ’datafouten’. In de volgende oefening zullen we illustreren dat, bij het oplossen van een stelsel Ax = b, een kleine verandering in de data grote gevolgen kan hebben voor de oplossing. Hierbij spelen intrinsieke eigenschappen van A een rol. Oefening 5.22 (ML) Beschouw twee stelsels vergelijkingen, nl. 10x1 + 7x2 + 8x3 + 7x4 = 32 7x1 + 5x2 + 6x3 + 5x4 = 23 8x1 + 6x2 + 10x3 + 9x4 = 33 7x1 + 5x2 + 9x3 + 10x4 = 31 en 2x1 + x2 + 5x3 + x4 x1 + x2 − 3x3 − 4x4 3x1 + 6x2 − 2x3 + x4 2x1 + 2x2 + 2x3 − 3x4
= 9 = −5 = 8 = 3
68
HOOFDSTUK 5. LINEAIRE ALGEBRA
a. Bepaal van beide stelsels de oplossingen. Door bij deze stelsels de rechterleden iets te veranderen en te kijken naar het effect daarvan op de oplossing is het mogelijk om een indruk te krijgen hoe gevoelig de matrices zijn voor kleine veranderingen in de data. b. Los het eerste stelsel ook op met rechterleden (32.1, 22.9, 32.9, 31.1) en (32.01, 22.99, 32.99, 31.01). Los het tweede stelsel ook op met rechterleden (9.1, −5.1, 7.9, 3.1) en (9.01, −5.01, 7.99, 3.01). Ga na wat het effect van deze veranderingen op de oplossing is. Door de co¨effici¨enten van de vergelijkingen iets te veranderen en te kijken naar het effect daarvan op de oplossing is het ook mogelijk om een indruk te krijgen van de gevoeligheid voor kleine veranderingen in de data. c. Verander de elementen van de oorspronkelijke co¨effici¨enten–matrices door er de matrix 0.1 ∗ rand(4) bij op te tellen (in dat geval wordt bij elk element een willekeurig getal tussen 0 en 0.1 opgeteld). Bepaal nu de oplossingen van deze aangepaste stelsels bij de oorspronkelijke rechterleden, respectievelijk (32, 23, 33, 31) en (9, −5, 8, 3). En ga na wat het effect van deze veranderingen op de oplossing is. d. Bij stelsels die gevoelig zijn voor kleine veranderingen in de data noemen we de co¨effici¨entenmatrix slecht-geconditioneerd. Welk van beide stelsels vergelijkingen heeft de best geconditioneerde co¨effici¨entenmatrix? Oefening 5.23 Met de MATLAB opdracht A = hilb(n) wordt A de n bij n Hilbert matrix, gedefinieerd door aij = 1/(i + j − 1). a. Bepaal de 5 bij 5 Hilbert matrix. b. Los op Ax1 = b1 = (2.2833, 1.4500, 1.0929, 0.8845, 0.7456)T . c. Los op Ax2 = b2 = (2.2834, 1.4501, 1.0928, 0.8844, 0.7457)T . d. Vergelijk b1 met b2 en x1 met x2 . De matrix A is blijkbaar slecht geconditioneerd.
4.3
Effici¨ entie
Het inrichten van algoritmen dient natuurlijk zo effici¨ent mogelijk te gebeuren. Als maat voor de effici¨entie van een algoritme wordt wel gebruikt het aantal bewerkingen (optellingen, vermenigvuldigingen etc.) dat voor het uitvoeren van het algoritme nodig is. Bij MATLAB heet een enkele bewerking een ’flop’ (floating point operation). Na afloop van iedere sessie geeft MATLAB het aantal gebruikte flops. De flop-teller kan op nul gezet worden met het commando
5. MATLAB EN LINEAIRE ALGEBRA
69
>> flops(0) Het aantal tot op dat moment gebruikte flops krijgt men met het commando >> flops Oefening 5.24 Hoeveel flops zijn er nodig voor het vermenigvuldigen van twee n × nmatrices? Controleer dit door bijv. twee random 5 × 5- en 10 × 10-matrices met elkaar te vermenigvuldigen. Zet voor iedere vermenigvuldiging de teller op nulen vraag na de vermenigvuldiging het aantal flops op. Oefening 5.25 Doe hetzelfde als in de voorgaande oefening maar nu voor het matrixproduct van een m × n matrix met een n × k-matrix. Oefening 5.26 Beschouw de matrix A11 ∅ A= ∅ A22 met A11 en A22 willekeurige 40 bij 40 matrices en ∅ een 40 bij 40 nulmatrix. Het berekenen van A ∗ A kan op twee manieren, namelijk rechtstreeks en gebruik makend van het feit dat A11 ∗ A11 ∅ A∗A=A= ∅ A22 ∗ A22 In het laatste geval is vermenigvuldiging van een nulmatrix met een andere matrix niet nodig. Bepaal voor beide manieren, m.b.v. MATLAB, het aantal benodigde flops en vergelijk de resultaten
5
MATLAB en Lineaire Algebra
We besluiten deze training met een aantal opdrachten die door toepassing van matrixrekening kunnen worden opgelost. MATLAB is hierbij het aangewezen gereedschap. Hoewel de opdrachten ook met de hand kunnen worden opgelost (eventueel ter controle) dient men bij iedere opdracht na te gaan hoe zoveel mogelijk van het rekenwerk op een effici¨ente manier met MATLAB uitgevoerd kan worden. Vaak zal de oplossing verkregen worden door het oplossen van een stelsel vergelijkingen in dat geval moeten alle mogelijke oplossingen bepaald worden.
5.1
Opdrachten
Oefening 5.27 Bepaal alle 2 × 2 matrices die commuteren met 1 0 . 0 3
70
HOOFDSTUK 5. LINEAIRE ALGEBRA
Oefening 5.28 Bepaal alle 2 × 2 matrices die commuteren met 1 3 . 0 2 Oefening 5.29 Beschouw n deeltjes in het vlak met plaatsvectoren r1 , r2 , · · · , rn en massa’s m1 , m2 , · · · , mn . Voor de plaatsvector rzp van het zwaartepunt van dit stelsel van n massa’s geldt dan rzp =
1 (r m1 + r2 m2 + · · · + rn mn ) , M 1
met M = m1 + m2 + · · · + mn . Beschouw een plaat in het vlak met hoekpunten (1,2), (5,1) en (4,4) waarvan de massa geheel geconcentreerd is in de hoekpunten. Hoe kunnen we een massa van 1 kg verdelen over de hoekpunten zodanig dat het zwaartepunt in het punt (3,3) ligt? Beschouw nu een plaat met vier hoekpunten waaronder in ieder geval de reeds gegeven drie punten (kies zelf een vierde punt). Wat zijn nu de mogelijke massaverdelingen? Oefening 5.30 De impulsvector P van een stelsel van n deeltjes met massa’s m1 , m2 , · · · , mn en snelheden v 1 , v 2 , · · · , v n is gegeven door P = m1 v 1 + m2 v 2 + · · · + mn v n . Beschouw twee deeltjes met snelheden v 1 = (2, 2, 2)T en v 2 = (5, 7, 9)T . De deeltjes botsen en hebben na botsing de snelheden w1 = (4, 6, 4)T en w2 = (4, 5, 8)T . Wat kunnen we, onder de veronderstelling van behoud van impuls, zeggen over de massa’s van de twee deeltjes? Oefening 5.31 Teken in het gebied −6 ≤ x ≤ 6, −6 ≤ y ≤ 6 de driehoek met hoekpunten 1 5 4 , p2 = , p3 = p1 = 2 1 3 Zij D de matrix 0 −1 D= 1 0 Teken de driehoek met hoekpunten bi = Dpi . Teken ook de driehoeken met hoekpunten D2 pi en D3 pi . Wat is het effect van de vermenigvuldiging met D? Met welke matrix moeten we vermenigvuldigen om te spiegelen in de x-as?
Hoofdstuk 6
Differentiaalvergelijkingen 1 1.1
M-files M-files maken
Tot nu toe is MATLAB hoofdzakelijk gebruikt in een commando-gestuurde mode; de commando’s worden regel voor regel ingevoerd, waarbij de commando’s op een regel onmiddelijk worden uitgevoerd. MATLAB kan echter ook een reeks van commando’s, die opgeslagen zijn in een file, uitvoeren. Dergelijke files heten M-files naar de extensie ‘.m’. Een M-file bestaat uit een serie opdrachten in de MATLAB taal. Andere programma’s in die taal kunnen een M-file aanroepen. Een M-file kan ook zichzelf aanroepen en zo een recursief proces vormen. Een M-file is een gewone ASCII-file die geschreven kan worden met iedere gewone teksteditor. Het ligt echter het meest voor de hand om de in MATLAB ingebouwde faciliteit te gebruiken. Deze wordt opgeroepen via het menu ’File/New/M-file’. Hiermee wordt de ’MATLAB Editor/Debugger’ geopend. Een M-file is een file met de extensie ’.m’. De file bevat een reeks matlabopdrachten. Ook alle standaard matlabcommando’s liggen vast in een file ’commando.m’. De opdrachten in de file commando.m worden uitgevoerd door binnen MATLAB de opdracht >> commando te geven.
1.2
Script- en functionfiles
Er bestaan twee soorten M-files, de zogenaamde scriptfiles en de functionfiles (zie ook §16). In de informatica is een script een commandoreeks die ge¨ınterpreteerd wordt in plaats van gecompileerd. Een scriptfile bestaat dus uit een aantal commandoregels die uitgevooerd 71
72
HOOFDSTUK 6. DIFFERENTIAALVERGELIJKINGEN
worden als commando’s die in de commando-gestuurde mode worden ingevoerd. De variabelen kunnen dus reeds eerder gedefinieerde variabelen zijn en nieuwe variabelen kunnen door de scriptfile worden ge¨ıntroducerd. De variabelen zijn zogenaamde globale variabelen. Bijvoorbeeld in de file ”rijom.m”, een file die de volgorde van de elementen in een rij omdraait, kunnen de volgende opdrachten voorkomen [r k]=size(x); y=x(k:-1:1) Deze scriptfile wordt aangeroepen met de opdracht ‘rijom’. Het programma neemt aan dat de variabele x bestaat. Na afloop van het programma staan de variabelen x, r, k en y in het geheugen. Een functionfile begint met een opdracht die de in- en uitvoer variabelen van de opdracht definieert. Een voorbeeld van een functionfile, die hetzelfde doet als bovenbeschreven scriptfile, is: function y=wissel(x) % De functionfile , % draait de volgorde van de elementen van de rij x om % en bergt het resultaat op in matrix y. [r k]=size(x); y=x(k:-1:1) De eerste regel begint met ‘function’ om aan te geven dat het niet om een ‘script’ gaat. De naam van de functie is ‘wissel’ en de naam van de file is wissel.m . De functie kent 1 parameter x. Deze variabelenaam is alleen van belang binnen de functiefile en bestaat niet buiten de functiefile. Het is een zogenaamde locale variabele. Het resultaat van de functie heet y. Ook y is een locale variabele. De variabelen die je in een functiefile gebruikt zijn dus afgeschermd van de andere variabelen en be¨ınvloeden de werking van de variabelen buiten de functie niet. Het %-teken geeft aan dat er een commentaarregel volgt. Het aaneengesloten commentaarblok na de eerste regel van een functionfile bevat doorgaans een beschrijving van de functie. Het commando ‘help wissel’ geeft de deze beschrijving op het scherm. De aanroep ‘a=wissel(b)’ zal nu de kolommen van een vooraf gedefinieerde matrix b in volgorde omwisselen en het resultaat opbergen in matrix a. Er is bij een functiefile sprake van locale variabelen. Wil men variabelen in een functiefile ook buiten de functiefile betekenis geven dan kan dit met het commando ‘global’ (zie ook help global). Oefening 6.1 Schrijf de bovenbeschreven scriptfile rijom.m Maak het MATLAB geheugen leeg met het commando ‘clear’ en definieer de rij x=[1 3 5 7 9 11] Voer vervolgens de scriptfile uit. Wat is de verwachte uitvoer van deze scriptfile ? Bekijk met het commando ‘who’ welke variabelen in het geheugen staan.
2. DIFFERENTIAALVERGELIJKINGEN
73
Oefening 6.2 Schrijf de bovenbeschreven functionfile wissel.m. Maak het MATLAB geheugen leeg met het commando ‘clear’ en definieer de rij a=[1 3 5 7 9 11]. Pas de functionfile toe op a en bekijk daarna met het commando ‘who’ welke variabelen in het geheugen staan.
2 2.1
Differentiaalvergelijkingen Het numeriek oplossen van differentiaalvergelijkingen
Voor het berekenen (eigenlijk ’benaderen’) van de oplossing van de differentiaalvergelijking wordt gebruik gemaakt van algoritmen die verwant zijn aan het Euler algoritme dat in het college Calculus is behandeld. Het Euler algoritme komt, in het kort, neer op het volgende. Hoewel men de oplossing van de differentiaalvergelijking niet expliciet kent, kent men wel de afgeleide van de oplossing in ieder punt (deze wordt gegeven door de diff.vgl.), m.a.w. men kent het richtingsveld. Dit richtingsveld wordt nu gebruikt om de oplossing te benaderen. Vanuit een gegeven beginwaarde wordt het volgende punt van de oplossing berekend door vanuit de beginwaarde een stap te zetten in de richting van het richtingsveld. Vanuit het zo verkregen punt doet men weer hetzelfde, etc. De benadering van de oplossing hangt dus af van de beginwaarde, het richtingsveld en de stapgrootte. In MATLAB kan men differentiaalvergelijkingen numeriek oplossen met behulp van de commando’s ‘ode23’ of ‘ode45’. De onderliggende algoritmen maken gebruik van Runge-KuttaFehlberg integratie met variabele stapgrootte, d.w.z. dat het algoritme de stapgrootte groter maakt wanneer de oplossing minder varieert. ‘ode23’ maakt gebruik van 2e en 3e orde formules, ‘ode45’ van 4e en 5e orde formules. Met ‘ode45’ wordt een hogere nauwkeurigheid bereikt. Beide commando’s zijn toepasbaar voor stelsels eerste orde differentiaalvergelijkingen van de vorm dy = f (t, y) dt Hierin is y de toestandsvector en f de vectorwaardige functie die de afgeleide van de toestand geeft als functie van t en y. De vorm waarin ‘ode23’ gebruikt wordt is: >> [t,x]=ode23(’filename’,[t0,tm],y0) Hierin is filename de naam van de .m file waar de differentiaalvergelijking is gedefinieerd, [t0,tm] het interval waarover ge¨ıntegreerd wordt en y0 de kolomvector met beginwaarden. Het algoritme in het ’programma’ ode23 bepaalt eerst in een punt de afgeleide (x˙ van de functie gedefinieerd in de functiefile ‘filename’. Hiermee wordt dan het volgende punt bepaald, etc. het algoritme bepaald zelf met welke stappen het interval [t0, tm] doorlopen wordt. Omdat ‘ode45’ met hogere orde formules werkt zijn hierbij veelal minder integratiestappen
74
HOOFDSTUK 6. DIFFERENTIAALVERGELIJKINGEN
nodig. Als gevolg hiervan geeft ‘ode45’ dus soms minder ‘gladde’ plaatjes dan ‘ode23’. De output van ‘ode’ geeft een vector t waarin de tijdstippen waarop de oplossing x is uitgerekend (de stappen waarmee het interval [t0, tm] doorlopen is) en een rij(vector) x waarin in de bijbehorende functiewaarden zijn uitgerekend. We zullen e.e.a. illustreren met een voorbeeld waarin ‘ode23’ wordt gebruikt. Voor ‘ode45’ gaat alles analoog. Voorbeeld 6.1 Beschouw de differentiaalvergelijking x˙ = 1 − tx Om te beginnen maken we de file diffvgl.m met als inhoud function xdot=diffvgl(t,x) xdot=1-t*x; De punt-komma voorkomt dat de output op het scherm komt. De opdracht >> [t,x]=ode23(’diffvgl’,[0,1],3) berekent nu voor het tijdsinterval [0, 1] waarden van x(t) en met beginwaarde x(0) = 3. Men kan de oplossingen plotten met >> plot(t,x)
In bovenstaande wordt de output van de ode-solver toegekend aan de vextoren t en x. Er bestaat echter ook de eigenschap ‘OutputFcn’ waarmee men de ode-solver kan aangeven wat er met de output moet gebeuren. De standaard (default) OutputFcn is ‘odeplot’. D.w.z. dat de berekende waarden worden geplot zodra zij zijn berekend, de plot wordt voortdurend aangepast. Het commando >> ode23(’diffvgl’,[0,1],3) levert dus uiteindelijke hetzelfde resultaat als ’plot(t,x)’. Oefening 6.3 Bereken met ‘ode23’ de oplossing van de differentiaalvergelijking x˙ = 1 − 2tx met beginwaarde x(0)=1, op het interval [0,3]. Plot het resultaat. Doe dit op de manier van voorbeeld 6.1. Wat is het resultaat indien men de toekenning [t,x]= weglaat?
2. DIFFERENTIAALVERGELIJKINGEN
75
Oefening 6.4 Beschouw het beginwaardeprobleem x˙ =
tx , x(0) = 3 . (1 + x2 )
a. Bereken met ode23 x(t) op het interval [0,2]. Hoeveel stappen zijn er gebruikt? Plot de grafiek. b. Wat is de benadering voor x(2) in onderdeel a. ?
2.2
Stelsels differentiaalvergelijkingen
De opdrachten ‘ode23’ en ‘ode45’ kunnen ook in veel algemenere situaties worden toegepast. Voorbeeld 6.2 Beschouw het stelsel differentiaalvergelijkingen x˙ 1 = 5x1 − 2x2 x˙ 2 = 7x1 − 4x2 door de beginwaarden x1 (0) = 2 en x2 (0) = 8 wordt precies ´e´en oplossing vastgelegd. Met de opdrachten ‘ode23’ en ‘ode45’ kunnen dergelijke differentiaalvergelijkingen worden aangepakt. Om te beginnen maken we de file diffvgl.m met als inhoud function xdot=diffvgl(t,x) xdot=[5*x(1)-2*x(2);7*x(1)-4*x(2)]; Merk op dat xdot een kolomvector moet zijn. De opdracht >> [t,x]=ode23(’diffvgl’,[0,1],[2,8]) berekent nu voor het tijdsinterval [0, 1]de waarden van x1 (t) en x2 (t) met beginwaarden x1 (0) = 2 en x2 (0) = 8. De variabele x krijgt als waarde een array met twee kolommen, de eerste kolom bevat de waarden van x1 (t) en de tweede de waarden van x2 (t). Het plotten gebeurt door >> plot(t,x) Omdat x uit twee kolommen bestaat, worden er twee grafieken getekend. Indien men alleen het commando De opdracht >> ode23(’diffvgl’,[0,1],[2,8])
76
HOOFDSTUK 6. DIFFERENTIAALVERGELIJKINGEN
geeft worden ook de grafieken van x1 (t) en x2 (t) getekend omdat per default de ‘OutptFcn’ odeplot wordt gebruikt. Merk op dat men de functionfile bijvoorbeeld ook alsvolgt kan programmeren function xdot=diffvgl(t,x) A=[5,-2;7,-4]; x=[x(1);x(2)]; xdot=A*x;
2.3
Het richtingsveld
De oplossingen van de differentiaalvergelijking in voorbeeld 6.2 zijn krommen in de (t, x1 , x2 )ruimte gegeven door (t, x1 (t), x2 (t). Het richtingsveld wordt gegeven door de raakvectoren aan de oplossingskrommen, dus door de vectoren (1, x˙ 1 (t), x˙ 2 (t)). Merk op dat men het richtingsveld kan afleiden met behulp van de differentiaalvergelijkingen, de oplossingen hoeft men niet expliciet te weten. In het algemeen tekent men de oplossingen en het richtingsveld alleen in het (x1 , x2 )-vlak. Men tekent feitelijk de projectie op het (x1 , x2 )-vlak. Men noemt dit vlak ook wel het fasevlak of toestandsruimte (de toestand van het systeem op tijdstip t wordt gegeven door (x1 (t), x2 (t))). De geprojecteerde oplossing wordt een fase-pad of ook wel (deel van) een integraalkromme genoemd. In een gegeven punt (x1 (t0 ), x2 (t0 ) van de toestandsruimte kan men de rictingsvector uitrekenen m.b.v. de gegeven differentiaalvergelijkingen. In het voorbeeld 6.2 volgt voor het punt (x1 (t0 ), x2 (t0 ) = (1, 2) de richtingsvector (x˙ 1 (t0 ), x˙ 2 (t0 ) = (5 ∗ 1 − 2 ∗ 2, 7 ∗ 1 − 4 ∗ 2) = (1, −1). Het richtingsveld krijgt men door in ieder punt de bijbehorende richtingsvector te bepalen. Om een richtingsveld te kunnen tekenen moet men eerst in het (x1 , x2 )-vlak de punten aangeven waar men de richtingsvectoren wil tekenen. Hiervoor kan men het commando ‘meshgrid’ gebruiken. Het commando >> [X,Y]=meshgrid(-1:.1:1 ,-2:.2:2) maakt een matrix X waarvan iedere rij gelijk is aan de vector -1:.1:1 en een matrix Y waarvan iedere kolom gelijk is aan de vector -2:.2:2. Het elementsgewijs combineren van deze twee matrices geeft de roosterpunten (x, y), waarbij x de vector -1:.1:1 doorloopt en y de vector -2:.2:2 doorloopt, d.w.z. de punten (Xij , Yij ). De vectoren X en Y kunnen nu gebruikt worden om berekeningen in de roosterpunten uit te voeren. Bijvoorbeeld U = 5 ∗ X − 2 ∗ Y en V = 7 ∗ X − 4 ∗ Y berekent voor het voorbeeld 6.2 in ieder roosterpunt de componenten (U, V ) van de richtingsvectoren. Om een richtingsveld in een serie roosterpunten te tekenen kent MATLAB het commando ‘quiver’. Het commando >> quiver(X,Y,U,V)
3. HET SYMBOLISCH OPLOSSEN VAN DIFFERENTIAALVERGELIJKINGEN
77
tekent de vectoren met componenten (Uij , Vij ) in de punten (Xij , Yij ). De vectoren worden automatisch geschaald. Hierdoor wordt het richtingsveld soms wat onduidelijk. Bij de differentiaalvergelijking uit voorbeeld 6.2 worden de vectoren U en V verkregen door x˙ 1 en x˙ 2 te berekenen in de punten van het rooster, d.w.z. U = 5 ∗ X − 2 ∗ Y en V = 7 ∗ X − 4 ∗ Y . Oefening 6.5 Teken het richtingsveld voor de differentiaalvergelijking uit voorbeeld 6.2 op het rooster gegeven door de vectoren x = [−1 : .1 : 1] en y = [−1 : .1 : 1].
2.4
Plotten van integraalkrommen
De output van ode23 voor een stelsel van twee vergelijkingen als in voorbeeld 6.2 kan ook gebruikt worden voor het tekenen van fase-paden. Hiertoe moet men de oplossingen x1 (t) en x2 (t) plotten in het fasevlak, d.w.z. het (x1 , x2 )-vlak. Dit kan door na ode23 de kolommen van de matrix x tegen elkaar uit te zetten met >> plot(x(:,1),x(:,2)) Dit geeft een figuur met ´e´en fasepad. Wil men verschillende fasepaden met verschillende plot-opdrachten in dezelfde figuur laten tekenen dan moet men na de eerste plot-opdracht het commando ‘hold on’ geven. (zie §14). Tekent men fasepaden en richtingsveld in ´e´en figuur dan zullen de richtingsvectoren raken aan de fasepaden. Men kan ook een fasepad laten tekenen door de ode-solver als eigenschap ’OutputFcn’ mee te geven ’odephas2’. Het commando >> ode23(’diffvgl’,[0,1],[2,8],odeset(’OutpuFcn’,’odephas2’)) resulteert in een 2D fasepad vanuit de gegeven beginvoorwaarde.
3
Het symbolisch oplossen van differentiaalvergelijkingen
3.1
Eerste-orde differentiaalvergelijkingen
Voor het symbolisch oplossen van differentiaalvergelijkingen kent de symbolic toolbox van MATLAB het commando ’dsolve’. >> dsolve(’Dx1=a*x1’) ans = exp(a*t)*C1 men kan hierbij eventueel een beginvoorwaarde toevoegen.
78
HOOFDSTUK 6. DIFFERENTIAALVERGELIJKINGEN
>> dsolve(’Dx1=a*x1’,’x1(0)=2’) ans = 2*exp(a*t) Merk op dat in de meeste gevallen de oplossing van een differentiaalvergelijking niet op te schrijven is als een nette formulefunctie. In deze gevallen heeft symbolisch oplossen geen zin en is men aangewezen op numerieke technieken. Oefening 6.6 Ga na of de differentiaalvergelijkingen uit voorbeeld 6.1 en oefening 6.4 symbolisch op te lossen zijn. Merk op dat nieuwe functies worden ge¨ıntroduceerd. Ga met ‘help’ na wat dit voor functies zijn.
3.2
Stelsels eerste-orde differentiaalvergelijkingen
Er kunnen zonodig meerdere vergelijkingen en beginvoorwaarden worden ingevoerd. In dit geval moet men aangeven wat de output variabelen zijn. >> [x1,x2]=dsolve(’Dx1=5*x1 - 2*x2’,’Dx2=7*x1 - 4*x2’,’x1(0)=0’,’x2(0)=5’) x1 = 2*exp(-2*t)-2*exp(3*t)
x2 = -2*exp(3*t)+7*exp(-2*t) Indien men de functies x1 en x2 wil plotten kan dit met ’ezplot(x1)’ resp. ’ezplot(x2)’. Wil men de toestand (x1, x2) plotten dan dient men eerst de oplossingen x1 en x2 voor een aantal waarden van t te berekenen. Bijvoorbeeld: >> >> >> >>
t=[0:.1:1]; X=subs(x1,t); Y=subs(x2,t); plot(X,Y)
Oefening 6.7 Werk voorbeeld 6.2 uit. Plot in het fasevlak (in ´e´en figuur) de oplossingskrommen met beginwaarden (x1 (0), x2 (0)) resp. (2,8), (2,7), (2,6), (0.1,0.1), (-2,-8), (-2,-7), (-2,-6), (-0.1,-0.1) voor t ∈ [0, 1]. (Het verdient aanbeveling hierbij met ‘axis’ de assen vast te kiezen, bijv. met x1 en x2 tussen -10 en 10.) Doe dit zowel numeriek als symbolisch.
4. MATLAB EN DIFFERENTIAALVERGELIJKINGEN - OPDRACHTEN
3.3
79
Hogere-orde differentiaalvergelijkingen
Men heeft binnen matlab ook de mogelijkheid om hogere orde differentiaalvergelijkingen symbolisch op te lossen zonder deze eerst om te schrijven tot een stelsel eerste orde differentiaalvergelijkingen. Bijvoorbeeld de diffrentiaalvergelijking my 00 (t) + cy 0 (t) + ky(t) = 0 met beginvoorwaarden y(0) = 1, y 0 (0) = 0 kan men oplossen met het commando >> dsolve(’m*D2y+c*Dy+k*y=0’,’y(0)=1’,’Dy(0)=0’) Oefening 6.8 Beschouw de differentiaalvergelijking y 00 (t) + 2ky 0 (t) + 25y(t) = 0 , y(0) = 1 , y 0 (0) = 1 . a. Bepaal de oplossing van deze vergelijking voor k=0. De oplossingsfunctie beschrijft een trilling. Je kunt dit visualiseren door de oplossing te plotten. Is de trilling gedempt? Is de trilling periodiek? Wat is de frequentie van deze oplossing? b. Bepaal de oplossing van deze vergelijking voor k=1. Wat is de eigenfrequentie van het systeem? De oplossingsfunctie beschrijft een trilling. Je kunt dit visualiseren door de oplossing te plotten. Is de trilling gedempt? Is de trilling periodiek? Wat is de frequentie van deze oplossing? Plot de oplossing en de functie e−kt in dezelfde figuur. c. Bepaal de oplossing van de vergelijking y 00 (t) + 2y 0 (t) + 25y(t) = cos(5t) y(0) = 1 , y 0 (0) = 1 . De oplossingsfunctie beschrijft een trilling. Je kunt dit visualiseren door de oplossing te plotten. Is de trilling gedempt? Is de trilling periodiek? Wat is (op den duur) de frequentie en ampitude van deze oplossing?
4
MATLAB en differentiaalvergelijkingen - Opdrachten
We besluiten deze training met een aantal opdrachten waarin de in voorgaande paragrafen behandelde matlabcommando’s een rol spelen. U dient zelf de te volgen strategie te bepalen en te beslissen of u numeriek dan wel symbolisch wilt werken.
4.1
Opdrachten
Oefening 6.9 Beschouw de differentiaalvergelijking van een harmonische oscillator gegeven door een gedempt massa-veer systeem (vgl. [4, hfdst 3]). my 00 (t) + cy 0 (t) + ky(t) = 0
80
HOOFDSTUK 6. DIFFERENTIAALVERGELIJKINGEN
We gaan er van uit dat we een veer hebben die door een gewicht van 98 Newton 20 cm wordt uitgerekt. Dus k =490[N/m] en m =10[kg]. De demping c is zo gekozen dat c2 < 4mk, d.w.z. onderkritische demping. 10y 00 (t) + 50y 0 (t) + 490y(t) = 0 Deze differentiaalvergelijking is desgwenst te schrijven als een stelsel eerste orde differentiaalvergelijkingen. Stel x1 (t) = y(t) en x2 (t) = y 0 (t). Dan krijgen we de vergelijkingen: x˙ 1 = x2 x˙ 2 = −5x2 − 49x1 Plot de oplossing y(t) met beginwaarde y(t) = 0.2, y 0 (t) = 0 voor een voldoende groot tinterval en plot deze oplossing ook in het fasevlak. Wat gebeurt er bij andere beginwaarden? Plot ook voor bovenkritische demping, c > 140, de oplossing met dezelfde beginvoorwaarden. Kijk of je bij kritische demping, c = 140, beginvoorwaarden kunt vinden zodat de oplossing tenminste ´e´en nulpunt heeft (tenminste ´e´en doorgang door de evenwichtsstand). Oefening 6.10 Beschouw de Lotka-Volterra vergelijkingen voor een prooi-roofdier model (voor meer informatie zie [5]) dp dt dr dt
= (a1 − b1 r)p = (−a2 + b2 p)r
Met a1 , a2 , b1 , b2 > 0. Hierin is p(t) het aantal prooidieren en r(t) het aantal roofdieren. De prooidieren hebben een onuitputtelijke voedselbron en er is geen onderlinge concurrentie. Wanneer er geen roofdieren zijn ontwikkelt de populatie zich exponentieel. De prooidieren vormen voor de roofdieren de enige voedselbron. Wanneer er geen prooidieren zijn sterft de roofdierpopulatie uit. Men kan aantonen dat p(t) en r(t) aan de volgende relatie moeten voldoen a2 log p − b2 p + a1 log r − b1 r = constant Deze vergelijking geeft voor iedere constante een gesloten kromme in het (p, r)-vlak. De populaties ontwikkelen zich dus volgens een periodiek patroon. a. Neem de waarden a1 = 1.00, b1 = 0.10, a2 = 0.50 en b2 = 0.02. Plot de oplossingen p(t), r(t) en de bijbehorende kromme in het fasevlak voor enige zelfgekozen beginwaarden en een zelf gekozen tijdsinterval. Merk op dat het punt p = ab11 = 25, r = ab22 = 10 een stationair punt is in het fasevlak (constante oplossing). Opmerking: Vanwege numerieke oorzaken kan het zijn dat de getekende fasekromme niet een gesloten kromme is maar een spiraalkromme. Dit komt alleen naar voren bij plotten met grote t-intervallen. ‘ode45’ geeft hier het beste resultaat. b. De parameters a1 en a2 bepalen de ontwikkeling wanneer er geen interactie is. De parameters b1 en b2 bepalen de interactie tussen roofdier en prooi. Bekijk voor een aantal verschillende parametercombinaties hoe de populaties zich ontwikkelen. Merk op dat alleen het eerste kwadrant van het (p, r)-valk van belang is.
4. MATLAB EN DIFFERENTIAALVERGELIJKINGEN - OPDRACHTEN
81
Oefening 6.11 Beschouw een aangedreven veer-demper systeem gegeven door de vergelijking m¨ y + cy˙ + ky = cx˙ + kx Hierin is x de aandrijving of input en y de responsie of output. m=870 kg, k=70000 Nm−1 , c=5000 Nsm−1 (men kan denken aan de vering van een auto). We veronderstellen voor de input x(t) = sin(t). Bereken (kies zelf beginvoorwaarden, niet (0,0)) en teken de output (kies een interval waarop het periodiek gedrag zichtbaar is). Oefening 6.12 Beschouw als in de voorgaande opdracht een aangedreven veer-demper systeem gegeven door de vergelijking m¨ y + cy˙ + ky = k ∗ sign(sin(2πt)) m=870 kg, k=70000 Nm−1 , c=5000 Nsm−1 . Ga na wat de functie in het rechterlid precies voorstelt. Merk op dat de matlab-functie ’sign’ geen symbolisch argument accepteert. Bereken (kies zelf beginvoorwaarden, niet (0,0)) en teken de output (kies een interval waarop het periodiek gedrag zichtbaar is).
82
HOOFDSTUK 6. DIFFERENTIAALVERGELIJKINGEN
Bibliografie [1] The Student Edition of MATLAB, Version 4, User’s Guide / The Mathworks Inc., Prentice Hall, Englewood Cliffs, N.J., 1995. [2] Studiehandleiding Lineaire Algebra voor W (2Y650), diktaat nr. 2532, Technische Universiteit Eindhoven, 1997. [3] B. Kolman, Introductory Linear Algebra with Applications, 6th ed.,Prentice-Hall internatinal Inc., Upper Sadle River, 1997. [4] Syllabus Calculus voor W (2Y130), Technische Universiteit Eindhoven, 1998. [5] E.C. Pielou, Mathematical Ecology, John Wiley & Sons, New York, 1977. [6] E. P¨art-Enander et al, The MATLAB handbook, Addison-Wesley, Harlow, 1996.
83
