Cursus Matlab voor Technische Aardwetenschappen. K. Dekker I.A.M. Goddijn P. Sonneveld

Cursus Matlab voor Technische Aardwetenschappen K. Dekker I.A.M. Goddijn P. Sonneveld

1

2

Inhoudsopgave 1 Introductie in Matlab 1.1 Het invoeren van opdrachten 1.2 De help-faciliteit . . . . . . 1.3 Rekenen met getallen . . . . 1.4 Toekenningen . . . . . . . . 1.5 Werken met matrices . . . . 1.6 Matrix-bewerkingen . . . . . 1.7 Functies in Matlab . . . . . 1.8 Plotten . . . . . . . . . . . . 1.9 Printen . . . . . . . . . . . . 1.10 Thuis werken met Matlab .

. . . . . . . . . .

1 1 2 2 4 6 8 10 11 15 15

2 Programmeren in Matlab (m-files) 2.1 Scripts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Function files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16 16 18

3 Control flow (for, while, 3.1 De for-loop . . . . . . 3.2 De while-loop . . . . 3.3 Het if-statement . . .

if ...) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21 21 26 29

4 Driedimensionale plots 4.1 Ruimte-krommen . . 4.2 Oppervlakken . . . . 4.3 Staven en taarten . . 4.4 Animatie . . . . . . .

. . . .

32 32 33 36 36

. . . .

. . . .

. . . .

. . . . . . . . . .

. . . .

. . . . . . . . . .

. . . .

. . . . . . . . . .

. . . .

. . . . . . . . . .

. . . .

i

. . . . . . . . . .

. . . .

. . . . . . . . . .

. . . .

. . . . . . . . . .

. . . .

. . . . . . . . . .

. . . .

. . . . . . . . . .

. . . .

. . . . . . . . . .

. . . .

. . . . . . . . . .

. . . .

. . . . . . . . . .

. . . .

. . . . . . . . . .

. . . .

. . . . . . . . . .

. . . .

. . . . . . . . . .

. . . .

. . . . . . . . . .

. . . .

. . . . . . . . . .

. . . .

. . . . . . . . . .

. . . .

. . . . . . . . . .

. . . .

. . . . . . . . . .

. . . .

. . . . . . . . . .

. . . .

. . . . . . . . . .

. . . .

. . . . . . . . . .

. . . .

. . . . . . . . . .

. . . .

. . . .

ii

Voorwoord Deze cursus 1 is ontstaan naar aanleiding van geconstateerde problemen bij het gebruik van Matlab door ouderejaars studenten Toegepaste Aardwetenschappen. Het doel van deze cursus is enige basisvaardigheden op te doen in het gebruik van Matlab zodat het in de toekomst eenvoudiger wordt hier bij verschillende vakken gebruik van te maken. De oefenstof is verdeeld in een viertal hoofdstukken. In ieder hoofdstuk staan voorbeelden en opgaven. Het is de bedoeling dat de voorbeelden worden uitgeprobeerd en de opgaven gemaakt gedurende 7 donderdagochtenden, namelijk 9, 16 en 23 februari, 8, 15, 22 en 29 maart (8.45 - 10.30 uur) in zaal 1.97, Faculteit Civiele Techniek en Geowetenschappen, Stevinweg 1, 2628CN Delft. Er wordt een presentielijst bijgehouden van de aanwezige studenten. Wil je met een voldoende resultaat de cursus kunnen afsluiten en hiervan zo veel mogelijk profijt hebben dan is het belangrijk het een en ander rustig door te nemen en uit te proberen. Lezen alleen is echt niet genoeg. Er zijn een docent en twee studentassistenten aanwezig op bovenstaande tijden om, indien nodig, assistentie te verlenen. Bij hen moeten ook de gemaakte opgaven worden afgetekend. Delft, februari 2012

1

Er is hierbij dankbaar gebruik gemaakt van een cursus die jaren geleden is opgezet door A.T. Hensbergen

iii

iv

1

Introductie in Matlab

Start vanuit Windows Matlab op en sluit vervolgens alle schermen op het Command Window na. te klikken in de rechter bovenhoek van de te sluiten schermen of Dit kan door op door gebruik te maken van het menu. In het laatste geval kan er achtereenvolgens op de menu-items Desktop, Desktop Layout en Command Window Only worden geklikt. Matlab is een interactief systeem: de ingetikte commando’s (opdrachten), afgesloten door ‘Enter’, worden direct uitgevoerd, en indien gewenst verschijnen de resultaten direct op het scherm. Om ‘uitvoer’ op het scherm te onderdrukken moet een commando met een puntkomma worden afgesloten.

1.1

Het invoeren van opdrachten

Het intikken van opdrachten vraagt wel enige nauwkeurigheid. Het vergeten van bijvoorbeeld een * of een ) kan ertoe leiden dat de opdracht niet of verkeerd wordt uitgevoerd. Hieronder volgt een lijst met links de in de wiskunde gebruikelijke notatie en rechts de vorm waarin het ingetikt moet worden (a, b, x en y zijn getallen). Wiskundige notatie Intikken a+b a+b a−b a-b ab a*b 3xy 3*x*y a a/b b ab a^b √ x sqrt(x) π pi e exp(1) 3 − 4i 3-4*i sin x, arctan x, . . . sin(x), atan(x), ... x e , ln x exp(x), log(x) log x ( = 10 log x ) log10(x) |x| abs(x) ∞ Inf (i) Houd er rekening mee dat Matlab ‘case sensitive’ is. Dat wil zeggen dat er verschil is tussen hoofd- en kleine letters. Het getal π bijvoorbeeld wordt binnen Matlab geschreven als pi, het invoeren van Pi zal een foutmelding opleveren. (ii) Opdrachten in Matlab worden afgesloten door het indrukken van de Enter of RETURN-toets. Het antwoord verschijnt direct op het scherm. Door het commando af te sluiten met een ; (puntkomma) wordt de uitvoer op het scherm onderdrukt. 1

1.2

De help-faciliteit

Matlab beschikt over een uitgebreide help-faciliteit. Er zijn verschillende manieren om hiervan gebruik te maken. (i) Klik achtereenvolgens op de menu-items Help en Full Product Help. (ii) Klik op knop

.

(iii) Maak gebruik van het commando help onderwerp Met behulp van dit commando verschijnt er een toelichting bij onderwerp. Het typen van bijvoorbeeld: help elfun

(Druk nu de Enter of RETURN–toets in)

geeft informatie over de elementaire functies die Matlab kent. Lees (via ↑ ↓ of Page Up en Page Down) wat er staat. Met alleen het commando help verschijnt een lijst van mogelijke onderwerpen. Opmerking Het zal duidelijk zijn dat alleen het commando na de Matlab-prompt ( ) moet worden getypt. Hij wordt, in sommige voorbeelden, voor de duidelijkheid weergegeven.

1.3

Rekenen met getallen

Matlab is in feite een zeer uitgebreide programmeerbare rekenmachine. Om met het pakket kennis te maken en te wennen aan het invoeren van commando’s zullen we een aantal korte opdrachten uitvoeren. Opgaven 1. Voer in: 12/3


Op het scherm verschijnt de uitvoer van het commando in de volgende vorm: ans = 4 Voer ook eens in: 12\3


Wat valt er op ? Herhaal beide opdrachten nog eens met andere getallen. 2. Vraag eens met het commando help informatie op over SIN, Sin en sin en probeer vervolgens eens sin( π2 ) uit te rekenen. Wat valt op ? 2

Matlab is een numeriek pakket, dat wil zeggen dat er geen exacte berekeningen uitgevoerd kunnen worden. Het hangt van de opgegeven nauwkeurigheid af in hoeveel decimalen er wordt gerekend. Matlab rekent in ongeveer 16 decimalen nauwkeurig (waarvan er standaard vijf op het scherm verschijnen). Voer in: 12/9 Als antwoord verschijnt op het scherm:


ans = 1.3333 Rekenen met complexe getallen is voor Matlab geen probleem. Het getal i kent Matlab ook als i. het getal a + bi voer je in als a+b*i. Met de commando’s abs en angle kun je de modulus resp. het argument van een complex getal vinden, en verder zijn er de commando’s real, imag, conj en exp met voor de hand liggende acties (zo niet: raadpleeg help). Opgave

√ 3. Bereken (1 + 2i)(5 − 3i) en (1 + i 3)5 , √ √ 3 en ga met Matlab na dat e 4 πi = − 12 2 + 12 2i. Standaard geeft Matlab op het scherm 5–cijferige getallen met een vaste komma. Met behulp van het format–commando kun je opgeven welke notatie–vorm je wilt gebruiken. In onderstaande tabel staan verschillende opties: In de rechterkolom staan de representaties van π, van 1.23456 · 10−3 en van 1.23456 · 10−6 . Commando

Resultaat

Voorbeeld

format short

5 cijfers, ‘scaled fixed point’

3.1416 0.0012 1.2346e–006

format short e

5 cijfers, ‘floating point’

3.14163e+000 1.2346e–003 1.2346e–006

format long

15 cijfers, ‘scaled fixed point’

