Introductie in MATLAB

Introductie in MATLAB Andr´ e Hensbergen Inhoudsopgave Introductie in MATLAB, deel 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9

Invoeren van opdrachten De help-faciliteit Rekenen met getallen Toekenningen Werken met matrices Matrixbewerkingen; oefeningen Functies in Matlab Plotten Printen

Introductie in MATLAB, deel 2 2 3 3.1 3.2 4 5

.m-files Control Flow for-loops if-else en while Algemene commando’s Opgaven

Opdrachten

Introductie in MATLAB, deel 1 Matlab (afkorting voor Matrix Laboratory) 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 de ‘output’ op het scherm te onderdrukken sluit je een commando af met een puntkomma.

1.1

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 a+b a−b ab 3xy a b b

a √

x π e 3 − 4i sin x, arctan x, . . . ex , ln x log x ( = 10 log x ) |x| ∞

Intikken a+b a-b a*b 3*x*y a/b a^b sqrt(x) pi exp(1) 3-4*i sin(x), atan(x), ... exp(x), log(x) log10(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-functie. Met behulp van het commando help onderwerp verschijnt er een toelichting bij onderwerp. Bijvoorbeeld:  help elfun geeft informatie over de elementaire functies die Matlab kent. Lees (via ↑ ↓ in de rechterbalk of Page Up/Down) wat er zoal in staat. Met alleen het commando help verschijnt een lijst van mogelijke onderwerpen.

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. Tik in:  12/3

(Druk nu de Enter of RETURN–toets in)

Op het scherm verschijnt de uitvoer van het commando in de volgende vorm: ans = 4

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 ongeveer 16 decimalen nauwkeurig (waarvan er standaard vijf op het scherm verschijnen). Voer in:  12/9

(Druk nu de Enter of RETURN–toets in)

Als antwoord verschijnt op het scherm: ans = 1.3333

2

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

√ 5 Bereken (1 + 2i)(5 − 3i) en (1 + i 3)√ , √ 3 1 πi en ga met MATLAB na dat e 4 = − 2 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.1416 1.2346e–003 1.2346e–006

format long

15 cijfers, als short

3.14159265358979 0.00000123456000 1.23456000000000e–006

format zonder toevoegingen is hetzelfde als format short. Voor overige mogelijkheden: raadpleeg help format. Stel Matlab in op long. Voer in:  pi  sin(pi)

Het antwoord is (bijna!) 0.

Ga nu weer over op de standaard-instelling .

3

1.4

Toekenningen

Als een matrix (of getal) meerdere keren gebruikt gaat worden kun je de waarde ervan toekennen aan een variabele. Je hebt hem dan als het ware een naam gegeven, en kunt hem via die naam weer oproepen.   1 2 1 2 . Als voorbeeld berekenen we de eerste vier machten van de matrix  2 −1 1 3 −2  A = [1, 2, 1; 2, -1, 2; 1, 3, -2] nemen)  A2 = A^2; A3 = A^3; A4 = A^4;

(i.p.v. komma’s kun je ook spaties

Vervolgens kun je de som hiervan berekenen via  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 ´e´en 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. N.B.1 Als je niet zelf een variabele introduceert om een waarde in op te slaan, 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. N.B.2 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:     

who size(x) clear x1 clear all clc

opvragen van alle variabelen met een waarde; geeft de afmetingen (aantal rijen, kolommen) van x; verwijder de variabele x1 ; verwijder de waarden van alle variabelen. veegt het scherm schoon (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; (hiermee wordt x ‘leeg’ gemaakt)  functie = x^2 + 4 * sin(x) Om een dergelijke functie in te voeren moet deze door de gebruiker zelf ‘ingeprogrammeerd’ worden. De gebruiker breidt als het ware de lijst Matlab-commando’s uit met een nieuwe functie. In deze fase gaan wij hier niet verder op in. Het is (nog) 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: 4

 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.) Met de Delete- en de Back Space-toets (dat is bij sommige toetsenborden de ← toets boven Enter) kun je symbolen wegvegen; 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 ´e´en 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] ( spaties tussen de getallen! ) Op het scherm verschijnt: A= 1 4 7

2 3 5 6 8 9

Voorbeeld van het invoeren van een kolomvector:  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 de mogelijkheden: A(i,j) A(:,j) A(i,:) A(k:l, m:n) v(i:j)

Aij (het element in de ie rij, je kolom) j e kolom van A ie rij van A deelmatrix van alle elementen Aij k ≤ i ≤ l, m ≤ j ≤ n (oftewel: k-de t.e.m. l-de rij, m-de t.e.m. n-de kolom) ‘deelvector’ (vi , vi+1 , . . . , vj ) van de vector v

5

En voor een vierkante matrix: diag(A) geeft de hoofddiagonaal van A (als kolom-vector); diag(A,k) geeft de k-de nevendiagonaal van A (als k < 0, dan wordt degene onder de hoofddiagonaal genomen) Probeer het een en ander uit: Voer in de 3 × 4-matrix: A = [1 2 3 4; 5 6 7 8; 1 1 1 1] (spaties!) Voorspel en controleer wat het resultaat is van de volgende commando’s: (N.B.: ´e´en van de commando’s zal een foutmelding geven. Welke, en waarom? ) • 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) 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

6

1.6

Matrix-bewerkingen