3.14159265358979 0.00123456000000 1.234560000000000e–006

3

format zonder toevoegingen is hetzelfde als format short. Voor overige mogelijkheden: raadpleeg help. Opgave 4. (a) Ga eens na hoeveel verschillende ‘formats’ Matlab kent. 77 (b) Geef weer in alle ‘formats’. 4 Ga nu weer over op de standaard-instelling .

1.4

Toekenningen

Als een matrix (of getal) meerdere keren gebruikt gaat worden kun je deze matrix toekennen aan een variabele. Je hebt hem dan als het ware een naam gegeven, en kunt hem via die naam weer oproepen. Als voorbeeld berekenen we de eerste vier machten van de matrix:   1 2 1  2 −1 2 . 1 3 −2 Voer in: A = [1, 2, 1; 2, -1, 2; 1, 3, -2] (i.p.v. komma’s kun je ook spaties nemen) Op het scherm verschijnt: A = 1 2 1

2 -1 3

1 2 -2

Voer in: A2 = A^2; A3 = A^3; A4 = A^4; Omdat de Matlab-commando’s zijn afgesloten met ; worden de resultaten onderdrukt. Vervolgens kun je de som hiervan berekenen door in te voeren: som = A + A2 + A3 + A4 (De spaties staan hier slechts voor de leesbaarheid. Matlab negeert ze.) De naam van een variabele hoeft natuurlijk niet beperkt te blijven tot één karakter. Het is aan te raden bij het maken van programma’s ‘logische’ variabele-namen te kiezen; dit zal de leesbaarheid en duidelijkheid zeer ten goede komen.

4

Opmerkingen (i) Als je niet zelf een variabele introduceert , kiest Matlab hiervoor de variable ans. Dit is hierboven al een paar keer gebeurd. Het is echter een goede gewoonte om de resultaten van elke berekening in een eigen variabele op te slaan, in dit geval bijv. in A2, A3, A4. (ii) Bij de eerstvolgende toekenning A = ... wordt de oude waarde van A weggegooid/vergeten. Een paar commando’s om de ‘actuele stand van zaken’ op te vragen of te be¨ınvloeden: Intikken who size(x) clear x clear all clc

de namen van alle variabelen die een waarde hebben worden opgevraagd; de afmetingen (aantal rijen, kolommen) van de variabele x worden opgevraagd; de variabele x wordt verwijderd; alle variabelen worden verwijderd; het scherm wordt leeg gemaakt (idem: Clear Command Window).

In Matlab is het alleen mogelijk aan een variabele een gedefinieerde waarde toe te kennen Het volgende geeft bijvoorbeeld problemen: clear x; functie = x^2 + 4 * sin(x) Om een functie in te voeren moet deze door de gebruiker eerst ‘geprogrammeerd’ worden. De gebruiker breidt dan de lijst Matlab-commando’s uit met een nieuw commando. In deze fase gaan wij hier niet verder op in. Het is niet mogelijk formulemanipulatie uit te voeren met Matlab, daartoe zijn andere pakketten beschikbaar (bijvoorbeeld: Maple, Mathematica). Het berekenen van de functiewaarde van boven bedoelde functie f in bijv. het punt π/3 kan als volgt gebeuren: x = pi/3 y = x^2 + 4*sin(x) Voer dit uit. Opmerking Met de ↑-toets kun je voorgaande commando’s terughalen en, al of niet gewijzigd, opnieuw uit laten voeren door op Enter te drukken. (Hiervoor hoeft de cursor niet aan het eind van de regel te staan.)

5

Met de Delete- en de Back-space-toets (dat is bij sommige toetsenborden de ← toets boven Enter) kun je symbolen weghalen; met de ← en → rechts van de Ctrl-toets kun je binnen een regel van links naar rechts lopen. Probeer dit uit!

1.5

Werken met matrices

Het basis-element van Matlab is de matrix. Speciale gevallen hiervan zijn een 1 × 1– matrix, d.w.z. een scalar, en een matrix bestaande uit slechts één rij of kolom, d.w.z. een vector. De grenzen van een matrix hoeven niet gedeclareerd te worden. (Ze kunnen opgevraagd worden met het commando size.) Bij het elementsgewijs invoeren van een matrix worden de verschillende rij-elementen gescheiden door spaties of door komma’s. De rijen worden gescheiden door puntkomma’s, of door naar een nieuwe regel te gaan. De getallen worden omsloten door twee rechte haken, [ en ].   1 2 3 De matrix A =  4 5 6  wordt als volgt ingevoerd: 7 8 9 A = [ 1, 2, 3 ; 4, 5, 6 ; 7, 8, 9] Op het scherm verschijnt: A = 1 4 7

2 5 8

3 6 9

Een (kolom)vector kan worden gezien als een matrix met één kolom kan kan dus als volgt worden ingevoerd : b = [ 1; 6; 7; 5; -3 ] Matrix-elementen kunnen afzonderlijk of in groepen opgeroepen of van een waarde voorzien worden, onderstaande tabel geeft een overzicht van een aantal mogelijkheden: Intikken A(i,j) A(:,j) A(i,:) A(k:l, m:n)

v(k:l)

het element in de i-de rij en j-de klom van de matrix A (aij ); de j-de kolom van de matrix A (aj ); de i-de rij van de matrix A; de ‘ondermatrix’ van de matrix A bestaand uit de k-de t/m de l-de rij en de m-de t/m de n-de kolom (aij met ( k ≤ i ≤ l en m ≤ j ≤ n)); de ‘deelvector’ van v met als elementen vi ( k ≤ i ≤ l).

6

En voor een vierkante matrix: diag(A) de hoofddiagonaal van A (als kolomvector); diag(A,k) de k-de nevendiagonaal van de matrix A. (Als k > 0 dan wordt de k-de nevendiagonaal boven de hoofddiagonaal genomen en als k < 0 dan de k-de nevendiagonaal onder de hoofddiagonaal.) Opgave 5. Voer de volgende 3 × 4-matrix in:   1 2 3 4 A =  5 6 7 8  1 1 1 1 Voorspel en controleer wat het resultaat is van de volgende commando’s: B = A(:, 3) C = A(1:3,2:4) D = diag(C,-1) A(1,1) = 9 A(2:3, 1:3) = [0 0 0; 0 0 0] A(1:3, 1) = [1, 1, 5] A(1:2, 1:2) = [-1, -1; 3, 3] A = A(3:3, 1:4) Opmerking Bij één van de voorgaande commando’s zou het voor de hand liggen als Matlab een foutmelding zou geven. Welke, en waarom ? Hoewel de laatste regel een voorbeeld is van ‘slordig programmeren’, kun je er wel duidelijk aan aflezen wat er in feite gebeurt bij een toekenning (Engels: assignment) van het type A = B. Eerst wordt de waarde van B uitgerekend, en vervolgens wordt deze waarde (een getal, een vector of een matrix) toegekend aan de variabele A. Als A reeds een waarde had dan wordt die ‘weggegooid/vergeten’. Het ‘=’teken betekent hier dus iets heel anders dan ‘is gelijk aan’. Je kunt A = B dan ook het best lezen als: ‘A wordt B’. In programmeertalen als Pascal, maar ook in computeralgebra pakketten als Maple en Derive, is hiervoor het symbool := gebruikelijk. Voorspel wat het resultaat zal zijn van de volgende commando’s: s = 1, s = s+1, s = s+5, s = s^2

7

1.6

Matrix-bewerkingen

Wat nu volgt is een greep uit het grote arsenaal aan matrix-bewerkingen waar de Matlabgebruiker standaard de beschikking over heeft. De commando’s met een worden daarna toegelicht. De commando’s met een hebben betrekking op begrippen die bij het vak Lineaire Algebra aan de orde zijn geweest. Voor alle bewerkingen geldt: is het je niet duidelijk wat er bedoeld wordt: raadpleeg help. C = A+B C = A-B C = A*B

de optelling van matrices; het verschil van matrices; de matrixvermenigvuldiging (Bij onjuiste afmetingen van de matrices A en/of B volgt een foutmelding.); C = A.*B de ‘elementsgewijze’ vermenigvuldiging (De matrices A en B moeten dezelfde afmetingen hebben.); C = A+k bij elk element van de matrix A wordt het getal k opgeteld; C = A^k de gewone machtsverheffing (k ∈ Z); bijv. te gebruiken voor berekening van A−1 ; C = A.^k de ‘elementsgewijze’ machtsverheffing (kan handig zijn bij het invoeren van een polynoom); C = A’ de getransponeerde AT ; C = A\B berekent voor een niet-singuliere, vierkante matrix A de unieke oplossing van de vergelijking A X = B en voor een overbepaald stelsel (meer vergelijkingen dan onbekenden) een kleinste kwadraten oplossing; C= B/A oplossing van X A = B, analoog aan vorige commando; [L,U,P] = lu(A) berekent de L U -ontbinding van P A; L = eig(A) L is een kolomvector met de (mogelijk complexe) eigenwaarden van de matrix A; [X,D] = eig(A) D is een diagonaalmatrix met op de diagonaal de eigenwaarden van de matrix A , de kolommen van X zijn bijbehorende eigenvectoren; [Q,R] = qr(A) produceert een Q R-ontbinding van matrix A; A = diag(v) geeft diagonaalmatrix A met diagonaal v(1), v(2), ..., v(n); A = zeros(n) A wordt de n×n nulmatrix; A = zeros(n,m) A wordt de n×m nulmatrix; A = eye(n) A wordt de n×n identiteitsmatrix; A = ones(n,m) A wordt de n×m matrix met louter 1-en; x = a:b x wordt de vector (a, a + 1, a + 2, . . . , a + n h) a ≤ b, waarbij n = bb − ac; x = a:h:b x wordt de vector (a, a + h, a + 2h, . . . , a + n h) a ≤ b, b−a waarbij n = b c; h

8

Opmerking De notatie bxc wordt gebruikt voor de ‘entier’ van x d.w.z. ‘het grootste gehele getal kleiner of gelijk aan x’. Een paar voorbeelden ter illustratie (verklaar waar nodig de foutmeldingen): (i) Voer in: x y z w

= = = =

[1,3,5,8] [2,2,0,6] x.*y x*y’ (wat gaat ‘voor’ ?

* of ’ ? )

(ii) Voer in: s = 1:10 s1 = s^2 s2 = s.^2 s3 = s.^3 (iii) Zie Lay [1, p. 381, Theorem 3]     1 1 2    Laat A = 1 0 en b = 1  0 1 5 Heeft de vergelijking A x = b oplossingen ? Bereken u = A\b, en onderzoek of Au = b. Merk op dat dit in overeenstemming is met de uitleg bij het commando u = A\b. In welke deelruimte van R3 ligt de vector Au ? En de vector b − Au ? Om het ‘denken in rijen en kolommen’ te bevorderen: Opgaven 6. (a) Voer, door gebruik te maken van het commando a:h:b de volgende matrix in: x = 1 3 5 · · · 15 (b) Zoek uit hoe het commando linspace(a,b,n) gebruikt kan worden om dezelfde matrix in te voeren en doe dit. (c) Wat is de relatie tussen h en n? (d) Voer, op twee verschillende manieren, de matrix y = x0 x1 x2 · · · xn in zodat x0 = 2, xn = 22, n = 10 en het verschil van twee opeenvolgende elementen steeds gelijk is.

9

7. Voer (op een  1 3  1 3   1 3   1 3  . . A =   .. ..  . .  .. ..   . .  .. .. 1 3       B =      

1 1 1 .. . .. . .. .

handige manier) de volgende twee matrices in:  5 · · · 15 5 · · · 15   5 · · · 15   5 · · · 15   .. .  . · · · ..  (A heeft 50 rijen. Gebruik de .. .  commando’s a:h:b en ones. ) . · · · ..   .. ..  . ··· .  5 · · · 15

1 2 4 .. . .. . .. .

1 2

50

1 3 9 .. . .. . .. . ···

1 4 16 ··· ··· ··· ···

··· ··· ··· .. . .. . .. .

1 9 92 .. . .. . .. .

      .     

· · · 950

Aanwijzing: maak eerst een matrix met 51 rijen van de vorm 1 2 3 4 · · · vervolgens een matrix waarvan de 9 kolommen de getallen 0 1 2 3 · · · 50 bevatten en tenslotte de matrix B.

9 en

Overigens zullen we later een andere manier zien om een ‘gestructureerde’ matrix in te voeren, namelijk met een for-loop. 8. Gegeven is de rijvector v = (1, 2, · · · , 100). Bereken de som van de elementen van v op twee verschillende manieren. (a) Door het produkt te nemen met een geschikte kolomvector en (b) door een geschikt Matlab-commando te gebruiken. (c) Bereken analoog de som van de eerste honderd kwadraten, dat wil zeggen: 1 + 4 + 9 + . . . + 992 + 1002 . (Antwoord: 338350)

1.7

Functies in Matlab

Binnen Matlab kunnen een aantal elementaire functies worden gebruikt. Als je wilt weten welke: vraag nog maar eens op met het Matlab-commando help elfun

10

Als je invoert: x = 0:0.1:pi; (Zie lijst met matrixbewerkingen.) y = sin(x) dan worden er twee variabelen x en y ‘aangemaakt’ met waarden (0, 0.1, 0.2, ... , 3.0, 3.1) resp. (sin(0), sin(0.1), sin(0.2), ... , sin(3.0), sin(3.1) ). Matlab neemt dus elementsgewijs de sinus van x, en slaat de ‘functie’ als het ware op in een tabel. Meer ‘weet’ Matlab niet van de ‘functie’ y. Het is zinloos om bijvoorbeeld met y(0.25) de sinus van 0.25 te willen berekenen. (Kijk maar wat er gebeurt.) Wat stelt y(5) voor ? De ‘elementsgewijze’ matrix-operaties .* en .^ zijn met name handig om ‘functies’ in te voeren: Voorbeeld We berekenen de functiewaarden f (x) = 13 x2 + x sin x voor x in : { 0, 0.1, 0.2, . . . , 2.8, 2.9, 3.0}. Eerst voeren we de x-waarden in : x = 0:0.1:3.0;

(Dit geeft een rijvector met 31 elementen.)

Dit kan ook als volgt : x = linspace(0,3.0,31) Vervolgens berekenen we f (x(k)) voor k = 1, 2, . . . , 30, 31. Met behulp van de operatoren .* en .^ lukt dit als volgt: y = (1/3)*x.^2 + x.*sin(x) Het resultaat van de toekenning y = (1/3)*x^2+x.*sin(x) is niet direct leesbaar te noemen. We krijgen een beter idee wanneer we x en y transponeren of grafisch weergeven.

1.8

Plotten

Met het plot–commando kun je vectoren (en matrices) grafisch laten weergeven. Als je een ‘functie’ invoert door een rij x-waarden met bijbehorende y-waarden te geven, krijg je met het commando plot(x,y) min of meer de grafiek op het scherm: Intikken plot(y)

tekent de punten (1, y(1)), (2, y(2)), . . . , (n, y(n)) en verbindt ze met rechte lijnen;

plot(x,y) doet hetzelfde voor de punten (x(1), y(1)), (2, y(2)) . . . , (x(n), y(n)) x en y mogen hier rij- of kolomvectoren met een gelijk aantal elementen zijn; figure er wordt een nieuw Plot Window geopend. 11

Voorbeelden Voer in:

x = [2,5,4,8,10,4] y = [6,3,9,1,4,0] plot(x,y)

Voer in:

figure; x = 0:0.1:3.2; y = sin(x); plot(x,y); (Geen paniek, zie eerstvolgende opmerking! )

Voer in:

figure; x = linspace(0,3.2,33); y = sin(x); plot(x,y);

Opmerking Bij een volgende plot-opdracht komt het plaatje niet altijd direct in beeld. Als een Plot Window zich achter het Command Window bevindt, kun je het op verschillende manieren naar voren halen: (i) Als een deel van het Plot Window wel zichtbaar is, klik je het simpelweg aan met de muis. (Je kunt een deel zichtbaar maken door het Matlab Window iets te verkleinen.) (ii) Door in het menu op Windows te klikken en vervolgens op A Figure 1, B Figure 2,· · · . (iii) Door te ‘bladeren’ via Alt-Tab. Nog een voorbeeld: Voer in:

x = 1:10; y = 20:30; plot(x,y); Wat gaat hier fout? Bij het plotten van meerdere grafieken in één plaatje is het soms handig om verschillende kleuren en tekens te gebruiken. De volgende tabel biedt voldoende variatiemogelijkheden; voor het volledige aanbod: raadpleeg help plot. Symbool Kleur y geel (yellow) r rood g groen b blauw w wit

Symbool . o -. *

12

Stijl ‘stip’ ‘rondje’ ‘lijn’ ‘streep-stip’ ‘ster’

Voorbeeld Voer in:

x = 0:0.1:3.2; y = sin(x); z = x.*sin(x); (Elementsgewijze matrixvermenigvuldiging!) plot(x,y,’y-’, x,z,’r-.’);

Er zijn twee manieren om meerdere grafieken in één plaatje te krijgen: met één plotopdracht, zoals boven, of met een aantal plot-opdrachten na elkaar, met tussendoor het commando hold on. Voorbeeld Voer in: x = 0:0.1:3.2; y = sin(x); x2 = 0:0.3:3; y2 = sin(x2); (Regel terughalen en wijzigen!) plot(x,y,’g’); Ga nu naar het Command Window terug. Het Plot Window is dan nog steeds aanwezig. Voer vervolgens in:

hold on; plot(x2,y2,’yo’);

en de twee grafieken staan in één plaatje (achter het Command Window ). Alle volgende plot-opdrachten geven grafieken in hetzelfde plaatje. Wil je een heel nieuw plaatje maken dan kun je dat doen door het commando hold off te gebruiken, of door het Plot Window te sluiten. Met onderstaande commando’s kun je de ‘layout’ van een bestaande plot be¨ınvloeden: Intikken axis([xmin xmax ymin ymax])