Wat nu volgt is een greep uit het grote arsenaal aan matrix-bewerkingen waar de Matlab-gebruiker standaard de beschikking over heeft. De commando’s met een worden toegelicht in het volgende punt. De commando’s met een  hebben betrekking op begrippen die pas later (bij lineaire algebra) aan de orde zullen komen. Voor alle bewerkingen geldt: is het je niet duidelijk wat er bedoeld wordt: raadpleeg help. gewone optelling van matrices gewone verschil van matrices gewone matrixvermenigvuldiging (bij onjuiste afmetingen volgt een foutmelding); C = A.*B ‘elementsgewijze’ vermenigvuldiging (A en B moeten dezelfde afmetingen hebben); C = A+k bij elk element van A wordt het getal k opgeteld; C = A^k gewone machtsverheffing (k ∈ Z); bijv. te gebruiken voor berekening van A−1 ; C = A.^k ‘elementsgewijze’ machtsverheffing; (kan handig zijn bij het invoeren van een polynoom) C = A’ berekent de getransponeerde AT ; C = A\B berekent voor een niet-singuliere, vierkante matrix A de unieke oplossing van het de vergelijking AX= B; berekent voor een overbepaald stelsel (meer vergelijkingen dan onbekenden) een kleinste kwadraten oplossing; C= B/A oplossing van XA = B, analoog aan vorige commando; [L,U,P] = lu(A) berekent de PLU-ontbinding van A;  L = eig(A) L is een kolomvector met de (mogelijk complexe) eigenwaarden van A;  [X,D] = eig(A) D is een diagonaalmatrix met op de diagonaal de eigenwaarden van A, de kolommen van X zijn de bijbehorende eigenvectoren;  [Q,R] = qr(A) produceert een QR-ontbinding van matrix A; A = diag(v) geeft diagonaalmatrix 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 = k:n x wordt de vector (k, k+1, k+2, ..., n); x = a:h:b x wordt de vector (a, a+h, a+2h, ..., b), aangenomen dat b–a deelbaar is door h C = A+B C = A-B C = A*B

7

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 

   1 1 2 (iii) Laat A =  1 0  en b =  1  0 1 5

Merk op:

Ax = b is strijdig.

Bereken u = A\b, en ga na of Au = b. N.B. dit is geheel in overeenstemming met de uitleg bij het commando U= A\B (meer hierover: straks bij lineaire algebra)

Om het ‘denken in rijen en kolommen’ te bevorderen: Opgave Voer (op een handige manier) de volgende twee matrices in:   1 3 5 7 9  1 3 5 7 9     : : : : :  (A heeft 18 rijen. Gebruik de   A=  commando’s a:h:b en ones. )  : : : : :   : : : : :  1 3 5 7 9    B=  

1 2 3 1 4 9 1 8 27 1 16 81 1 32 ...

4 16 ... ... ...

... ... ... ... ...

9 92 93 94 95

     

Overigens zullen we later een andere manier zien een ‘gestructureerde’ matrix in te voeren, nl. met een for-loop.

8

Opgave Gegeven een rijvector v van lengte n. Hoe kun je de som van de elementen van v berekenen? Aanwijzing: neem het product met een geschikte kolomvector. Bereken op deze manier de som van de eerste honderd kwadraten, dat wil zeggen: 1 + 4 + 9 + . . . . + 992 + 1002 . (Antwoord: 338350)

1.7

Functies in Matlab

Zoals gezegd: Matlab giet alle structuren in de vorm van matrices. Zo √ ook functies. α x Matlab kent de elementaire functies x , sin, cos, tan, e , ln ( log !), x ( sqrt(x) ), en zo nog wat meer. (Als je wilt weten welke: vraag het op met help elfun.) 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.) Vraag: wat stelt y(5) voor ? De ‘elementsgewijze’ matrix-operaties .* en .^ zijn met name handig om ‘functies’ in te voeren: Voorbeeld De functie f (x) = 13 x2 + x sin x op het interval [0, 3] met ‘stapgrootte’ 0.1. Ten eerste de x-waarden:  x = 0:0.1:3.0; (dit geeft een vector van lengte 31) Nu moeten we voor k = 1, 2, 3, . . . , 31 berekenen y(k) = f (x(k)). Met behulp van de .* en de .^ operator 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. Het volgende punt betreft het grafisch weergeven van gegevens.

9

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: 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)), ... , (x(n), y(n)). x en y mogen hier rij- of kolomvectoren van gelijke lengte zijn. Voorbeelden Voer in x = [2,5,4,8,10,4] y = [6,3,9,1,4,0] plot(x,y) Voer in (bij het invoeren van tekst springt Matlab automatisch naar het Command Window terug): x = 0:0.1:3; y = sin(x); plot(x,y) (geen paniek, zie eerstvolgende opmerking! ) Opmerking Bij een volgende plot-opdracht komt het plaatje niet altijd direct in beeld. Als het Plot Window zich achter het Command Window bevindt, kun je het op verschillende manieren naar voren halen: A) als een deel van het Plot Window wel zichtbaar is, click je het simpelweg aan met de muis. (Je kunt een deel zichtbaar maken door het Command Window iets te verkleinen.) of B) door op de menubalk ”Windows” aan te clicken en vervolgens ”Figure No. 1”. of C) door te ‘bladeren’ via Alt-Tab.

Nog een voorbeeld: Voer in

x = 1:10; y = 20:30; plot(x,y)

Wat gaat hier fout?

10

Bij het plotten van meerdere grafieken in ´e´en 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 Symbool Stijl y geel (yellow) . ‘stip’ r rood o ‘rondje’ g groen ‘lijn’ b blauw -. ‘streep-stip’ w wit * ‘ster’ Voorbeeld Voer in: x = 0:.1:3.1; y = sin(x); z = x.*sin(x); (elementsgewijze matrixvermenigvuldiging!) plot(x,y,’y-’, x,z,’r-.’) Er zijn twee manieren om meerdere grafieken in ´e´en plaatje te krijgen: met ´e´en plot-opdracht, 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); plot(x,y,’g’)