de ‘x-coördinaat’ loopt van xmin tot xmax, de ‘y-coördinaat’ van ymin tot ymax;

axis(’off’)

de assen worden niet getekend;

title(’tekst’)

geeft de titel tekst (tekst moet tussen quotes gezet worden.);

xlabel(’tekst’)

er wordt tekst langs de x-as gezet;

ylabel(’tekst’)

analoog aan xlabel;

gtext(’tekst’)

doet het (evt. nog lege) Plot Window in beeld verschijnen, en zet tekst op de positie die je met de muis (linkerknop) aanklikt.

Ook kun je de ‘layout’ van een bestaande plot be¨ınvloeden door op het item Insert in het menu van het Plot Window te klikken en vervolgens op X Label, Y Label, Title etc.

13

Voorbeeld Voer in:

p = [2, 5, 5, 2, 2]; q = [2, 2, 5, 5, 2]; plot(p,q); title(’vierkant ?’)

Het resultaat is het ‘vierkant’ met hoekpunten (2, 2), (5, 2), (5, 5) en (2, 5). Door op de x-as en y-as het interval [0,7] te nemen krijg je een duidelijker beeld: axis([0,7,0,7]); Wil je tenslotte dat het vierkant er als een vierkant uitziet, geef dan de opdracht: axis(’square’); Probeer tenslotte het volgende uit: gtext(’Dit is het punt (5, 5)’); Opmerking De zelf gedefinieerde instellingen m.b.v. axis blijven geldig zolang je met hetzelfde plaatje werkt. Opgaven 9. Teken de (gladde) grafiek van f (x) = sin x op het interval [0, 4π]. Verdeel hierbij [0, 4π] in 100 deelintervallen. kπ Geef daarna in de grafiek de punten kπ , sin , k = 0, 1, . . . , 24 met een * aan. 6 6 10. De kromme gedefinieerd door: x = x(t), y = y(t), a ≤ t ≤ b kun je met Matlab eenvoudig op het scherm krijgen. We nemen eerst: x = 3 cos(t), y = sin(t), t ∈ [0, 2π]. 2π 2π , 2 100 , . . . . en vervolgens rijvectoren Voer eerst een rijvector t in met als elementen: 0, 100 x en y met de bijbehorende x en y-coördinaten. Plot vervolgens de kromme met bovenstaande parametrisering. 11. Maak mooie plaatjes van een aantal Lissajoux-figuren. Bepaal zelf een geschikt interval voor t en een geschikte interval-verdeling (zodat de volledige kromme glad op het scherm verschijnt). (a) x(t) = cos(2t), y(t) = sin(5t) (b) x(t) = cos(2t + π/3), y(t) = sin(5t) (c) x(t) = cos(5t + π/4), y(t) = sin(4t) 12. Teken de kromme z = t eit , 0 ≤ t ≤ 4π. Aanwijzing: je kunt t en z vrijwel letterlijk invoeren. z wordt dan een vector van complexe getallen. Geef ook eens de stijl-optie ‘punt’ of ‘ster’ mee; dan kun je zien hoe de ‘snelheid’ verandert.

14

13. Zie Lay [1, §6.6, Exercise 12] Tussen de systolische bloeddruk p (in milimeters kwik) en het gewicht w (in Engelse ponden) van een gezond kind bestaat , bij benadering, de volgende relatie: β0 + β1 ln(w) = p In de volgende tabel zijn experimentele data opgenomen : w 44 p 91

61 98

81 113 131 103 110 112

(a) Bereken, door gebruik te maken van de ‘kleinste kwadratenmethode’, een benadering van β0 en β1 . (b) Teken in één plaatje p als functie van w ( 40 ≤ w ≤ 135 ) en de punten uit de tabel. Geef de punten de stijl-optie ‘ster’ mee. (c) Geef een schatting van de systolische bloeddruk van een kind met een gewicht van 100 pond. 14. Zie Lay [1, §6.6, Exercise 13] Om inzicht te krijgen in de wijze waarop een vliegtuig start, is de horizontale positie van een vliegtuig gemeten en wel elke seconde tussen t = 0 en t = 12. De gemeten posities (in feet) zijn: 0, 8.8, 29.9, 62.0, 104.7, 159.1, 222.0, 294.5, 380.4, 471.1, 571.7, 686.8, 809.2. (a) Bepaal de ‘kleinste kwadratenkromme’ met als vergelijking: y = β0 + β1 t + β2 t2 + β3 t3 . (b) Teken in één plaatje de onder (a) genoemde kromme en de metingen waarop deze kromme is gebaseerd. (c) Geef een benadering van de snelheid van het vliegtuig op t = 4.5 seconde .

1.9

Printen

Als je het plaatje op het scherm er acceptabel uit vindt zien, dan kun je het als volgt vastleggen op papier (= printen): klik in het menu van het Plot Window achtereenvolgens op de items File en Print... en tenslotte in het Print Window wat dan zichtbaar wordt op de OK–knop. Als de opdracht is verwerkt kom je weer terug in het Plot Window. (Printen vanuit Matlab geeft nogal eens complicaties. Voordat je bij weigering om te printen hulp inroept van een begeleider of systeembeheerder, kijk dan eerst of het drukken op de ‘Start’–knop van de printer soelaas biedt.)

1.10

Thuis werken met Matlab

Omdat de TUDelft een campuscontract heeft afgesloten met de firma ‘The Mathworks’ kan een installatiepakket van Matlab 2008b worden gedownload via blackboard. Daar staan ook systeemvereisten en installatie-instructies vermeld. 15

2

Programmeren in Matlab (m-files)

In de eerste instantie heb je kennis gemaakt met een aantal basiscommando’s van het pakket Matlab, alsmede met enkele van de plot-mogelijkheden. Je gaat nu , een klein beetje, leren programmeren in Matlab. Laten we om te beginnen eens kijken naar de structuur van Matlab. Matlab bestaat uit een kleine kern, met daarin basiscommando’s als +,

.*,

:,

rand, who, size, plot, det, sin, exp

en daaromheen gebouwd een scala aan al dan niet specialistische procedures en functies, opgeslagen in zogenoemde .m-files. Voorbeelden hiervan zijn: sum, mean, hist, quad, fzero, trace De (implementaties van de) basis-commando’s zijn voor de gebruiker niet zichtbaar en niet te be¨ınvloeden. De gebruiker kan zelf m-files (of M-files) toevoegen. Die optie maakt dragelijk dat het Command Window nogal beperkt is in z’n mogelijkheden. Commando’s kunnen niet (opnieuw) worden uitgevoerd door hier op te klikken en het telkens terughalen van eerder ingetypte commando’s gaat bij een wat ingewikkelder probleem al snel irriteren.

2.1

Scripts

Het is mogelijk een aantal commando’s in een zogenaamde script zetten (ofwel een programma schrijven) en die vervolgens te ‘runnen’. Ook kunnen er procedures en functies, zogenaamde function files, worden gemaakt, die vanuit het Command Window of vanuit een script kunnen worden aangeroepen. Voor het maken van een script file kun je op de volgende twee manieren te werk gaan: (i) Klik achtereenvolgens op de menu-items File, New en M-File of klik op de knop . De M-file Editor/Debugger opent een nieuw document, waarin een programma kan worden geschreven. Ben je klaar met schrijven dan moet je het ‘programma’ (dat nu alleen nog maar op het scherm staat) opslaan in een daarvoor geschikte map. Klik de menu-items File en Save As aan en maak eventueel eerst een map aan waarin je het script wilt opslaan. Geef vervolgens de file een naam, bijvoorbeeld: naam.m en bewaar haar. De extensie .m is belangrijk! (Als jij hem niet geeft, dan voegt de editor die standaard toe.) (ii) Maak met een simpele ‘ASCII-editor’, bijvoorbeeld Notepad een file aan met Matlab-commando’s. Sla deze file file op dezelfde manier op als onder (i) staat beschreven. 16

Voor het ‘runnen’ van het script kun je nu op de volgende twee manieren te werk gaan: (i) Voer in het Command Window in: naam (ii) Klik in het menu van de M-file Editor/Debugger op Debug en Run of klik op de knop . Voorbeelden (Zie ook § 1.8 ) Maak een .m-file waarin de volgende commando’s staan : % Plotten van de sinus op het interval [0,4pi] x = 0:0.1:4*pi; y = sin(x); plot(x,y); hold on; (toegelicht in § 1.8: ‘plotten’) x = 0:pi/6:4*pi; y = sin(x); plot(x,y,’r*’); Bewaar dit onder de naam plotsin.m (of iets anders, mits nieuw en met extensie .m), ga naar het Command Window en voer in: plotsin Het voordeel moge duidelijk zijn: mocht je achteraf besluiten dat je toch liever groene rondjes in plaats van rode sterretjes hebt, dan wijzig je de laatste regel van de .m-file. Voer iets dergelijks uit. Maak een .m file waarin de volgende commando’s staan : % Het berekening van het kwadraat van de natuurlijke % getallen tussen 90 en 100 x = 90:100; y = x.^2; disp([x;y]’); Bewaar dit bijvoorbeeld onder de naam kwadraat.m en klik op de knop

.

Opmerkingen (i) Het is een goede gewoonte om een programma te beginnen met één of enkele regels commentaar over de werking van het programma (Een commentaarregel wordt voorafgegaan door %.). Door in te voeren: help naam 17

krijg je dit commentaar later snel in beeld. Commentaarregels na een blanco regel of na een Matlab-opdracht worden niet getoond. Zij worden opgevat als intern commentaar. (ii) Krijg je na een poging naam uit te voeren de melding: ??? Undefined function or variable ’naam’. dan heb je een tikfout gemaakt of Matlab kent het pad naar de file naam.m niet. Door de ‘Current Directory’ te wijzigen in het pad naar de file naam.m kun je het te klikken en vervolgens laatste probleem opheffen. Dit kan door op de knop het juiste pad te kiezen. (iii) Zoals de naam al aangeeft kun je de M-file Editor/Debugger ook gebruiken om fouten in je programma’s op te sporen en te verwijderen ( = ‘debuggen’) Kijk eens onder de menu-items Debug en Breakpoints en let eens op de kleuren die de editor gebruikt. Gebruik je een gewone ‘ASCII-editor’ voor het maken van m-files dan mis je natuurlijk deze extra’s. (iii) Laat de naam van een .m file nooit met cijfers beginnen! Matlab herkent scripts met dergelijke namen niet. (iv) In een script moeten alle Matlab commando’s worden afgesloten met ; (puntkomma). Resultaten kunnen zichtbaar worden gemaakt door het commando disp te gebruiken.

2.2

Function files

Zoals gezegd, er kunnen twee soorten .m-files onderscheiden worden: scripts en function files. In scripts worden simpelweg een aantal commando’s in een .m-file gegroepeerd. Door middel van function files kun je zelf functies en procedures aan Matlab toevoegen. Scripts kun je ‘runnen’, function files roep je aan. De eerste regel van een function file luidt in het algemeen als volgt: function y = functienaam(x) Ook hier is het een goede gewoonte om te te beginnen met een of meer regels commentaar over de werking van de functie.

18

Met help functienaam kun je dit commentaar later eenvoudig opvragen. Voorbeelden (i) De functie f : R → R met functievoorschrift f (x) = x2 e−x kun je toevoegen als m-file newfunc:

2

function y = newfunc(x) % newfunc(x) berekent (elementsgewijs) x^2 * exp(-x^2) y = x.^2 .* exp(-x.^2); Doe dit, en maak een plot van deze functie op het interval [0, 5]. (ii) De volgende functie berekent het grootste verschil tussen twee opeenvolgende elementen van een vector: function result = maxdiff(x) % maxdiff(x) berekent het grootste verschil |x(i) - x(i-1)| % van vector x L = length(x); v = x(2:L) - x(1:L-1); result = max(abs(v)); Als je deze m-file opslaat als maxdiff.m, dan kan bijvoorbeeld het ‘maximale verschil’ tussen twee opeenvolgende elementen van de vector y = (f (0), f (0.1), f (0.2), · · · , f (3.0)) worden berekend door in te voeren: x = 0:0.1:3; y = newfunc(x); z = maxdiff(y) of door het volgende script te schrijven en te ‘runnen’: x = 0:0.1:3.0; y = newfunc(x); z = maxdiff(y); disp(z); Opmerkingen (i) De variabelen result en x in de regel function result = functienaam(x) zijn, hoe kan het anders in Matlab, matrices. x heet de invoervariabele, result de uitvoervariabele. In plaats van een enkele variabele, kun je ook géén of meerdere variabelen als invoer en/of uitvoer nemen. Enkele voorbeelden van een mogelijke aanheffen zijn:

19

function result = functienaam (geen getalsmatige invoer) function [res1, res2, ..., resk] = functienaam (geen getalsmatige invoer) function functienaam(x) (geen getalsmatige uitvoer) function functienaam(x1, x2, ..., xn) (geen getalsmatige uitvoer) function [res1, res2, ..., resk] = functienaam(x1, x2,..., xn) (meerdere in-en uitvoervariabelen) (ii) Alle Matlab-opdrachten in function files behalve de aanhef moeten worden afgesloten met een ; (punt-komma). Doe je dit niet kun je ongewenste uitvoer krijgen. Kijk ook nog eens terug naar de voorbeelden van function files. (iii) Een belangrijk verschil tussen scripts en function files is het karakter van de variabelen. De variabelen die in een function file gebruikt worden zijn lokaal. Ze leven als het ware alleen binnen de procedure. We gaan hier niet verder in op dit punt. Tot slot: (iv) Voor één-regelige functies als in voorbeeld (i) kent Matlab (vanaf versie 5) het commando inline. De functie newfunc is daarmee in te voeren als newfunc = inline(’x.^2 .*exp(-x.^2)’, ’x’); waarna je bijvoorbeeld newf unc(5) kunt berekenen via newfunc(5) Vraag desgewenst hulp op over inline. Opgaven 15. Maak een function file verwissel.m die als invoer een matrix heeft en twee indices i, j en als uitvoer een matrix die gelijk is aan de ingevoerde matrix op de i-de en j-de rij na; de i-de en j-de rij zijn verwisseld. Maak hierbij gebruik van zo min mogelijk hulpvariabelen. Als je handig bent kan het geheel zonder. 16. Maak een function file inwprod.m die als invoer twee kolomvectoren heeft en als uitvoer het inwendig product van deze vectoren. 17. Maak een function file stat.m die als invoer een vector x heeft, die bijvoorbeeld tentamenresultaten bevat, en als uitvoer het gemiddelde en de standaarddeviatie van deze resultaten. 18. Maak een function file plaatje.m die als invoer a, b ∈ R en n ∈ N\{0} heeft en als uitvoer een plaatje van de grafiek van de functie f : [α, β] → R gegeven door f (x) = xn . Hierbij α = min(a, b) en β = max(a, b). 20

3

Control flow (for, while, if ...)

3.1

De for-loop

Met de commando’s: for k = 1:n ; ; ... ... ; end; wordt de serie opdrachten ,, ... , n keer uitgevoerd. Het verdient aanbeveling om elke opdracht af te sluiten met een puntkomma, om onnodige output op het scherm te onderdrukken. In het algemeen zullen de uit te voeren opdrachten afhangen van de zogenoemde loopvariabele k. Het is een goede gewoonte om de loop-opdrachten in te laten springen, en niet meer dan één opdracht per regel te tikken. Voorbeelden (i) Methode van Newton Zie: Stewart [2, §4.8] √ 2 is een nulpunt van de functie f : R → R gegeven door f (x) = x2 − 2. √ Als een benadering x0 van 2 bekend is kan een, hopenlijk betere, benadering x1 worden gevonden door de raaklijn aan de grafiek van f in (x0 , f (x0 )) te snijden met de x-as. Deze raaklijn heeft als vergelijking y − f (x0 ) = f 0 (x0 )(x − x0 ) en dus x1 = x0 −

f (x0 ) ofwel f 0 (x0 )

(x20 − 2) 2 x0 Met de volgende commando’s kun je dus de wortel uit 2 benaderen (maak een m-file squareroot.m) : √ x = 1; (De eerste benadering, x0 , van 2.) for k = 1:50 x = x - (x^2 - 2)/(2*x); end; disp(x);

x1 = x0 −

21

De opdracht die 50 keer herhaald wordt is: bereken x −

(x2 − 2) en geef de 2x

variabele x deze waarde. Voor wie nog nooit geprogrammeerd heeft: kijk eens wat er gebeurt als in de derde regel de puntkomma weggelaten wordt. (ii) De rij van Fibonacci 0 1, 1, 2, 3, 5, 8, 13, . . . . . is vastgelegd door de eerste twee waarden en het voorschrift fk+1 = fk + fk−1 Het volgende ‘programma’ genereert de eerste vijftig Fibonacci-getallen. format short e; f = zeros(50,1); (Preallocatie !) f(2) = 1; (Merk op dat f (1) = 0.) for k = 3:50 f(k) = f(k-1) + f(k-2); end; disp(f); format; Na het uitvoeren van de for-loop is f een vector met 50 kentallen, en het k-de element van f is precies het k-de Fibonacci-getal. Als je alleen ge¨ınteresseerd bent in de waarde f50 , dan spaar je op de volgende manier enige geheugenruimte uit (maak een m-file fibo.m): format short e; f1 = 0; f2 = 1; for k = 1:48 f3 = f1+f2; f1 = f2; f2 = f3; end; disp(f2); Er zijn nu slechts drie variabelen f1, f2 en f3 in het spel; na elke stap wordt de waarde die niet meer nodig is ‘weggegooid’. Ook weer voor wie nog weinig of geen programmeerervaring heeft: na de volgende wijzigingen zijn de eerste acht stappen goed te volgen. Laat k lopen van 1 t/m 8, en voeg de volgende commando’s als zesde regel (direct boven end) toe: disp([f1, f2, f3]), pause(2); (iii) Het is ook mogelijk om binnen een for-loop een andere for-loop te plaatsen. We spreken dan van een geneste loop.

22