(haal de vorige regel terug en wijzig)

Ga nu naar het Command Window terug. Het Plot Window is dan nog steeds aanwezig. Geef vervolgens: hold on plot(x2,y2,’yo’) en de twee grafieken staan in ´e´en 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 hold off te geven, of door het Plot Window te verlaten door de Control Menu Box (het hokje linksboven) te double-clicken.

11

Met onderstaande commando’s kun je de ‘layout’ van een bestaande plot be¨ınvloeden: axis([xmin xmax ymin ymax]) ‘x-co¨ordinaat’ loopt van xmin tot xmax, de ’y-co¨ordinaat’ van ymin tot ymax. axis(’off’)

tekent de assen niet

title(’tekst’)

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

xlabel(’tekst’)

zet tekst langs de x-as

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.

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. Oefeningen 1. Construeer een (gladde)

grafiek van sin x op het interval [0, 4π], waarin de kπ , sin , k = 0, 1, . . . , 24 met een * zijn aangegeven. punten kπ 6 6 2. (Krommen in het (complexe) vlak. Zie moduul 1, § 1.4) De kromme gedefinieerd door x = x(t), y = y(t), a ≤ t ≤ b, of door Φ : [a, b] → C kun je met Matlab eenvoudig op het scherm krijgen. We nemen als voorbeeld x = 3 cos(t), y = sin(t), t ∈ [0, 2π]. Voer in

12

   

t = 0:0.1:2*pi; x = 3*cos(t); y = sin(t); plot(x,y)

(een redelijk fijne interval-verdeling)

Nu zelf: (a) 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). (i) x(t) = cos(2t), y(t) = sin(5t) (ii) x(t) = cos(2t + π/3), y(t) = sin(5t) (iii) x(t) = cos(5t + π/4), y(t) = sin(4t) (b) Teken de kromme z = t eit , 0 ≤ t ≤ 4π. (In feite opgave 2 van § 4 uit moduul 1.) Aanwijzing: je kunt t en z vrijwel letterlijk invoeren. z wordt dan een vector van complexe getallen. Je krijgt de kromme in beeld met  plot(real(z), imag(z)) Geef ook eens de stijl-optie ‘punt’ of ‘ster’ mee; dan kun je als het ware zien hoe de ‘snelheid’ verandert

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 de menubalk van het Plot Window ”File” aan, en in de Dialog Box die verschijnt ”Print”. In het ”Afdrukken” Window tenslotte click je ”OK” (okay) aan. Als de opdracht is verwerkt kom je weer terug in het Plot Window. Een van de printers zal een plaatje uitspuwen (eventueel nadat jij op de knop continue reset hebt gedrukt, of vandaag of morgen weer een andere knop). [Printen vanuit Matlab geeft nogal eens complicaties. Voordat je bij weigering om te printen hulp inroept van een begeleider of systeembeheerder, kun je het volgende nog proberen: verlaat (bv. met Alt+Tab) Matlab, click in de Program Manager het Turbo Pascal ikoon aan, en vervolgens het ikoon Recapture(dosprint). Vorig jaar lukte het printen dan soms wel.]

13

Introductie in MATLAB, deel 2 2

.m-files

In de eerste introductie heb je kennisgemaakt met een aantal basiscommando’s van het pakket Matlab, alsmede met enkele van de plot-mogelijkheden. In deze tweede introductie leer je (een klein beetje) 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 met daaromheengebouwd een scala aan al dan niet specialistische procedures en functies, opgeslagen in zogenoemde M-files (of .m-files). Voorbeelden hiervan: 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 toevoegen. Die optie maakt dragelijk dat het Matlab Command Window nogal beperkt is in z’n mogelijkheden. In feite is Matlab nog steeds een DOS-programma met een Windows-jasje. Commando’s (opnieuw) laten uitvoeren door ze aan te clicken werkt niet, en het telkens terughalen van eerder ingetypte commando’s gaat bij een wat ingewikkelder probleem al snel irriteren. Een Matlab-sessie opslaan om er op een later moment mee verder te gaan heeft ook geen zin.

Script Files Ten eerste kun je een aantal commando’s in een zogenaamde script file zetten (oftewel een programma’tje schrijven) en die vervolgens ‘runnen’. Ten tweede kun je function files maken met procedures en functies, die je vanuit het command window aan kunt roepen. Werkend op een computer in de pc-zaal ga je eerst via het commando cd a: naar de a-drive (nadat je hebt gezorgd dat zich daarin een schijfje bevindt). Het is namelijk niet de bedoeling allemaal prive-spullen op het netwerk op te slaan. Voor een script file ga je als volgt in te werk (eerst a.u.b. cd a:) — Click op de menubalk aan: File

New

M-file

— De Matlab Editor/Debugger opent een nieuw document, waarin jij een programma kunt schrijven. (In oudere versies van Matlab wordt de notepad-editor gebruikt.) — Het is een goede gewoonte om het programma te beginnen met een of enkele regels commentaar over de werking van het programma. (Een commentaar wordt voorafgegaan door %). Met help naam (de nog te verzinnen naam) krijg je dit commentaar later snel in beeld. 14