clear all;clc; A = zeros(8,8); (Preallocatie !) for k = 1:8 for l = 1:8 A(k,l) = k+l-1; end; end; (Twee keer end !) disp(A); (Voor het resultaat !) Hiermee wordt de 8×8 matrix A met Akl = k + l − 1 geconstrueerd. (Nu wordt het wel echt vervelend de puntkomma in de loop weg te laten.) Opmerking Vermijd for-loops als er een (efficiënter) alternatief voorhanden is. Onder andere door handig gebruik van het Matlab commando : is dit vaker mogelijk dan je (vooral de wat meer ervaren programmeur !) zou verwachten, en het kan de rekentijd aanzienlijk bekorten. Voorbeeld De som van de eerste N kwadraten kun je berekenen met de commando’s (Maak een function-file somkwad.m waarbij N in- en het resultaat uitgevoerd kan worden.): s = 0; for k = 1:N s = s + k^2; end; maar ook met de commando’s (Maak opnieuw een function file waarbij N in- en het resultaat uitgevoerd kan worden.): v = 1:N; enen = ones(N,1); kwad = v.^2; somkwad = kwad*enen of zo je wilt (maar leesbaarder wordt het er niet op) als: somkwad = ((1:N).^2)*ones(N,1) Met de commando’s tic en toc kun je de verstreken tijd in beeld krijgen. Voeg in de twee bovengenoemde function files als eerste regel toe: tic en als laatste regel: toc en roep ze aan voor N = 10, 100, 1000 en 10000.

23

Opgaven  1  1   1   1  . 19. A =   ..  .  ..   .  .. 1

3 3 3 3 .. . .. . .. .

5 5 5 5 .. . .. . .. .

··· ··· ··· ··· ··· ···

··· 3 5 ···

15 15 15 15 .. . .. . .. .

        (A heeft 50 rijen.)      

15

Voer matrix A op drie verschillende manieren in. (a) Eerst door twee geneste loops te gebruiken en daarna (b) door één loop te gebruiken en (c) door geen enkele loop te gebruiken. (Zie ook 7.)   1 1 1 1 ··· 1  1 2 3 4 ··· 9    2   1 4 9 16 · · · 9   .. .. ..   .. .. . . . ··· . . . 20. B =   . . .. .. ..   .. .. . · · · . .     . .  . . . . . . . .  . . . ··· . .  1 250 · · · · · · · · · 950 Beantwoord dezelfde vragen als voor de matrix A uit voorgaande opgave. 21. n! (spreek uit ‘n faculteit’) is gedefinieerd door: n! = n(n − 1)(n − 2) · · · · 3 · 2 · 1 (en apart: 0! = 1 ). 20 X 1 1 : S= Bereken de som van de eerste eenentwintig getallen n! n! n=0 Waarvan is deze som een benadering ? 22. Construeer de 30 bij 30 matrix A met op de hoofddiagonaal 2-en, direct boven en onder de diagonaal 1-en, en verder overal 0-en. Doe dit op twee verschillende manieren: (a) Met behulp van het commando diag(.., ..): als v een (rij- of kolom-)vector is van lengte m, dan wordt met V = diag(v, k) een vierkante matrix V geconstrueerd met de vector v als k-de nevendiagonaal, en met verder overal 0-en.

24

(b) Via enkele loops. (Begin met het invoeren van een 30 × 30 nulmatrix.) Bekijk ter controle de 5 × 5 rechter-onder-matrix, en bereken de determinant van A. (Antwoord: 31) 23. (a) Maak een function file met als invoer m ∈ N en een rijvector x = (x1 , x2 , · · · , xn ) en als uitvoer de (m + 1) × n matrix   1 1 1 ··· 1  x1 x2 x3 · · · xn   2   x1 x22 x23 · · · x2n    .. .. ..   .. . . . ··· . . B =   . . .. .   .. .. . · · · ..     .  . . . .. .. · · · ..   .. xm xm · · · · · · xm 1 2 n Aanwijzing: maak eerst een matrix met m + 1 rijen van de vorm x1 x2 x3 · · · xn en vervolgens een matrix waarvan de n kolommen de getallen 0 1 2 · · · m bevatten en tenslotte de matrix B. Kijk hiervoor ook eens terug naar de opgaven 7 en 20. (b) Laat p(x) = a0 + a1 x + a2 x2 + · · · + am xm . Maak een function file met als invoer de rijvectoren a = (a0 , a1 , · · · , am ) en x = (x1 , x2 , · · · , xn ) en als uitvoer de rijvector (p(x1 ), p(x2 ), · · · , p(xn )). Pas hiertoe de function file uit onderdeel (a) aan. (c) Maak een function file met als invoer α, β ∈ R en de rijvector a = (a0 , a1 , · · · , am ) en als uitvoer een plaatje van de grafiek van het polynoom p(x) = a0 + a1 x + a2 x2 + · · · + am xm op het interval [α, β]. Pas hiertoe de function file uit onderdeel (b) aan.

25

3.2

De while-loop

We bekijken in deze paragraaf nogmaals iteratieprocessen maar nu maken we gebruik van een while-loop in plaats van een for-loop. Met het commando: while cond ; ; ... ... ; end wordt de serie opdrachten ,, ..., uitgevoerd zolang voldaan is aan de voorwaarde cond. De voorwaarden zullen zijn van het type n < N, n = p, n 6= m, of combinaties ervan. De volgende tabel geeft een overzicht van de ‘logische operatoren’: n n n n n n

n n n n n n

> m < m >= m <= m == m ~= m

A & B of and(A,B) A && B A | B of or(A,B) A || B xor(A,B) ~A of not(A,B)

is is is is is is

groter dan m ; kleiner dan m ; groter dan of gelijk aan m ; kleiner dan of gelijk aan m ; gelijk aan m ; niet gelijk aan m ;

A en B ; A en B ; (Zie onderstaande opmerkingen.) A of B (evt. beiden) ; A of B (evt. beiden) ; (Zie onderstaande opmerkingen.) A of B, maar niet beiden ; (exclusieve ‘of’) ; ontkenning van A (‘niet A’).

Opmerkingen (i) De logische operatoren, die elementsgewijs werken, kunnen twee waarden opleveren: 0 en 1, 1 staat voor ‘waar’ en 0 staat voor ‘onwaar’. (ii) Wanneer A en B getallen zijn (0 of 1), dus géén matrix, en de situatie kan optreden dat A & B onafhankelijk is van B wordt A && B gebruikt. Zie ook de volgende tabel: Invoer A B 0 0 0 1 1 0 1 1

and A&B 0 0 0 1

or A|B 0 1 1 1

not ~A 1 1 0 0

xor xor(A,B) 0 1 1 0 26

Iets dergelijks geldt voor A | B. Zie ook de opmerkingen bij de volgende voorbeelden. Probeer eens uit: 2 == 3 waarde = [1 3 5] <= [1 2 6] x = sqrt(10); x > 2 & x < 3 (Kennis van B is nodig om iets te kunnen zeggen over het resultaat.) x = 1:10, w = ( sin(x) > 0.5 | sin(x) < -0.5 ) (A en B zijn matrices.) Opmerkingen (i) Een while-loop wordt gebruikt als je niet weet hoe vaak een loop doorlopen moet worden. ( Bouw de eerder gemaakte m-file om tot een function file met invoervariabele eps en uitvoervariabele x.) Voorbeeld x = 1; (De eerste benadering van while abs(x^ 2 - 2) > 1e-8 (eps = 10−8 ) x = x - (x^2 - 2)/(2*x); end;

√

2.)

(iii) Kom je er tijdens het rekenen achter dat er iets misgaat (met name bij een while-loop ligt het gevaar van een ‘oneindige’loop op de loer), dan kun je het rekenproces kun je stopzetten met de toetsen Ctrl + C of Ctrl + Break. Voorbeeld n = 1; while n > 0 n = n+1; end; Roep bovenstaande function file eens aan met eps = 10−20 . Je zult merken dat je wel erg lang moet wachten. Wat is de reden dat je in een ‘oneindige loop’ terecht bent gekomen ? Opgave 24. Hoeveel van de getallen sin(1), sin(2), sin(3), · · · sin(1000) zijn groter dan 0.5 ? (Antwoord: 332 )

27

Behalve met de zojuist genoemde ‘noodstop’, kun je met de volgende commando’s het rekenproces tijdens een loop in de gaten houden en be¨ınvloeden: toont de tekst tekst op het scherm; toont de inhoud (zonder de naam) van matrix X (X zal in het algemeen een ‘tekstvector’ zijn.); x = input(’tekst’) toont tekst op het scherm en wacht op een waarde van de gebruiker (gevolgd door Enter), en kent deze waarde aan x toe. (Indien de gebruiker direct Enter geeft dan wordt x de lege matrix.); str = input(’tekst’,’s’) toont tekst op het scherm en wacht op een antwoord van de gebruiker, bijvoorbeeld ja of nee (gevolgd door Enter) en kent dit antwoord aan str toe; [x,y] = ginput(n) (‘graphical input’) vraagt n punten (xi , yi ), die worden opgeslagen in de vectoren x en y, de gebruiker klikt de punten aan met de muis; [x,y] = ginput gaat door met punten op te slaan tot de gebruiker Enter geeft; pause onderbreekt een loop tot de eerstvolgende toetsaanslag (bijv. handig na disp ); pause(n) onderbreekt de loop n seconden (bijvoorbeeld om even een tussenresultaat te bekijken); break breekt een loop af; error(’tekst’) breekt een loop af en geeft tevens de (fout)melding tekst. disp(’tekst’) disp(X)

Hoe via het Command Window gecommuniceerd kan worden maakt het volgende voorbeeld van een m-file duidelijk: % De som van twee getallen ! a =’’; while (isempty(a))||(length(a) > 1) disp(’Geef a een waarde !’) a = input(’a = ’); end; b = ’’; while (isempty(b))||(length(b) > 1) disp(’Geef b een waarde !’) b = input(’b = ’); end; som = a + b; disp([’De som van a en b is :

’,num2str(som)]);

Als isempty(a) = 0 of isempty(b) = 0 hoeft de tweede voorwaarde niet meer gecontroleerd te worden!

28

3.3

Het if-statement

Met het commando: if cond ; else ; end; wordt eerst gekeken of aan de voorwaarde cond is voldaan. Zo ja, dan worden de opdrachten A uitgevoerd, zo nee, opdrachten B. Het gedeelte else ... mag ontbreken, en het mag ook op zijn beurt een ‘vertakking’ geven. De constructie gaat er dan als volgt uitzien: if condA ; elseif condB ; else ; end; Voorbeelden (i) Gegeven is een vector met 100 willekeurige kentallen in (0, 1 ). Het volgende ‘programma’ onderzoekt hoeveel van deze kentallen kleiner zijn dan 0.5. x = rand(1,100); for i = 1:100 if x(i) < 0.5 y(i) = 1; else y(i) = 0; end; end; s = sum(y); disp(s); Welk antwoord verwacht je ongeveer ?

29

(ii) En met het volgende ‘programma’ wordt gekeken of de elementen van vector x naar grootte geordend zijn : L = length(x); k = 1; Geordend = 1; while (Geordend && k < L) if x(k) > x(k+1) Geordend = 0; else k = k+1; end; end; disp(Geordend);

(x is geordend ⇔ Geordend = 1)

Opmerkingen (Voor diegenen die geprogrammeerd hebben in Pascal, Algol, ...) (i) In for- en while-loops geen do gebruiken, en bij de if ... else - constructie geen then. (ii) Bij een voorwaarde van het type ‘A is gelijk aan B’ A == B gebruiken, en niet A = B. Het laatste geeft de nogal vage foutmelding: missing variable or function (iii) Matlab kent ook een switch-constructie. Zo je wilt: raadpleeg help switch. Opgaven 25. Bekijk nogmaals het eerste voorbeeld uit §3.3. Herschrijf dit zodat gebruik gemaakt wordt van maximaal 4 programmeerregels. (Hint: bekijk de voorbeelden bovenaan blz. 27 nog eens goed.) 26. Bekijk nogmaals het eerste voorbeeld uit §3.2. Herschrijf dit zodat in plaats van een while-loop gebruik gemaakt wordt van een for-loop en een if-statement. 27. Bekijk de opgaven 15, 16 en 17 nog eens. Breid de eerder gemaakte m-files uit door er voor te zorgen dat bij onzinnige invoer foutmeldingen worden gegeven. 28. Schrijf een function-file kwadverg.m (met als invoervariabelen a, b, c) die de kwadratische vergelijking a x2 + b x + c = 0 voor alle waarden van a, b en c ∈ R oplost. 29. (Nulpuntsbepaling via bisectie) De functie f : R → R, met f (x) = x + exp(x) heeft een nulpunt p op het interval [ −2, 0 ]. Gevraagd wordt p te benaderen tot op 5 decimalen achter de komma, door herhaalde halvering van het interval [ −2, 0 ] tot een interval [ a, b ] van lengte kleiner dan 2·10−5 (a + b) wijkt dan ten hoogste 10−5 van p af. dat nulpunt p bevat. Het midden c = 2 30

Ga bijvoorbeeld als volgt te werk: stel om te beginnen a = −2, b = 0 en err = 2. Laat c het midden zijn van [ a, b ], bereken f (c), stel a = c dan wel b = c, al naar gelang f van teken wisselt tussen c en b, of tussen a en c. Halveer err en herhaal deze stappen totdat err kleiner is geworden dan 2·10−5 . 30. Neem een willekeurig positief geheel getal p groter dan 1. Als p even is: deel p door 2, als p oneven is: vermenigvuldig p met 3 en tel er 1 bij op. Herhaal dit proces tot het resultaat 1 is. Het is een open vraag uit de getaltheorie of dit proces altijd eindigt, c.q. het getal 1 altijd bereikt wordt. Twee commando’s voor ‘integer division’ (n en k zijn gehele getallen): r = rem(n,k) geeft de rest r bij deling van n door k: n = q k + r q = floor(n/k) geeft bovengenoemde q. (a) Na hoeveel stappen wordt 1 bereikt als p = 99 ? (Maak een M-file grow.m.) (b) Wat is de hoogste tussenwaarde die voor een startwaarde p = 99 bereikt wordt? (c) Wat is de hoogste tussenwaarde die voor een startwaarde p kleiner dan 100 bereikt wordt, en voor welke p gebeurt dit? 31. (‘bubble sort’) Schrijf een function file sorteer.m die de elementen van een vector x naar oplopende grootte sorteert, en wel op de volgende wijze: - ‘doorloop’ de vector x van voor naar achter, - en verwissel ‘buren’ die niet in de juiste volgorde staan, - houdt met een variabele change bij of er verwisselingen plaatsvinden. Zo ja, herhaal dan de stappen en zo nee, dan is x geordend. 32. Gegeven is het volgende iteratieproces :   f1 = 1 en f2 = 2 f + fn−1  fn+1 = n (n = 2, 3, · · · ) 2 Maak gebruik van de Matlab-commando’s disp en input bij het maken van een m-file iteration.m waarbij de invoer van f1 en f2 plaats vindt via het Command Window en waarbij na elke iteratie wordt gevraagd of de volgende iterand moet worden uitgerekend. Zo ja, dan moet dit worden berekend en zo nee, dan moet het programma worden afgebroken.

31

4

Driedimensionale plots

In dit hoofdstuk beschouwen we enkele commando’s om drie-dimensionale figuren te maken. Evenals bij het plot-commando (Zie § 1.8 ) worden daarbij vectoren (en matrices) als argument gebruikt. Het is daarbij belangrijk er op te letten dat de afmetingen goed (gelijk aantal elementen) zijn.

4.1

Ruimte-krommen

Met het plot3-commando kun je ruimte-krommen tekenen. Intikken plot3(x,y,z) tekent de punten (x(1), y(1), z(1)), (x(2), y(2), z(2)), . . . , (x(n), y(n), z(n)) en verbindt ze met rechte lijnen; x, y en z mogen hier rij- of kolomvectoren met een gelijk aantal elementen zijn. Voorbeeld Voer in:

t = 0:0.1:10*pi; x = cos(t); y = sin(t); z = t; plot3(x,y,z);

Het resultaat is een opwaartse spiraal. Door in het Plot menu op Tools te klikken, en vervolgens op Rotate3D, kun je het plaatje roteren m.b.v. de muis. Net als bij het plot-commando kun je de layout beinvloeden. Het instellen van kleuren en tekens gaat als beschreven in § 1.8, waar ook diverse hulpcommando’s staan vermeld. Speciaal voor 3D-plots zijn nog: Intikken axis([xmin xmax ymin ymax zmin zmax])

de ‘x-coördinaat’ loopt van xmin tot xmax, de ‘y-coördinaat’ van ymin tot ymax; de ‘z-coördinaat’ van zmin tot zmax;

zlabel(’tekst’)

er wordt tekst langs de z-as gezet.

Voorbeeld Voer in:

hold on plot3(x(1),y(1),z(1),’o’);

Het beginpunt van de spiraal wordt nu gemarkeerd met een ’o’. Nog een

32

Voorbeeld Voer in:

t = 0:0.1:5; x = 2*t; y = 3*t; z = 4*t - t.^2; plot3(x,y,z,’ro’); axis([0 10 0 15 0 10]);

De positie van de kromme in nog onduidelijk. Om een betere indruk te krijgen voeren we vervolgens in: grid Er worden nu in de figuur gestippelde roosterlijnen getekend. Opgaven 33. Construeer de (gladde) ruimte-kromme met parametervoorstelling x = (4 + sin(20t))cos(t), y = (4 + sin(20t))sin(t), z = cos(20t),

t ∈ [0, 2π].

Merk het beginpunt met een ‘*’. 34. Construeer de ruimte-kromme met parametervoorstelling x = cos(20t)sin(t), y = sin(20t)sin(t), z = cos(t),

t ∈ [0, π].

Merk op dat je voor t een geschikte interval-verdeling moet nemen om een gladde kromme te krijgen.

4.2

Oppervlakken