Ben je klaar met schrijven dan moet je het ‘programma’ (dat nu alleen nog maar op het scherm staat) ergens opslaan, in jargon: ‘saven’. — Click aan: File Save As Geef het beest een naam: naam.m gevolgd door Enter of Ok De extensie .m is belangrijk! (Als jij hem niet geeft, dan voegt de editor die standaard toe.) — Verlaat de .m-file. (Alt-F4, of Alt-Tab, of File Exit) — Laat vanuit het Matlab Command Window de commando’s uitvoeren, door simpelweg naam te geven (zonder .m), of via: File Run naam.m — Wil je het programma wijzigen, dan roep je de M-file op door: File Open naam.m (of door te bladeren met Alt-Tab, als je de editor aldus verlaten had.) Om de gewijzigde versie te bewaren geef je Save in plaats van Save As (tenzij je de oude versie ook nog wilt bewaren natuurlijk). — Zoals de naam al aangeeft kun je de Matlab-editor ook gebruiken om fouten in je programma’s op te sporen en te verwijderen (= ‘debuggen’ ). Indien nodig moet je dit t.z.t. zelf maar zien uit te vogelen. Voorbeeld (Tevens oprakelen van plotten) Maak een .m-file met de volgende commando’s: 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*’) Save het geheel (op de a: onder de naam plotsin.m (of iets anders, mits nieuw en eindigend op .m), ga naar het Command Window en geef het commando 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.

15

Function Files Deze paragraaf is ter extra informatie aangeboden. Voor de opdrachten hoef je niet per se zelf function-files te maken (hoewel het handig kan zijn). Sla hem in eerste instantie over, en lees verder bij: Control Flow. Zoals gezegd, er kunnen twee soorten .m-files onderscheiden worden: script files en function files. In script files wordt simpelweg een aantal commando’s in een .m-file gegroepeerd. Met function files kun je zelf functies en procedures aan Matlab toevoegen. Script files 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. Met help functienaam kun je dit commentaar later eenvoudig opvragen. Enkele voorbeelden: 2

1. De functie f : R → R met functievoorschrift f (x) = x2 e−x kun je toevoegen als .m-file newfunc: 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]. 2. 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 kun je voortaan het ”maximale verschil tussen twee opeenvolgende elementen van een vector y” laten berekenen via  s = maxdiff(y) Doe dit voor de vector z gedefinieerd door  x = 0:0.1:1.5*pi;  z = sin(x.^2);

16

Opmerkingen: 1. Bij ‘zinvoller’ function files zul je vaak gebruikmaken van de commando’s for, if, while, . . . . 2. Voor ´e´en-regelige functies als in voorbeeld 1. 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 newfunc(5) kunt berekenen via  newfunc(5) Raadpleeg desgewenst help inline. 3. De variabelen result en x in de regel function result = functienaam(x) zijn, hoe kan het anders in Matlab, matrices. x heet de input, result de output. In plaats van een enkele variabele, kun je ook meerdere variabelen als input en/of output nemen. De aanhef wordt dan function [res1, res2, ..., resk] = functienaam(x1, x2, ..., xn) Tot slot: 4. Een belangrijk verschil tussen script files 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.

17

3

Control Flow

3.1

(for, while, if...)

for-loops

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 loop-variabele k. (Spreek uit: ”loep-variabele”) Het is een goede gewoonte om de loop-opdrachten in te laten springen, en niet meer dan ´e´en opdracht per regel te schrijven.

Voorbeelden 1) Met de volgende commando’s kun je de som van de eerste 50 kwadraten berekenen (maak een .m-file som.m): s = 0; for k = 1:50 s = s + k^2; end; s De opdracht die 50 keer wordt uitgevoerd is: tel k^2 op bij s. Het laatste commando, zonder puntkomma, doet het antwoord op het scherm verschijnen. Voor wie nog nooit geprogrammeerd heeft: kijk eens wat er gebeurt als in de derde regel de puntkomma weggelaten wordt.

18

2) De rij van Fibonacci 1, 1, 2, 3, 5, 8, 13, . . . . . is vastgelegd door de eerste twee waarden en het voorschrift fn+1 = fn + fn−1 Het volgende ‘programma’ genereert de eerste vijftig Fibonacci-getallen. f(1) = 1; f(2) = 1; for k = 3:50 f(k) = f(k-1) + f(k-2); end; Aan het eind is f een vector van lengte 50, 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): f1 = 1; for k = f3 = f1 = f2 = end; f2

f2 = 1; 1:48 f1+f2; f2; f3;

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) tussen: disp([f1, f2, f3]), pause(2)

Het is ook mogelijk om binnen een for-loop een andere for-loop te plaatsen. We spreken dan van een geneste loop. Voorbeeld clear all for k = 1:8 for l = 1:8 A(k,l) = k+l-1; end; end; A

(twee keer end!) (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.)

19

Opmerking: Vermijd for-loops als er een (effici¨enter) 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 Zoals in de eerste introductie beschreven, kun je de som van de eerste N kwadraten ook berekenen met de commando’s: 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. Neem om de gedachte te bepalen N = 10000. Voeg in de m-file som.m (van voorbeeld 1) als eerste regel toe: tic en als laatste regel toc en verander de waarde 50 door de waarde 10000. Voer het programma uit. Vergelijk dit met het volgende ‘programma’: N= 10000; tic; somkwad = ((1:N).^2)*ones(N,1), toc Voor de volledigheid (en misschien begin je er ook wel lol in te krijgen) noemen we de while- en de if-statement. Heb je er geen lol in c.q. tijd voor, dan is dit voor jou het eind van de introductie. Blader door naar de pagina met kop foutmeld.m, en zet vervolgens je tanden in de opgaven en opdrachten.

20

3.2

while-loops en if ...

else

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

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