Met het surf-commando kunnen we ruimtelijke oppervlakken tekenen. Intikken surf(Z)

tekent het oppervlak bepaald door de punten (1, 1, Z(1, 1)), . . . , (1, m, Z(m, 1)), (2, 1, Z(1, 2)), . . . , (2, m, Z(m, 2)), . . . , (n, 1, Z(1, n)), . . . , (n, m, Z(m, n)) Z is hier een matrix met afmeting m × n;

surf(x,y,Z) tekent het oppervlak bepaald door de punten (x(1), y(1), Z(1, 1)), . . . , (x(1), y(m), Z(m, 1)), (x(2), y(1), Z(1, 2)), . . . , (x(2), y(m), Z(m, 2)), . . . , (x(n), y(1), Z(1, n)), . . . , (x(n), y(m), Z(m, n)) x en y moeten hier rij- of kolomvectoren zijn, en Z matrix met afmeting length(y) × length(x); surf(X,Y,Z) tekent het oppervlak bepaald door de punten (X(1, 1), Y (1, 1), Z(1, 1)), . . . , (X(m, 1), Y (m, 1), Z(m, 1)), (X(1, 2), Y (1, 2), Z(1, 2)), . . . , (X(m, 2), Y (m, 2), Z(m, 2)), . . . , (X(1, n), Y (1, n), Z(1, n)), . . . , (X(m, n), Y (m, n), Z(m, n)) X, Y en Z zijn matrices met gelijke afmeting. 33

Voorbeeld Voer in:

x = 0:0.1:2; y = 0:0.2:2; for i = 1:length(y) for j = 1:length(x) Z(i,j) = x(j)^2 + y(i)^2; end; end; surf(x,y,Z);

Het resultaat is een parabolo¨ıde. Opmerking De for-loop kun je ook vervangen door de opdracht Z = ones(length(y),1) * x.^2 + (y.^2)’ * ones(1,length(x)); Omdat zowel de for-loop als de laatste opdracht nogal wat hoofdbrekens kosten, is er het hulpcommando meshgrid, waarmee het automatisch goed gaat. Nogmaals het Voorbeeld Voer in:

x = 0:0.1:2; y = 0:0.2:2; [X,Y] = meshgrid(x,y); Z = X.^2 + Y.^2; surf(x,y,Z);

Een oppervlak dat gedefinieerd wordt door een parameter-voorstelling, waarbij de ‘zcoördinaat’ niet geschreven wordt als functie van x en y, wordt getekend met de matrixversie van surf. Voorbeeld Voer in:

phi = [0:100]*pi/100; theta = [0:200]*pi/100; X = sin(phi’)*cos(theta); Y = sin(phi’)*sin(theta); Z = cos(phi’)*ones(1,201); surf(X,Y,Z);

Het resultaat is een bol met straal 1. Merk op dat de 101 × 201-matrices X, Y en Z verkregen worden als product van een kolomvector met lengte 101 en een rij-vector met lengte 201. Opgaven 35. Teken het oppervlak gedefinieerd door z =

sin x sin y , xy 2 −y 2

36. Teken het oppervlak gedefinieerd door z = −xye−x 34

−10 ≤ x ≤ 10, −10 ≤ y ≤ 10. ,

−3 ≤ x ≤ 3, −3 ≤ y ≤ 3.

Met het commando subplot kun je meerder plaatjes in één window zetten. Dit gebruiken we om het verschil tussen surf en mesh te laten zien. Voorbeeld Voer in:

x = -1:0.1:1; y = x; [X,Y] = meshgrid(x,y); Z = sqrt(2 - X.^2 - Y.^2); subplot(1,2,1),surf(x,y,Z); subplot(1,2,2),mesh(x,y,Z);

Mesh tekent alleen de roosterlijnen, terwijl surf het oppervlak ook inkleurt. Vaak, zoals in de eerste voorbeelden, kunnen we het oppervlak opvatten als een functie van twee variabelen, z = f (x, y). Voor het opsporen van de (locale) extremen van zo’n functie zijn contourlijnen getekend in het (x, y)-vlak behulpzaam. Deze verkrijgen we met de commando’s surfc en meshc: Intikken surfc(x,y,Z) tekent het oppervlak bepaald door de vectoren x en y en de of meshc(x,y,Z) matrix Z, met tevens de contourlijnen in het (x, y)-vlak. Voorbeeld Voer in:

x = -1:0.1:1; y = x; [X,Y] = meshgrid(x,y); Z = 4-1./(1 + X.^2 + Y.^2); meshc(x,y,Z);

Het extreem (minimum) bevindt zich binnen het middelste cirkeltje. Indien de contourlijnen niet helemaal duidelijk zijn, doordat ze gedeeltelijk achter het oppervlak verscholen liggen, kunnen ze ook getekend worden met het commando Intikken contour(x,y,Z) tekent een 2D-plot van de contourlijnen van het oppervlak bepaald door de vectoren x en y en de matrix Z. Opgaven 37. Maak een plot van de contourlijnen van de functie uit de vorige opgave. Leid daaruit af in welke punten de functie een maximum aanneemt. Aanwijzing: Met het commando ginput(1) kun je de coördinaten verkrijgen van een punt in de grafiek dat je met de muis aanklikt. 38. Teken het oppervlak en de contourlijnen voor de functie z = x2 − 2xy + 4y 2 , waaruit blijkt dat het minimum optreedt in x = 0, y = 0.

35

4.3

Staven en taarten

Met de commando’s bar3 en pie3 kan de grootte van de elementen van matrices en vectoren grafisch worden weergegeven als staafdiagrammen of cirkelplots. Intikken bar3(Z)

geeft de grootte van de elementen Z(1, 1), Z(2, 1), . . . , Z(m, 1), Z(1, 2), Z(2, 2), . . . , Z(m, 2), . . . , Z(m, n) weer als staafdiagrammen, waarbij de elementen uit elke kolom van de matrix Z langs de y-as worden uitgezet.

pie3(x)

verdeelt een taart in stukken waarvan de grootte evenredig is met de elementen x(1), x(2), ..., x(n) van de vector x.

Voorbeeld Voer in:

4.4

Z = [1 2 3 6 4 1 2 3 4 5 6 7]; x=Z(2,:); subplot(1,2,1),bar3(Z); subplot(1,2,2),pie3(x);

Animatie

Een verkregen reeks plotjes kun je met de commando’s getframe en movie afdraaien als een filmpje. Om een goed resultaat te krijgen moet je er wel op letten dat de plotjes gelijke afmetingen hebben. Hiervoor kun je zorgen met het commando axis. Intikken F(j) = getframe slaat het huidige plotje op als j-de element van F movie(F,N)

vertoont de plotjes opgeslagen in F als filmpje, en draait dit filmpje N keer af.

Opmerking Het opslaan van F kost vrij veel geheugenruimte. Je dient hiervoor dan ook over een krachtige PC te beschikken, of anders data-compressie technieken toe te passen.

36

Voorbeelden (i) Voer in:

phi x = y = for

= 0:pi/100:2*pi; cos(phi); sin(phi); j = 1:20 r = j/2; x = r*cos(phi); y = r*sin(phi); plot(x,y); hold on; F(j) = getframe; end; movie(F);

Het resultaat is een filmpje waarbij steeds meer meer cirkels zichtbaar worden met de oorsprong als middelpunt. Matlab geeft een waarschuwing omdat niet vooraf is aangegeven hoeveel geheugenruimte er nodig is voor de opslag van de plaatjes uit het filmpje. Een ander voorbeeld maar nu wel met reservering van geheugenruimte ziet er als volgt uit : (ii) Voer in:

phi = 0:pi/100:2*pi; x = cos(phi); y = sin(phi); F(20).image = struct(’image’,[]); for j = 1:20 shiftx = 3*cos(j*pi/10); shifty = 3*sin(j*pi/10); plot(x+shiftx,y+shifty,’r’,... x/2+shiftx,y/2+shifty,’b’,3*x,3*y,’:k’); axis([-5 5 -5 5]); F(j).image = getframe; end; movie([F.image],10);

Het resultaat is een rode cirkel met een blauwe kern die 10 rondjes om de oorsprong draait.

37

Uiteraard kun je ook van 3D-plotjes een filmpje maken. (iii) Voer in:

x = -1:0.1:1; y = x; [X,Y] = meshgrid(x,y); Z = exp(-X.^2-Y.^2); F(20).image= struct(’image’,[]); for j = 1:20 shiftx = cos(j*pi/10); shifty = sin(j*pi/10); surf(x+shiftx,y+shifty,Z); axis([-2 2 -2 2 0 2]); F(j).image=getframe; end; movie([F.image]);

Het resultaat is een heuveltje dat om de oorsprong draait.

38

39

Referenties [1] D.C. Lay. Linear Algebra and its Applications, third edition, Update. Addison Wesley, U.S.A. 2003. [2] J. Stewart. Calculus, Early Transcedentals (International Student Edition), sixth edition. Thomson, U.S.A. 2008.

40

Cursus Matlab voor Technische Aardwetenschappen. K. Dekker I.A.M. Goddijn P. Sonneveld

Recommend Documents