A & B A | B xor(A,B) ~A

n n n n n n

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 of B (evt. beiden) (logische ‘of’) A of B, maar niet beiden (exclusieve ‘of’) ontkenning van A (‘niet A’)

De logische operatoren, die elementsgewijs werken, kunnen twee waarden opleveren: 0 en 1. 1 staat voor ‘waar’ en 0 staat voor ‘onwaar’. Probeer eens uit: 2 == 3 waarde = [1 3 5] <= [1 2 6] x = sqrt(10); x > 2 & x < 3 x = 1:10, w = ( sin(x) > 0.5 | sin(x) < -0.5 ) Oefening Hoeveel van de getallen sin(1), sin(2), sin(3), . . . . . , sin(1000) zijn groter dan 0.5 ? (Antwoord: 332 )

21

Opmerking Met een while-loop kan een for-loop nagebootst worden: k = 1; while k < n ; k = k+1; end; heeft dezelfde werking als: for k = 1:n ; end;

N.B. 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 stopzetten met de toetsen [Ctrl]+C of [Ctrl]+[Break]. (Voorbeeld: n = 1; while n > 0 n = n+1; end; )

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: disp(’tekst’) toont de tekst tekst op het scherm. disp(X) 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.) [x,y] = ginput(n) (‘graphical input’) vraagt n punten (xi , yi ), die worden opgeslagen in de vectoren x en y. De gebruiker clickt 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.

22

if ... else 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: cond.A elseif cond.B else end; if

; ; ;

Voorbeelden Het volgende ‘programma’ slaat de kwadraten onder de duizend op in de vector kw. n = 1; while n^2 < 1000 kw(n) = n^2; n = n+1; end; kw En met het volgende 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; (x is geordend ⇐⇒ Geordend = 1)

23

Opmerkingen (Voor diegenen die geprogrammeerd hebben in Pascal, Algol, ...) 1. In for- en while-loops geen do gebruiken, en bij de if ... else - constructie geen then. 2. 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 3. Matlab kent ook een case-constructie. Zo je wilt: raadpleeg help case.

24

4

Algemene commando’s

Nog even op een rijtje gezet een aantal commando’s om een soepele omgang met Matlab te bewerkstelligen.

maakt het command window leeg (maar de toegekende waarden blijven actief) clg leegt alle plot-windows diary bewaart de tekst van een Matlab-sessie save fname var slaat de waarde van de variabele var ‘extern’ op in de file fname of fname.mat load data.dat laadt variabele data (of data.mat) cd naam directory naam wordt de nieuwe working directory dir geeft de inhoud van de directory waarin u werkt help comm geeft hulp bij het commando comm what geeft de M-, Mat- en MEX-files uit de huidige directory lookfor naam doorzoekt de eerste ‘comment’-regels van alle M-files. Voorbeeld: lookfor zeros who geeft lijst met actieve variabelen in de huidige Matlab-sessie whos als who maar met extra informatie clear variabele verwijdert genoemde variabele clear all verwijdert alle variabelen disp ‘tekst’ toont tekst op het scherm (de commando’s T = ’tekst’; disp(T) geven hetzelfde resultaat) format stelt de vorm van de output in (opties: long, short, bank, compact, loose) clc

Verder kun je vanuit Matlab externe commando’s uitvoeren door ze vooraf te laten gaan door een uitroepteken. Voorbeeld (Als je in een DOS-omgeving werkt)

!copy file1 file2

!cd en !dir kunnen wel eens problemen geven, vermoedelijk omdat cd en dir ook Matlab-commando’s zijn. Tenslotte: is de hier gegeven uitleg te summier, dan kan de help-functie wellicht uitkomst bieden.

26

5

Oefeningen

Voor oefening 1 t.e.m. 4 die van harte aanbevolen zijn, heb je alleen for-loops nodig. Bij oefening 6 t.e.m. 8, bedoeld voor liefhebbers, is het gebruik van if en while aan te raden.

1. Voer de volgende matrices A en B in met behulp van ´e´en of meer loops.       A=      

1 1 1 1 : : : 1

3 3 3 3 : : : 3

5 5 5 5 : : : 5

...... ...... ...... ...... ...... ...... ...... ......

15 15 15 15 : : : 15

             (15 rijen),           

      B=     

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

2 4 8 16 32 ... ... ...

3 9 27 81 ... ... ... ...

4 16 ... ... ... ... ... ...

... ... ... ... ... ... ... ...

9 92 93 94 95 ... ... 99

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

2. n! (spreek uit ”n faculteit”) is gedefinieerd via: n! = n(n−1)(n−1) . . . . ·3·2·1 (en apart: 0! = 1 ). 20 X 1 1 : S= Bereken de som van de eerste eenentwintig getallen n! n! n=0 (Antwoord (via format long ): S = 2.71828182845905 )

3. 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 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 v als k-de nevendiagonaal, en met verder overal 0-en. (b) Via enkele loops. (Begin met het invoeren van een 30×30 nulmatrix.) (Een derde mogelijkheid is: met een dubbele loop waarin via if verschillende gevallen onderscheiden worden.) Bekijk ter controle de 5 × 5 rechter-onder-matrix, en bereken de determinant van A. (Antwoord: 31)

27

4. In hoofdstuk 3 moduul 1 is als beste eerste-orde-benadering van een functie de linearisering ingevoerd. In moduul 2 bekijken we hogere orde benaderingen, en dan blijkt de beste benadering van sin x met polynomen van graad ≤ 2n+1 gegeven te worden door het polynoom 3 5 2n+1 t2n+1 (x) = x − x + x − ... + (−1)n x 3! 5! (2n + 1)! In deze opgave bekijken we grafieken van tn voor verschillende n. Neem om te beginnen n = 4 en interval I = [−5, 5].

Voer in:

x = -5:.1:5; t9 kun je natuurlijk invoeren via: t9 = x - x.^3/6 + x.^5/120 - x.^7/(120∗6∗7) + x.^9/(120∗6∗7∗8∗9) maar voor grotere n wordt zoiets al snel ondoenlijk. We kiezen een slimmere aanpak, gebruik makend van de ’gelijkvormigheid’ van de termen die gesommeerd moeten worden: term term term term term

= = = = =

x; t9 = term; -term.∗x.^2/(2∗3); -term.∗x.^2/(4∗5); -term.∗x.^2/(6∗7); -term.∗x.^2/(8∗9);

t9 t9 t9 t9

= = = =

t9 t9 t9 t9

+ + + +

term; term; term; term;

De laatste vier regels lenen zich prima voor het gebruik van een for-loop! (i) Voer bovenstaande uit. Let goed op hoe je de loop-variabele, zeg k gebruikt (met name wat de begin- en eindwaarde moeten zijn). (ii) Bekijk sin x en t15 (x) op het interval [-5, 5] en op het interval [-8, 8]. Bekijk ook de verschilfunctie om te zien hoe groot de afwijking maximaal is. (iii) Bekijk sin x en t25 (x) op het interval [-10, 10] en op het interval [-12, 12]. Mooi h`e!

Terzijde Een handig commando om polynomen in te voeren is polyval: Laat v = [v1, v2, v3, ... , vk] en x vectoren zijn. Dan is y = polyval(v,x) equivalent met y = v(k) + v(k-1)∗x + v(k-2)∗x.^2 + ... + v(1)∗x.^(k-1) Hetgeen overigens veel effici¨enter berekend wordt via: y = (( ... (((v(1)*x+ v(2)).*x + v(3)).∗ .....

).*x + v(k)

(Ter vergelijking: 2x3 + 5x2 − 3x + 7 = ((2x + 5)x − 3)x + 7 )

28

5. (Lagrange-polynomen) Deze opgave geeft nog wat extra oefenmateriaal. Op de volgende pagina vind je een korte uitwerking. Gegeven n + 1 verschillende getallen x0 , x1 , x2 , ..., xn in R, en n + 1 niet noodzakelijk verschillende y-waarden y0 , y1 , y2 , ..., yn . We willen een polynoom van zo laag mogelijke graad berekenen dat door alle n + 1 punten (x0 , y0 ), (x1 , y1 ), . . . . , (xn , yn ) gaat. We gaan daartoe als volgt te werk: Het gezochte polynoom p moet voldoen aan de n + 1 eisen p(xi ) = yi , dus het ligt voor de hand om te proberen: p(x) = a0 + a1 x + a2 x2 + . . . + an xn . Invullen van de n + 1 gegeven punten levert voor de co¨effici¨enten a0 , ... , an het volgende stelsel vergelijkingen:  a0 + a1 x0 + a2 x20 + ... + an xn0 = y0     a0 + a1 x1 + a2 x21 + ... + an xn1 = y1    : : : : : :     : : :    a0 + a1 xn + a2 x2n + ... + an xnn = yn 

Oftewel:

      

1 1 1 : : 1

x0 x1 x2 : : xn

x20 x21 x22 : : x2n

    a0 y0 ..... xn0     ..... xn1    a1   y 1  n    ..... x2   a2   =  y2 .     :  :   :   :  :   :  ..... xnn an yn

Aangetoond kan worden dat de co¨effici¨entenmatrix van dit stelsel niet singulier is, als de ‘steunpunten’ x0 , x1 , x2 , ..., xn verschillend zijn. Neem bijv. de zes punten (-2, 3), (-1, 5), (0, 2), (1, 2), (2, 0) en (4, -2). Gezocht is nu het polynoom p5 (x) van graad ≤ 5 door deze punten. (i) Sla de x-waarden op in x, en de y-waarden in y en plot de zes punten met de plot-optie *. (ii) Voer de co¨eff. matrix A op een handige manier in. Controleer voordat je gaat rekenen of de gegevens goed zijn ingevoerd! (iii) Bereken de co¨effici¨enten a0 , ...., a5 via het commando u = A\y. (y moet een kolom zijn.) (iv) Vorm nu het polynoom p5 , waarvan de co¨effici¨enten in u zijn opgeslagen. (Zie opg. 4) (v) Bekijk het resultaat: hold on; plot(t,p5)

29

6. (‘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 nee, dan is x geordend, zo ja, herhaal de stappen. Overigens, Matlab bevat een veel sneller sorteercommando: sort.

7. (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 dat nulpunt p bevat. Het midden c = (a + b)/2 wijkt dan hoogstens 10−5 van p af. Ga bijvoorbeeld als volgt te werk: stel om te beginnnen a = -2; b = 0; 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 . Ook hier: Matlab heeft ingebouwde commando’s roots (voor nulpunten van polynomen) en fzero.

8. 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?

30

6

Opdrachten

In het verleden was er een verplicht Matlab-practicum waarin iedere student (onder anderen) vijf van de volgende opdrachten moest beantwoorden en inleveren. Ze zijn toegevoegd als extra (en eenvoudig) oefenmateriaal. In ´e´en opdracht wordt een print (of twee prints) van een kromme gevraagd. Bij de andere opdrachten is het voldoende de antwoorden in ´e´en .m-file te zetten en die vervolgens net als een plot uit te printen (nl. door aanclicken van File Print ). Hiertoe copieer je eerst de relevante zaken vanuit het Command Window naar de ‘verslag .m-file’. Dit gaat als volgt: — zet de mouse-pointer aan het begin van het te copi¨eren gedeelte; — houd de linkerknop van de muis ingedrukt terwijl je de mouse-pointer naar het einde van het te copi¨eren gedeelte schuift. Laat daar los; Het ‘ge-highlighte’ (zwartgemaakte) gedeelte kan gecopi¨eerd worden door — Edit Copy aan te clicken (of Edit Cut als je de tekst op de oorspronkelijke positie weggegooid wilt hebben.), — de muis naar de positie (eventueel in een ander window) te verplaatsen waar je de tekst wilt hebben, en vervolgens — Edit Paste te geven. (De zogenaamde shortcut keys Ctrl-C, Ctrl X en Ctrl-V hebben hetzelfde effect als Edit Copy, Edit Cut resp. Edit Paste.) In elke opdracht geldt: c1 c2 c3 c4 c5 c6 is je zes-cijferige studienummer. Met een parameter als c2 c3 c4 wordt bedoeld het getal met de drie cijfers c2 , c3 en c4 . Op de volgende wijze bepaal je jouw vijf opdrachten: Laat vk = ck + 2ck+1 , k = 1, 2, 3, 4, 5 en rk = vk (mod 4) (de rest bij deling door 4; dus rk = 0, 1, 2, of 3). Het k-de opdrachtnummer wordt: 5rk + k. Ingewikkeld? De volgende commando’s doen precies wat nodig is:  c = [c1 c2 c3 c4 c5 c6];  v = c(1:5) + 2*c(2:6); r = rem(v,4);  nrs = 5*r + (1:5) en de vijf opdrachtnummers verschijnen in beeld.

32

Opdrachten 1. A is de 24 × 24 matrix met in de eerste rij en in de eerste kolom c1 c2 c3 c4 c5 c6 c1 c2 ..... c6 , op de diagonaal c1 , 1, −1, 1, . . . . , −1, 1 en verder overal 0-en. (Handig invoeren die hap!) Bereken de determinant van A. (Vind zelf het commando om de determinant te berekenen. Vergeet bij dit commando niet de haakjes.)

2. Laat a1 = c1 c2 c3 en a2 = c4 c5 c6 . Laat verder, voor n > 2 an+1 = 32 an + 25 an−1 . Bereken a99 en a100 tot drie cijfers achter de komma. Let bij gebruik van een loop op de onder- en de bovengrens!

3. Construeer als volgt vector v: laat n = 10 · (c1 + c2 + c3 + c4 ), h = π/n, dan wordt v = (sin(h), sin(2h), sin(3h), . . . . . , sin(nh) ). Bereken p = v1 · v2 + v2 · v3 + . . . . + vn−1 · vn en ook de maximale co¨ordinaat van w = ( |v1 · v2 |, |v2 · v3 |, . . . , |vn−1 · vn | ).

4. Neem p = 1.c1 c2 , q = 2.c3 c4 . (Bijv. bij studienr. 123456: p = 1.12, q = 2.34 ) Maak een tabel van de functies f (x) = sin(px) en g(x) = x cos(qx) op het interval [0, 6] met stapgrootte 0.3. M.a.w.:

0.0 f (0.0) g(0.0) 0.3 f (0.3) g(0.3) .. .. .. .. .. .. .. .. .. .. .. .. 6.0 f (6.0) g(6.0)

5. Laat a het grootste cijfer zijn van je studienummer en b het kleinste. Laat c = a − 2, en d = b + a + 1. Geef in ´e´en plaatje de volgende twee Lissajouxkrommen weer:   x(t) = cos(at) x(t) = cos(at + π/d) en y(t) = sin(ct) y(t) = sin(ct) Neem het interval zo groot dat er de krommen gesloten zijn, en gebruik verschillende lijn-stijlen.

33

6. Gevraagd wordt van de volgende matrix A de determinant te berekenen. A is een 24 × 24 matrix met op de diagonaal c1 c2 c3 c4 c5 c6 c1 c2 . . . . . c6 , direct boven en onder de diagonaal 1-en en verder overal 0-en. (Handig invoeren die hap!) (Vind zelf het commando om de determinant te berekenen. Vergeet bij dit commando niet de haakjes.) 7. Neem α = 3.c2 . De rij hxn√ in is gedefinieerd door x0 = c3 c4 en x1 = c5 c6 , en voor n ≥ 1: xn+1 = α · 3xn + 2xn−1 . Bereken x14 en x15 tot op 14 decimalen nauwkeurig. Let bij gebruik van een loop op de onder- en de bovengrens!

8. Construeer als volgt vector v: laat n = 10 · (c1 + c2 + c3 + c4 ), h = π/n, dan wordt v = (cos(h), cos(2h), cos(3h), . . . . . , cos(nh) ) Bereken A = (v1 − v2 )2 + (v2 − v3 )2 + . . . . , (vn−1 − vn )2 , alsmede de maximale co¨ordinaat van w = ( |v1 − v2 |, |v2 − v3 |, . . . , |vn−1 − vn | ).

9 Neem p = 1.c1 c2 , q = 0.c3 c4 c5 . (Bijv. bij studienr. 123456: p = 1.12, q = 0.345 ) Maak een tabel van de functies f (x) = xe(−px) + x en g(x) = eqx /x op het interval [1, 5] met stapgrootte 0.2. M.a.w.:

1.0 f (1.0) g(1.0) 1.2 f (1.2) g(1.2) .. .. .. .. .. .. .. .. .. .. .. .. 5.0 f (5.0) g(5.0)

10. Laat a het grootste cijfer zijn van je studienummer en b het kleinste. r en β = R + r . Laat R = a+b, r = a–b–2 , α = R − r r (Als r toevallig kleiner zou worden dan 3: neem r = 3.) Geef in ´e´en plaatje de volgende ’Spirograaf-figuren’ weer  x(t) = (R − r) cos(t) + r cos(αt) K1 : y(t) = (R − r) sin(t) + r sin(αt)  x(t) = (R + r) cos(t) + r cos(βt) en K2 : y(t) = (R + r) sin(t) + r sin(βt) Neem het interval zo groot dat de krommen gesloten zijn, zorg dat de verhouding x:y gelijk is aan 1:1, en gebruik verschillende lijn-stijlen.

34

11. A is de 24 × 24 matrix met op de diagonaal c1 c2 c3 c4 c5 c6 c1 c2 . . . . . c6 , direct boven en onder de diagonaal 2-en en verder overal 1-en. (Handig invoeren die hap!) Bereken de determinant van A. (Vind zelf het commando om de determinant te berekenen. Vergeet bij dit commando niet de haakjes.) 12. Laat a1 = 0.c3 c4 c5 , en an+1 = 3.5an (1 − an ), n ≥ 1. Bereken a74 en a75 tot 10 cijfers achter de komma. Let bij gebruik van een loop op de onder- en de bovengrens! 13. Construeer als volgt vector v: laat n = 10 · (c1 + c2 + c3 + c4 ), h = 3/n, dan wordt v = (cos(h), cos(2h), cos(3h), . . . . . , cos(nh) ) v2 v3 vn Laat w = ( , , . . . , ). v1 v2 vn−1 Bereken A = (w1 )2 + (w2 )2 + . . . . . (wn−1 )2 en ook het minimum M van {w1 , w2 , w3 , . . . , wn−1 } 14. Neem p = 1.c1 c2 , q = 0.c3 c4 c5 . Maak een tabel van de functies f (x) = x2 e(−px) + x en g(x) = eqx /(x + 1) op het interval [0, 4] met stapgrootte 0.2. M.a.w.:

0.0 f (0.0) g(0.0) 0.2 f (0.2) g(0.2) .. .. .. .. .. .. .. .. .. .. .. .. 4.0 f (4.0) g(4.0)

15. We bekijken de kromme die een punt doorloopt dat meebeweegt met een wiel dat over een horizontale weg loopt. Het wiel heeft straal R, en het punt begint zich op afstand d van de as. (d kan groter zijn dan R, denk bijvoorbeeld aan een treinwiel over een rail, en als d = R dan doorloopt het punt een cyclo¨ıde.) Een  parametrisering voor de kromme x(t) = Rt − d sin(t) y(t) = R − d cos(t) Laat R = 1, d1 = 34 R en d2 = b ·R, waarbij b = 1 + (c1 + c2 )/40. Geef in ´e´en plaatje de twee krommen weer met de gegeven waarden voor R en d, waarbij het wiel twee volledige omwentelingen maakt. Gebruik verschillende lijn-stijlen. Maak ook een plot waarbij de assen zodanig gekozen zijn dat de verhouding x:y ongeveer 1:1 is.

35

16. A is de 3×24 matrix met als eerste rij vier keer (naast elkaar) c1 c2 c3 c4 c5 c6 , als tweede rij [ 1 2 3 . . . . 24 ], en als derde rij [ 25 26 . . . 48 ]. (Handig invoeren die hap!) Bereken de determinanten van de matrices AAT en (AT A − I24 ). (Vind zelf het commando om de determinant te berekenen. Vergeet bij dit commando niet de haakjes.) 17. Laat a1 = c1 c2 c3 en a2 = c4 c5 c6 . Laat verder, voor n > 2 an+1 = 21 an + 12 an−1 + 10. Bereken a99 en a100 tot drie cijfers achter de komma. Let bij gebruik van een loop op de onder- en de bovengrens!

18. Construeer als volgt vector v: laat n = 2(c1 + c2 + c3 + c4 + c5 + c6 ), n 2 ) ) dan wordt v = (cos((1.1)2 ), cos((1.2)2 ), cos((1.3)2 ), . . . . . , cos((1 + 10 Laat w = (w1 , w2 , . . . , wn−1 ) = ( |v1 − v2 |, |v2 − v3 |, . . . , |vn−1 − vn | ). Bereken A = (w1 )2 + (w2 )2 + . . . . . + (wn−1 )2 en het maximum M van {w1 , w2 , w3 , . . . , wn−1 }.

19. Neem p = 1.c1 c2 , q = 0.c3 c4 c5 . (Bijv. bij studienr. 128756: p = 1.12, q = 0.875 ) Maak een tabel van de functies f (x) = x sin(px) en g(x) = cos(qx2 ) op het interval [0, 6] met stapgrootte 0.3. M.a.w.:

0.0 f (0.0) g(0.0) 0.3 f (0.3) g(0.3) .. .. .. .. .. .. .. .. .. .. .. .. 6.0 f (6.0) g(6.0)

c4 20. Laat a = 14 + 20 en b = 3 + 12 (c5 + c6 ). Teken de kromme gedefinieerd door z = e(a+bi)t , waarbij −2 ≤ t ≤ 4. Bereken de minimale waarde r en de maximale waarde R van |z|, en teken in hetzelfde plaatje de cirkels met middelpunt (0, 0) en stralen r resp. R.

36