Handleiding Matlab voor W en BMT

Handleiding Matlab voor W en BMT bij de eerstejaars-training Matlab 2000/2001

Voorwoord Voor u ligt de trainings-handleiding Matlab voor het eerste studiejaar W en BMT. De handleiding bestaat uit vier delen, geschreven door de verschillende betrokken docenten. De bijbehorende training omvat zes dagdelen gedurende het eerste studiejaar. Doel van deze dagdelen is het aanleren van basisvaardigheden met betrekking tot het gebruik van Matlab als rekengereedschap bij het oplossen van technische problemen. De opgedane vaardigheden worden reeds bij een aantal vakken en cassussen in het eerste jaar benut. Zo heeft b.v. de training in trimester 1.1 als belangrijke “afnemer” het vak Calculus, terwijl de training in trimester 1.2 vooral het vak Lineaire Algebra bedient. Hieronder volgt een overzicht van de trainings-dagdelen. Dagdeel 1 2 3 4 5 6W 6 BMT

Trimester 1.1 1.1 1.2 1.2 1.3 1.3 1.3

Docent Martens Martens Van der Meer Van der Meer Van Essen Haitjema Van den Bosch

Handleiding Deel 1, H1,H2,H3 Deel 1, H4 Deel 1, H5 Deel 1, H6 Deel 2 Deel 3 Deel 4

Het eerste dagdeel vormt een inleiding in het gebruik van Matlab en een verkenning van de beschikbare numerieke gereedschappen. Het tweede dagdeel is een kennismaking met de symbolische rekenmogelijkheden van Matlab. Gedurende het derde dagdeel wordt geoefend in het oplossen van lineaire-algebra-problemen. Het vierde dagdeel besteedt aandacht aan het numeriek oplossen van differentiaalvergelijkingen. Het vijfde dagdeel is gericht op het aanleren van basisvaardigheden in het programmeren. Ten slotte maken de W-studenten gedurende het zesde dagdeel kennis met de statistische functies van Matlab t.b.v. het verwerken van meetgegevens en maken de BMTstudenten kennis met het simuleren van dynamische systemen. De dagdelen hebben het karakter van een Begeleide Zelfstudie. Het materiaal in deze handleiding kan zelfstandig worden doorgewerkt. Daarbij is het efficiënt om in groepen van twee personen te werken. Gedurende de dagdelen zal steeds gelegenheid zijn om vragen te stellen aan de docent. Aan het einde van de dagdelen twee, vier en zes vindt een korte individuele toets plaats, waarvan de uitwerking wordt ingeleverd bij de docent. Deze uitwerkingen dienen als basis voor een individuele beoordeling van de training als geheel (zie de regeling Matlabvaardigheid). Deelname aan alle dagdelen is een verplicht onderdeel van het eerstejaars curriculum. Gestreefd wordt naar een optimale afstemming van de trainingen op de curricula W en BMT. Hierbij is uw ervaring met de verschillende dagdelen in combinatie met deze handleiding van groot belang. Opmerkingen en suggesties t.a.v. de training kunt u richten tot ondergetekende, bij voorkeur via E-mail: [email protected] René van de Molengraft Matlab-coördinator W/BMT augustus 2000

Regeling Matlabvaardigheid De training omvat zes dagdelen (ochtenden of middagen). Op dagdeel 2, 4 en 6 vindt gedurende het laatste half uur een korte individuele toets plaats m.b.t. de stof van de twee voorgaande middagen. De student levert een schriftelijke uitwerking van de toets in bij de docent. De docent classificeert deze als voldoende/onvoldoende. Uiteraard is het gebruik van notebook en Matlab-handleiding tijdens de toets toegestaan. Niet toegestaan is gebruik van het computer-netwerk t.b.v. icq, irc, netmeeting, email e.d. Indien de student voor elk van de drie deeltoetsen een voldoende scoort, is zijn vaardigheid Matlab voldoende. Indien één of meerdere toetsen onvoldoende zijn, is zijn vaardigheid Matlab onvoldoende. Deeltoetsen kunnen niet afzonderlijk worden herkanst. Aan het eind van het jaar vindt een herkansingstoets over alle zes dagdelen plaats. Deze is bestemd voor studenten met een onvoldoende eindresultaat en voor studenten die niet aan elk van de drie deeltoetsen hebben meegedaan. Bij voldoende resultaat van de herkansingstoets wordt de vaardigheid Matlab van de student alsnog als voldoende aangemerkt. Het gezamenlijke resultaat van de drie deeltoetsen alsook het resultaat van de herkansingstoets behoren tot OGO-1.3. De beoordelingen van de eerste twee deeltoetsen maken geen deel uit van OGO-1.1 en OGO-1.2 en blijven derhalve buiten beschouwing bij de toelating tot OGO-1.2 en OGO-1.3. Indien de vaardigheid Matlab als onvoldoende beoordeeld is, zullen de studiepunten van OGO 1.3 niet worden toegekend. De student zal dan ook niet worden toegelaten tot OGO-2.1.

Deel 1 Handleiding Matlab voor W en BMT trimester 1.1 en 1.2, dagdeel 1 t/m 4 F.J.L. Martens J.C. van der Meer Faculteit Wiskunde en Informatica 2000/2001

Handleiding MATLAB W en BMT trimester 1.1., 1.2.

F.J.L. Martens J.C. van der Meer Faculteit Wiskunde en Informatica

2000/2001

2

Inhoudsopgave 1 Deze handleiding

5

2 Basiselementen van MATLAB

7

1

Het pakket MATLAB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

2

Het starten en stoppen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

3

Opdrachten in MATLAB . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

4

Het afbreken van berekeningen of van uitvoer . . . . . . . . . . . . . . . . . .

10

5

Onvolledige opdrachten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

6

Variabelen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

7

Het opslaan van databestanden . . . . . . . . . . . . . . . . . . . . . . . . . .

12

8

Het opslaan van de tekst van een MATLABsessie . . . . . . . . . . . . . . . .

12

9

Het beheer van les en het zoekpad . . . . . . . . . . . . . . . . . . . . . . . .

13

10

Wiskundige expressies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14

11

De arraystructuur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

12

Arrays met een rij of een kolom . . . . . . . . . . . . . . . . . . . . . . . . . .

17

13

Array-operaties in MATLAB . . . . . . . . . . . . . . . . . . . . . . . . . . .

18

13.1

Algebrasche array-operaties . . . . . . . . . . . . . . . . . . . . . . . .

18

13.2

Relationele array-operaties . . . . . . . . . . . . . . . . . . . . . . . .

21

13.3

Andere array-operaties . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

13.4

Transponeren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24

14

Gra eken . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

15

Script les . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

16

Eigen functies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

3

4

INHOUDSOPGAVE

17

De helpfaciliteiten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

18

Opgaven . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

3 De numerieke gereedschapskist

33

1

Inleiding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

2

De gra ek . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

3

De nulpunten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

4

De extrema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

5

Numerieke integratie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

6

Opgaven . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

4 De symbolische gereedschapskist

39

1

Inleiding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

39

2

Expressies met variabelen . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

40

3

Substitutie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41

4

Dierentieren en integreren . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

5

Numerieke waarden . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

43

6

Het manipuleren van expressies . . . . . . . . . . . . . . . . . . . . . . . . . .

44

7

Symbolen, strings en getallen . . . . . . . . . . . . . . . . . . . . . . . . . . .

45

8

Opgaven . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

47

5 Lineaire Algebra

51

1

Inleiding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

51

2

Matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

51

2.1

Wat is een matrix? . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

51

2.2

Het invoeren van matrices . . . . . . . . . . . . . . . . . . . . . . . . .

52

2.3

Matrices met symbolische elementen . . . . . . . . . . . . . . . . . . .

53

2.4

Het invoeren van vectoren . . . . . . . . . . . . . . . . . . . . . . . . .

53

2.5

Bijzondere matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . .

55

2.6

Indices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

56

2.7

Matrixbewerkingen in MATLAB . . . . . . . . . . . . . . . . . . . . .

57

INHOUDSOPGAVE

5

2.8

Matrixfuncties in MATLAB . . . . . . . . . . . . . . . . . . . . . . . .

60

Oplossen van stelsels lineaire vergelijkingen . . . . . . . . . . . . . . . . . . .

61

3.1

Het commando A n b . . . . . . . . . . . . . . . . . . . . . . . . . . . .

61

3.2

De rijgereduceerde trapvorm . . . . . . . . . . . . . . . . . . . . . . .

62

3.3

Het oplossen van stelsels met de Symbolic Toolbox . . . . . . . . . . .

63

Numerieke aspecten bij gebruik van MATLAB . . . . . . . . . . . . . . . . .

64

4.1

Getalrepresentatie en afrondfouten, een voorbeeld . . . . . . . . . . .

64

4.2

Datafouten, een voorbeeld . . . . . . . . . . . . . . . . . . . . . . . . .

65

4.3

EÆcientie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

66

MATLAB en Lineaire Algebra . . . . . . . . . . . . . . . . . . . . . . . . . .

67

5.1

67

3

4

5

Opdrachten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6 Dierentiaalvergelijkingen 1

2

3

4

69

M- les . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

69

1.1

M- les maken . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

69

1.2

Script- en function les . . . . . . . . . . . . . . . . . . . . . . . . . . .

69

Dierentiaalvergelijkingen . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

71

2.1

Het numeriek oplossen van dierentiaalvergelijkingen . . . . . . . . . .

71

2.2

Stelsels dierentiaalvergelijkingen . . . . . . . . . . . . . . . . . . . . .

73

2.3

Het richtingsveld . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

74

2.4

Plotten van integraalkrommen . . . . . . . . . . . . . . . . . . . . . .

75

Het symbolisch oplossen van dierentiaalvergelijkingen . . . . . . . . . . . . .

75

3.1

Eerste-orde dierentiaalvergelijkingen . . . . . . . . . . . . . . . . . .

75

3.2

Stelsels eerste-orde dierentiaalvergelijkingen . . . . . . . . . . . . . .

76

3.3

Hogere-orde dierentiaalvergelijkingen . . . . . . . . . . . . . . . . . .

77

MATLAB en dierentiaalvergelijkingen - Opdrachten . . . . . . . . . . . . . .

77

4.1

Opdrachten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

77

Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

80

Referenties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

80

6

INHOUDSOPGAVE

Hoofdstuk 1

Deze handleiding Deze handleiding is bedoeld om u snel met enkele aspecten van MATLAB vertrouwd te maken. Het is de bedoeling dat u, verdeeld over twee trimesters in een viertal dagdelen met dit materiaal oefent. Doe dit samen met iemand anders, dat werkt eÆcienter. Als u weinig met computers gewerkt heeft, dan zult u misschien niet alles in deze twee dagdelen af kunnen werken. Haal dit snel op een geschikt moment in! Het programma is als volgt:

Trimester 1.1. dagdeel 1 De hoofdstukken \De handleiding", \Basiselementen van MATLAB" en \De numerieke gereedschapskist". dagdeel 2 De hoofdstukken \De symbolische gereedschapskist".

Trimester 1.2. dagdeel 1 Het hoofdstuk \Lineaire Algebra". dagdeel 2 Het hoofdstuk \Dierentiaalvergelijkingen". De trainingen aan het begin van trimester 1.1 sluiten aan bij het college Calculus (2Y130) en de basisvaardigheden uit deze trainingen zullen in dit college, met name bij de oefeningen, verder worden ontwikkeld. Het eerste dagdeel in het tweede trimester sluit aan op het eerste college Lineaire Algebra (2Y650) en is de aanzet tot het gebruik van MATLAB bij lineaire algebra. Het tweede dagdeel in het tweede trimester sluit aan op het zevende college Lineaire Algebra (2Y650) en de colleges Dynamische systemen (4Q220) en Dynamica (4A230) en is een eerste aanzet voor het gebruik van MATLAB bij colleges in het derde trimester. MATLAB speelt een zeer belangrijke rol bij het OGO, het Onwerp Georienteerd Onderwijs. Hierbij zult u MATLAB moeten gebruiken. Eventueel moet u zichzelf met behulp van de 7

8

HOOFDSTUK 1. DEZE HANDLEIDING

MATLABdocumentatie commando's eigen maken om bepaalde berekeningen uit te kunnen voeren. Daarom is een goede beheersing van de elementaire MATLABcommando's essentieel. Iedere sectie is voorzien van oefeningen. Voer deze oefeningen nauwgezet uit en schrijf

uw uitwerkingen op. Vermeld bij uw uitwerkingen de gebruikte MATLABopdrachten. Als u geen aantekeningen maakt, zijn de oefeningen volstrekt zinloos!

Het is verstandig om uw aantekeningen te bewaren zodat u in de toekomst nog eens kunt naslaan hoe iets ook alweer moet. Aan het eind van een trimester moet u zo bedreven zijn in het gebruik van MATLAB, dat u het deel van de handleiding, dat bij de genoten trainingen hoort, niet meer hoeft in te kijken. De hulpfaciliteiten van MATLAB bieden ook alle informatie. Een belangrijk onderdeel van deze trainingen is het leren gebruiken van deze faciliteiten. We gaan ervan uit dat u met Windows95/98 vertrouwd bent. U moet in staat zijn om ` les' te beheren. U moet ze kunnen aanmaken, vinden, veranderen en weggooien. Op uw machine moet MATLAB inclusief de Èxtended Symbolic Math Toolbox' of de `Symbolic Math Toolbox' en een Internetbrowser staan. De versies van deze programma's doen minder ter zake. Een commando kan dus net iets anders werken dan in deze handleiding staat, maar met wat experimenteren komt u er vast wel uit. Er zijn enkele redenen om niet al te gedetailleerd te omschrijven hoe u MATLAB op moet starten, welke menukeuzes u moet maken, wat u moet invoeren en hoe de resultaten eruit zien. De eerste reden is dat op uw machine MATLAB op een andere manier genstalleerd kan zijn en dat uw versie net even anders is dan de versie die uitgangspunt is geweest bij deze handleiding. De tweede reden is dat bij iedere nieuwe versie op detailpunten de programma's veranderen. De helpfaciliteiten helpen u wel op de goede weg.

Hoofdstuk 2

Basiselementen van MATLAB 1 Het pakket MATLAB MATLAB is een pakket voor technische berekeningen. Met het pakket kunnen berekeningen, visualisaties en en zelfgeschreven programma's gecombineerd worden. MATLAB kan worden uitgebreid met zogenaamde `toolboxes'. In het algemeen zijn dit collecties MATLABfuncties die geschikt zijn voor bijzondere klassen van problemen. De \(Extended) Symbolic Math Toolbox" wijkt af. Deze `toolbox' bestaat voor een groot deel uit ingredienten van het pakket Maple dat sterk is in formulemanipulatie. MATLAB voert berekeningen uit met behulp van matrices. Matrices zijn rechthoekige schema's van getallen en worden ook arrays genoemd. De \Symbolic Math Toolbox" wijkt nogal af van MATLAB zelf. Symbolische berekeningen maken nauwelijks gebruik van arrays. MATLABfuncties maken in tegenstelling van functies uit de \Symbolic Math Toolbox" gebruik van berekeningen met matrices. Om deze reden is het van belang om te weten of een gebruikte functie uit MATLAB of uit de \Symbolic Math Toolbox" is. Op voorhand zal dit niet duidelijk zijn, maar na geregeld gebruik van MATLAB zal dit duidelijk worden.

2 Het starten en stoppen Starten kan op verschillende manieren:

Klik op `Start' en `Run' en type in de box de naam van de directory, waarin \matlab.exe" is opgeslagen, gevolgd door \matlab . Bij versie 5.3 van MATLAB is dit matlabr11\bin\matlab .

Klik het MATLABicoon aan.

In beide gevallen verschijnt het `Command Window' waarin opdrachten kunnen worden ingetypt na de MATLABprompt 9

10

HOOFDSTUK 2. BASISELEMENTEN VAN MATLAB

>>

Opdrachten worden afgesloten met de `return'-toets. De opdracht wordt dan door MATLAB uitgevoerd. Als MATLAB klaar is, dan verschijnt na eventuele uitvoer een nieuwe prompt. Er kan dan weer een opdracht gegeven worden. Het verlaten van MATLAB kan ook op verschillende manieren:

Door het commando >> quit

Via `FilenExit MATLAB'

Oefening 2.1 Start MATLAB, bekijk de mededelingen op het scherm en verlaat het pro-

gramma.

3 Opdrachten in MATLAB Na de MATLABprompt >> met de knipperende verticale streep, de cursor |, kan een opdracht worden ingetikt. Na indrukken van de `Return'-toets wordt de opdracht uitgevoerd en verschijnt het resultaat eventueel op het scherm. Een opdracht en zijn resultaat ziet er als volgt uit: >> a = 1 + 2 + 3 a = 6

Het resultaat van 1 + 2 + 3 wordt aan de variabele a toegekend. Met deze variabele kan men verder rekenen. >> b = 2*a + a/3 b = 14

In MATLAB worden variabelen gentroduceerd door er een waarde aan toe te kennen. Men kan dan met de variabele verder rekenen en er later eventueel weer een nieuwe waarde aan toekennen. Het is niet mogelijk om berekeningen uit te voeren met een variabele waaraan

3. OPDRACHTEN IN MATLAB

11

geen waarde is toegekend. Een opdracht hoeft niet met een toekenning van de vorm variabele = te beginnen. In zo'n geval wordt het resultaat automatisch aan de variabele ans toegekend. Men kan dan met ans verder rekenen. >> 4*a + 1 ans = 25 >> ans*ans ans = 625

Wilt u nu weten wat de waarde van a is, dan volstaat de volgende opdracht: >> a a = 6

MATLAB zet normaal de uitvoer onder de opdracht. De uitvoer wordt onderdrukt door achter de opdracht een puntkomma te plaatsen. Na uitvoering van de opdracht s = 1 + 2; ziet het scherm er als volgt uit >> s = 1 + 2; >>

De variabele s heeft de waarde 3 gekregen. Als u de waarde van deze variabele wil weten, dan kunt u s opnieuw opvragen. >> s s = 3

Als een opdracht meer dan een regel beslaat, dan kunt u voor het eind van de regel drie punten neerzetten en de `Return'-toets indrukken. U kunt dan op de volgende regel verder gaan. Het scherm ziet er na het voltooien van de opdracht als volgt uit:

12


>> s = 1 + 2 + ... 3 + 4 s = 10

Oefening 2.2 Het is de bedoeling dat u onderstaande opdrachten in de aangegeven volgorde uitvoert. Bedenk, alvorens een opdracht uit te voeren, wat er op het scherm verschijnt en wat het resultaat van de opdracht zal zijn. a. b. c. d. e.

>> a = 1 + 2 + 3 + 4 + 5 + 6 >> b = 1 + 2 + 3 + 4 + 5 + 6; >> b + 3 >> 1 + 2 + 3 + 4 + 5 + 6 >> c = 10*ans

4 Het afbreken van berekeningen of van uitvoer Een langdurige berekening van MATLAB kan met de toetscombinatie `Control-C' worden afgebroken. Na deze combinatie komt er een nieuwe prompt. Nieuwe opdrachten kunnen weer worden uitgevoerd.

Oefening 2.3 Evaluatievan de opdracht pause(30) heeft als eect dat MATLAB 30 secon-

den pauze houdt. Voer deze opdracht uit en pas snel de toetscombinatie `Control-C' toe. Wat is het eect?

5 Onvolledige opdrachten Als een opdracht onvolledig is en u drukt op de `Return'-toets, dan zijn er twee mogelijkheden:

Er verschijnen foutmeldingen en een nieuwe MATLABprompt. U kunt weer opnieuw aan de slag.

De knipperende cursor staat helemaal links aan het begin van de regel onder de opdracht. Er zijn dan twee mogelijkheden.

{ De opdracht kan worden afgemaakt en vervolgens worden geevalueerd. { De toetscombinatie `Control-C' zorgt ervoor dat de invoering van de opdracht wordt afgebroken. Er verschijnt ook een nieuwe MATLABprompt.

13

6. VARIABELEN

Oefening 2.4 Door de opdracht a

toegekend.

(a) Voer a

= [2,3,4] wordt aan de variabele a het array (2; 3; 4)

= [2,3,4 in en druk op de `Return'-toets.

Maak de opdracht af.

(b) Voer a

= [2,3,4 in en druk op de `Return'-toets. Onderbreek de invoering van de opdracht.

6 Variabelen De naam van een variabele moet met een letter beginnen. Daarna mag de naam uit een willekeurig aantal cijfers en/of letters bestaan. MATLAB maakt onderscheid tussen grote en kleine letters. MATLAB gebruikt enkele speciale variabelen: ans eps i j pi Inf NaN

Deze variabele bevat de uitkomst van de laatste berekening die niet aan een andere variabele is toegekend. De waarde van deze variabele is ongeveer 2:2204 10 16 en wordt intern door MATLAB gebruikt. Verwar eps niet met de e-macht. Want e = exp(1) 2:7183 . Het complexe getal i met de eigenschap i2 = 1. Het complexe getal i. De notatie j voor het getal i komt veel voor. 3:14159 : : : De waarde is oneindig. Delen van 1 door 0, 1/0, levert Inf op. Deze variabele is een representatie van \Not a Number". Een berekening met NaN resulteert altijd in NaN. De opdracht 0/0 levert NaN op.

Het commando workspace resulteert in een `window' met de naam \Workspace Browser" waarin alle door u zelf gebruikte variabelen vermeld staan. Er staat ook aangegeven tot welke klasse de variabelen behoren. Eventueel kunt u via dit `window' enkele variabelen verwijderen.

Oefening 2.5 Gebruik in enkele opdrachten de variabelen u en v correct en verwijder deze

vervolgens.

Enkele andere belangrijke commando's voor het variabelenbeheer zijn: who whos clear clear x y

geeft lijst met gebruikte variabelen. geeft lijst met gebruikte variabelen en extra informatie. verwijdert alle variabelen. verwijdert de variabelen x en y.

Waarschuwing: Het is mogelijk om aan de MATLABvariabelen zelf een andere waarde toe te kennen. Berekeningen kunnen hierdoor benvloed worden.

Geef nooit een andere waarde aan een MATLABvariabele.

Als u per ongeluk aan een MATLABvariabele een andere waarde geeft, dan kunt u deze waarde verwijderen met behulp van clear of de de \Workspace Browser".

14


7 Het opslaan van databestanden Resultaten van berekeningen van MATLAB worden vaak elders weer opnieuw gebruikt. De waarden van variabelen kunnen bewaard worden in een `.mat le'. Dit gebeurt met de opdracht >> save filenaam

waardoor alle gebruikte variabelen worden bewaard in de le met naam ``filenaam.mat''. De extensie `.mat' hoeft niet opgegeven te worden. Hoeven er maar een aantal variabelen te worden opgeslagen, dan dienen die gespeci ceerd te worden na de lenaam. De opgeslagen le kan in MATLAB worden ingelezen met de opdracht >> load filenaam

waarbij de extensie `.mat' mag worden weggelaten. De variabelen met hun waarde kunnen nu weer gebruikt worden.

Oefening 2.6 Voer de twee opdrachten >> f = 2

en

>> g = 3 + f

uit. Sla de variabelen f en g op in de le \probeer.mat". Verlaat MATLAB, start MATLAB opnieuw op en laad de le \probeer.mat". Kijk of de variabelen f en g bekend zijn. Waar is de le \probeer.mat" terecht gekomen? Verwijder de le \probeer.mat".

8 Het opslaan van de tekst van een MATLABsessie Het is ook mogelijk om alle in- en uitvoer van een MATLABsessie in een tekst le op te slaan. U kunt deze tekst niet voor een nieuwe sessie gebruiken, maar u kunt wel zien wat u allemaal gedaan heeft. De opdracht >> diary filenaam

zorgt ervoor dat alle in- en uitvoer naar de le \ lenaam" wordt weggeschreven. Met de opdracht >> diary off

wordt het opslaan beeindigd. Na het beeindigen kunt u de le \ lenaam" bekijken.

9. HET BEHEER VAN FILES EN HET ZOEKPAD

15

Oefening 2.7 Open een tekst le \afschrift" om in- en uitvoer op te slaan. Voer de opdrachten >> f = 2*3

en

>> g = f+2

uit en beeindig het opslaan van in- en uitvoer. Bekijk de inhoud van de le \afschrift". Waar is de le terecht gekomen? Verwijder deze le.

9 Het beheer van les en het zoekpad Als MATLAB draait, dan is er een directory actief. Aangemaakte les in een MATLABsessie komen hier terecht tenzij er expliciet directories worden vermeld. In MATLAB kunt u opvragen welke directory actief is. >> pwd ans = C:\MATLABR11\work

De opdracht is een afkorting van \path working directory". Het antwoord is natuurlijk afhankelijk van de installatie van MATLAB op uw machine en van de manier van opstarten. Er kan een andere directory, bijvoorbeeld \C:ntmp", gekozen worden waarin de aangemaakte les terechtkomen. >> cd C:\tmp

Bij het zoeken naar les doorloopt MATLAB een zoekpad. De directories waarin zelf gemaakte les staan, moeten in het zoekpad staan, want anders kan MATLAB de zelf gemaakte les niet terugvinden. Maak een directory \D:nmatlabtr" op uw D-schijf. De opdracht >> matlabpath

laat een lijst van alle directories zien waar MATLAB naar les zoekt. Het toevoegen gaat als volgt: >> addpath D:\matlabtr

Het opnieuw opvragen van het zoekpad laat zien of het toevoegen gelukt is. Het opstarten van MATLAB in de directory \D:nmatlabtr" en het toevoegen van deze directory aan het zoekpad kan worden geautomatiseerd door een le \startup.m" in die directory

16


te plaatsten, waar ook al de le \matlabrc.m" staat. Deze le wordt dan bij het opstarten van MATLAB automatisch gelezen en de opdrachten erin worden uitgevoerd. >> which matlabrc C:\MATLABR11\toolbox\local\matlabrc.m

In dit geval moet \startup.m" in \C:nMATLABR11ntoolboxnlocal" worden gezet. Maak een tekst le \startup.m" met de tekst % Voeg mijn directory met eigen M-files toe aan zoekpad. addpath D:\matlabtr % Ga naar mijn eigen directory. cd D:\matlabtr % Zet formaat voor papier op A4. set(0,'DefaultFigurePaperType','a4') % Gebruik centimeters in plaats van inches. set(0,'DefaultFigurePaperUnits','centimeters')

en zet deze le op de juiste plaats. Een regel achter een %-teken is commentaar. Het gebruik van deze le zorgt ervoor dat ook meteen enkele papierinstellingen worden veranderd. Verlaat MATLAB en start MATLAB opnieuw op. Controleer uw instellingen met onderstaande opdrachten. >> pwd ans = D:\matlabtr >> get(0,'DefaultFigurePaperType') ans = A4

Vanaf nu worden tijdens een MATLABsessie zelf gemaakte les zonder plaatsaanduiding automatisch naar \D:nmatlabtr" geschreven en ook daar gezocht!

10 Wiskundige expressies Het invoeren van wiskundige expressies gaat niet letterlijk. De volgende tabellen maken duidelijk welke vertaalslag u moet uitvoeren. De variabelen in deze tabellen stellen getallen voor.

17

10. WISKUNDIGE EXPRESSIES

Bewerkingen standaard MATLAB a+b a + b a b a - b ab a * b a b ab

a / b a ^ b

Opmerking De volgorde van de bewerkingen in MATLAB is de gewone volgorde: Eerst machtsverheen, dan vermenigvuldigen en delen en vervolgens optellen en aftrekken. Constanten Standaard MATLAB e exp(1) pi i i of j 1 inf of Inf Functies Standaard MATLAB Standaard MATLAB px sin(x) sin(x) sqrt(x) x cos(x) cos(x) e exp(x) tan(x) tan(x) ln(x) log(x) 10 arcsin(x) asin(x) log(x) log10(x) arccos(x) acos(x) jxj abs(x) arctan(x) atan(x) sign(x) sign(x)

Opmerking De ronde haken `(' en `)' worden gebruikt om de volgorde van de berekeningen vast te leggen.

Oefening 2.8 Bereken in MATLAB de volgende expressies: (a)

3 + 22 16 + 52

(b) 42=3 (c) log(e) (d) sin( ) 4

Oefening 2.9 Geef de variabele x de waarde 2. Voer de expressies aan de linkerkant in MATLAB in en bereken deze. Het resultaat van de berekeningen staat aan de rechterkant.

18


Expressie x3 6 x p 2 1+x 1 23 1+x2

e

x3 sin(x2 ) arctan(x) 1 + x2

Resultaat 1.3333 0.8944 1.2599 148.4132 -6.0544 0.2214

11 De arraystructuur MATLAB gebruikt als basiseenheid voor het rekenen matrices ofwel arrays. Een array is een rechthoekige schema van getallen, die elementen worden genoemd, gerangschikt in m rijen en n kolommen: 0 B B B @

a11 a12 : : : a1n a21 a22 : : : a2n .. .. .. .. . . . . am1 am2 : : : amn

1 C C C A

:

Een array kan op verschillende manieren worden ingevoerd. 1.

Tussen \[" en \]" waarbij de elementen gescheiden worden door spaties of komma's en de rijen gescheiden worden door een puntkomma. De opdrachten zijn >> a = [1 2 3;4 5 6;7 8 9]

of

>> a = [1,2,3;4,5,6;7,8,9]

2.

Tussen \[" en \]" waarbij de elementen gescheiden worden door spaties of komma's en de rijen steeds op een nieuwe regel met behulp van de `Return'-toets worden ingevoerd. De opdracht luidt nu

>> a=[1 2 3 4 5 6 7 8 9]

Het resultaat is dan a =

EN RIJ OF E EN KOLOM 12. ARRAYS MET E

1 4 7

2 5 8

19

3 6 9

Met arrays kan men als volgt manipuleren:

MATLAB

Betekenis van opdracht

a(1,2)

Geef het element uit de eerste rij en de tweede kolom van a

a(1,:)

Geef de eerste rij van a.

a(:,3)

Geef de derde kolom van a.

a(:,2)=[10;10;10]

Verander de tweede kolom van a in een kolom met alleen maar tienen. Geef het 6e element van de rij, ontstaan uit het aan elkaar plakken van de kolommen van a.

a(6)

12 Arrays met een rij of een kolom Een array (matrix) met een rij wordt ook wel een rij of een rijvector genoemd, en een array (matrix) met een kolom een kolom of een kolomvector. Als het verschil tussen de opeenvolgende elementen van de rij steeds hetzelfde is, kan de rij ook gevormd worden door `variabele' = `beginwaarde' : `stapgrootte' : èindwaarde' Hierbij is `beginwaarde' het eerste element van de rij, èindwaarde' het laatste element van de rij en `stapgrootte' het verschil tussen de opeenvolgende elementen. Bijvoorbeeld >> x = 1:-.2:0

geeft de rij x = 1.0000

0.8000

0.6000

0.4000

0.2000

0

Als het verschil tussen de opeenvolgende elementen van de rij een is kan de `stapgrootte' worden weggelaten. Bijvoorbeeld >> x=1:5

geeft als resultaat de rij

20


x = 1

2

3

4

5

Men kan arrays (matrices, vectoren) samenvoegen tot nieuwe arrays. Rangschikt men de arrays in een rij, dan moet het aantal rijen in de arrays overeenkomen. Rangschikt men de arrays in een kolom, dan moet het aantal kolommen in de arrays overeenkomen. Het aaneenplakken van de arrays a = [2 , 3 , 4] en b = [6 , 5] gaat als volgt: >> [a,b] ans = 2

3

4

6

5

Elementen van rijen en kolommen kunnen worden aangegeven door te verwijzen naar de index. In dit geval kan men met een index volstaan.

MATLABopdracht

Resultaat

b(1)

6, het eerste element van b.

a([3,1])

[ 4 2 ], een array bestaande uit het 3e en 1e element van a.

Oefening 2.10 Maak de rijvector met beginwaarde -1, eindwaarde 9 en stapgrootte 0.5 . Oefening 2.11 Maak de rijvector met beginwaarde 9, eindwaarde 0 en stapgrootte 1 .

13 Array-operaties in MATLAB Een array is een verzameling getallen gerangschikt in rijen en kolommen. Een \arrayoperatie" is een operatie (zoals de optelling of de vermeingvuldiging) die men op overeenkomstige elementen van twee arrays van dezelfde vorm toepast. Wanneer men een bewerking op alle elementen van een array toepast, dan spreekt men ook wel van een \array-operatie".

13.1 Algebrasche array-operaties De bewerkingen optellen en aftrekken worden automatische als arraybewerkingen opgevat. Na de toekenningsopdrachten >> a = [ 1 , 2, 3]; b = [ 6 , 5, 4];

levert het optellen van de arrays

13. ARRAY-OPERATIES IN MATLAB

21

>> a+b ans = 7

7

7

en het aftrekken van de arrays >> a-b ans = -5

-3

-1

De arraybewerkingen vermenigvuldigen (*), delen (/) en machtsverheen (^) geeft men aan door voor het bewerkingsteken een punt ( . ) te plaatsen. Dus >> a.*b ans = 6

10

12

>> a./b ans = 0.1667

0.4000

0.7500

>> a.^b ans = 1

32

81

Bij al deze array-operaties gelden de volgende conventies: (1) Als bij een array-operatie een van de arrays gelijk is aan een getal, dan wordt dit getal gezien als een array waarvan ieder element gelijk is aan dat getal en waarvan de afmetingen gelijk zijn aan de afmetingen van het andere array. (2) Als ieder element van een array met hetzelfde getal vermenigvuldigd moet worden dan hoeft men geen punt voor het vermenigvuldigingsteken te zetten. Voorbeelden: >> 2.^b

22


ans = 64

32

16

>> [2 2 2].^b ans = 64

32

16

>> a./2 ans = 0.5000

1.0000

1.5000

3*a ans = 3

6

9

Informatie over deze bewerkingen kan worden opgevraagd met het commando help arith.

Het maken van een array van functiewaarden

x3 . 1 + x2 We willen aan de variabele y de rijvector (f (0); f (0:2); : : : ; f (1:8); f (2)) toekennen. Eerst maken we een array met de argumenten (0; 0:2; : : : ; 1:8; 1).

Beschouw de functie f (x) =

>> x = 0:0.2:2 x = Columns 1 through 7 0

0.2000

0.4000

0.6000

1.8000

2.0000

0.8000

1.0000

Columns 8 through 11 1.4000

1.6000

Het maken van de rij van functiewaarden gaat met array-operaties! >> y = x.^3 ./ (1 + x.^2)

1.2000

23


y = Columns 1 through 7 0

0.0077

0.0552

0.1588

1.3755

1.6000

0.3122

0.5000

0.7082


1.1506

In feite wordt hier een rij van tellers componentsgewijs gedeeld door een rij van noemers.

Oefening 2.12 (a) Voer de opdracht t = x.^3 / (1 + x.^2) uit. U krijgt geen foutmelding omdat u iets heel anders heeft uitgerekend!

Opmerking: Ook al krijgt u uitvoer, wees altijd kritisch!

(b) Voer de opdracht t = x^3 ./ (1 + x.^2) uit. Wat betekent de foutmelding?

Oefening 2.13 Maak voor de volgende functies f het array (f (0); f (0:1); : : : ; f (2)). (a) f (x) = x2 + 2x + 1. 3x . 1 + 3x x p . (c) f (x) = 1+ x

(b) f (x) =

13.2 Relationele array-operaties Relationele array-operaties zijn operaties waarbij arrays elementsgewijs worden vergeleken. Bijvoorbeeld voor de matrices A=

3 5 20 11 9 1

en B =

4 3 11 8 2 4

is de opdracht om te kijken welke elementen van A groter dan die van B zijn: >> A > B ans = 0 1

1 1

1 0

24


Een 1 geeft aan welke elementen van A groter zijn dan de overeenkomstige van B . Het vergelijken van elementen kan alleen als de matrices A en B dezelfde afmetingen hebben of een van de matrices een getal is. Zo geeft >> A > 2 ans = 1 1

1 1

1 0

De relationele operaties staan in de volgende tabel. Operatie Betekenis > groter dan >= groter dan of gelijk aan < kleiner dan <= kleiner dan of gelijk aan == gelijk aan ~= niet gelijk aan Relationele operaties worden vaak gebruikt in combinatie met de functie find. Bekijk het resultaat van help find.

Oefening 2.14 Voer de opdracht

a = rand(5,5) uit. Deze opdracht genereert een array met 5 rijen en 5 kolommen waarbij ieder getal willekeurig gekozen wordt tussen 0 en 1. Maak een array met al die elementen van a, die kleiner dan 0.3 zijn. Hint: Gebruik bij deze opdracht find.

13.3 Andere array-operaties In MATLAB werken de wiskundige standaardfuncties sin, cos, sqrt, atan, asin etc. automatisch op array's. >> a = [1 2 3]; >> sin(a) ans = 0.8415

0.9093

0.1411

1.4142

1.7321

>> sqrt(a) ans = 1.0000


25

Beschouw de functie f (x) = x sin(x). Een rij van functiewaarden wordt met de volgende opdracht gegenereerd >> x = 0 : 0.25 : 4*pi; >> y = x .* sin(x)

Het resultaat is een array met 51 functiewaarden. t p . 1+ t Maak m.b.v. array-operaties en de functie sqrt de rij (f (0); f (0:1); f (0:2); : : : ; f (1)).

Oefening 2.15 Beschouw de functie f (t) =

Enkele andere operaties zijn >> size(A)

Het resultaat geeft de afmetingen van A aan, dat wil zeggen het aantal rijen en het aantal kolommen. >> sum(A)

Deze functie berekent de kolomsommen van een rechthoekig array A en geeft deze sommen weer in een rijvector. (Een kolomsom is de som van de elementen van een kolom.) Als A een rij is, dan wordt de som van de elementen van A berekent. De functie >> prod(A)

doet hetzelfde voor het product.

Oefening 2.16 Geef een eenregelige opdracht die bij een array a het aantal elementen van

a bepaalt die groter dan 5 zijn. Test uw opdracht uit met onderstaande arrays.

(a) a = [1, 2, 3; 4, 5, 6; 7, 8, 9] (b) a = 1:10 (b) a = 10*rand(6,6) Dezelfde opdracht moet in alle drie de gevallen werken!

Oefening 2.17 Geef een eenregelige opdracht die bij een array a de gemiddelde waarde van de elementen van a berekent. Test uw opdracht uit met onderstaande arrays.

26


(a) a = [1, 2, 3; 4, 5, 6; 7, 8, 9] (b) a = 1:10 (b) a = rand(6,6) Dezelfde opdracht moet in alle drie de gevallen werken!

13.4 Transponeren In een array (matrix) kunnen de rijen en kolommen verwisseld worden. Dit heet transponeren. Het resultaat wordt de getransponeerde array genoemd. Het transponeren gebeurt in MATLAB met de functie transpose. Na de opdrachten >> a = [1 2 3] a = 1

2

3

>> A = [1 2;3 4;5 6] A = 1 3 5

2 4 6

gaat het transponeren als volgt: >> transpose(a) ans = 1 2 3

en >> transpose(A) ans = 1 2

3 4

5 6

27

14. GRAFIEKEN

Het transponeren kan ook met behulp van de combinatie van een punt en een accent ( .' ). De opdrachten transpose(A) en A.' leveren hetzelfde resultaat op.

14 Gra eken Gra eken worden getekend met de `plot'-opdracht. Het commando >> plot(v)

waarbij v een rijvector is die de getallen v1 t/m vN bevat, genereert een plaatje waarin de punten (1; v1 ), (2; v2 ) t/m (N; vN ) achtereenvolgens met lijnstukken verbonden zijn. Het commando >> plot(v,w)

met v en w twee rijvectoren met evenveel elementen, tekent de punten (vi ; wi ) en verbindt deze achtereenvolgens met lijnstukken. Het commando >> plot(v,w,x,y)

tekent in dezelfde guur de punten (vi ; wi ) en (xi ; yi ). sin(x) + 2x Veronderstel dat we de gra ek van de functie f (x) = willen tekenen. Het tekenen 1 + x2 van de gra ek van f (x) met MATLAB gebeurt door een groot aantal punten (xi ; f (xi )) achtereenvolgens met rechte lijnstukken te verbinden. De vloeiendheid van de gra ek wordt natuurlijk door het aantal punten bepaald. Het plotten van de gra ek kan met de volgende opdrachten: >> x=0:pi/100:2*pi; >> y=(sin(x)+2*x)./(1+x.^2); >> plot(x,y)

Bij het tekenen van gra eken worden de assen automatisch gekozen. Men kan zelf de lengte van de assen kiezen met >> axis([XMIN XMAX YMIN YMAX])

waarbij de gra ek getekend wordt in de rechthoek XMIN YMAX.

x XMAX en YMIN y

Na iedere nieuwe `plot'-opdracht verdwijnt de informatie van de oude `plot'. Wil men verschillende gra eken met verschillende `plot'-opdrachten in dezelfde guur laten tekenen, dan moet men na de eerste `plot'-opdracht het commando

28


>> hold on

geven. MATLAB onthoudt dan de gra eken. Het commando >> hold off

herstelt de toestand waarin iedere volgende `plot'-opdracht de voorgaande plot doet verdwijnen. Voor het omgaan met het gra sche scherm zijn ook de volgende commando's van belang: Het tonen van een gra ek gaat door activeren of met de opdracht >> shg

Het gra sche scherm blijft bewaard tot een nieuw `plot'-commando wordt gegeven of tot het wordt leeggemaakt. De opdracht voor het leegmaken is >> clf

Voor het plaatsen van tekst en dergelijke in of bij een tekening zijn de volgende commando's van belang: >> grid

voorziet de gra ek van een raster. >> title('titel')

zet een titel boven een `plot', en wel de string aangegeven tussen de aanhalingstekens, in bovenstaand voorbeeld het woord \titel". >> text(x,y,'string')

zet in het punt (x; y) van de laatste `plot' de tekst \string". >> xlabel('tekst')

schrijft informatie langs de x-as. >> ylabel('tekst')

doet hetzelfde voor de y-as. Het opslaan van een guur in een le gebeurt met het commando

29

15. SCRIPT FILES

>> print -depsc filename

MATLAB maakt een le aan volgens het opgegeven `format', hier aangegeven met de optie -depsc ( Level 1 color Encapsulated PostScript (EPS) ). De le wordt weggeschreven naar de le \ lename.eps". Figuren in EPS-formaat kunnen in documenten van allerlei soort worden toegevoegd.

Oefening 2.18 Teken de functies sin(t) en cos(t) op het interval 0 t 2 in een guur.

p

Oefening 2.19 Teken de functie t sin(2t) op het interval [0; ]. Zet bij deze gra ek de

namen van de assen.

Oefening 2.20 Teken de functie e

t2

op het interval [ 2; 3].

15 Script les Als het maken van een plaatje uit meerdere opdrachten bestaat, dan kan men beter de opdrachten bij elkaar in een le zetten en deze door MATLAB in een keer laten uitvoeren. De le moet een naam met achtervoegsel `.m' hebben. Zo'n le waarin opdrachten staan die moeten worden uitgevoerd heet een `script le'. Het tekenen van de functie f (x) = jx sin(2x)j op het interval [0; 2] met een `script le': Via het menu \FilenNewnM- le" wordt de \MATLAB Editor/Debugger" geopend. Tik de volgende regels in. x = 0:0.1:2; y = abs(x .* sin(2*pi*x)); plot(x,y) title('f(x) = |x sin(2 pi x)| op interval [0,2]') axis([0 2 0 2]) shg

Sla de le op via het menu \FilenSave" onder de naam \plaat.m". Als u in MATLAB de opdracht plaat geeft, dan zal de gra ek van f op het interval [0; 2] verschijnen. Het is beter om de opdrachten, die een guur genereren, te bewaren dan de guur zelf.

Oefening 2.21 Waarom kan men beter de opdrachten voor een guur dan de guur zelf bewaren?

Oefening 2.22 Beschrijf de betekenis van de regels in de le \plaat.m".

30


16 Eigen functies In MATLAB kan men ook eigen functies de nieren. MATLAB gaat er standaard van uit, dat alle functies op arrays werken. Men moet bij het maken van een eigen functie rekening met array-operaties houden. Men kan dan de eigen functies met MATLABfuncties combineren. Een functiede nitie moet worden opgeslagen in een `fuction le', een le met de extensie `.m'. De naam van de le moet gelijk zijn aan de naam van de functie met de extensie `.m'. Beschouw de functie f (x) = x2 + ex . In MATLAB maken we een functie met dezelfde naam. Deze de nitie moet worden opgeslagen in de le \f.m". Via het menu `FilenNewnM- le" wordt de \MATLAB Editor/Debugger" geopend. Als men de volgende twee regels intikt function y = f(x) y = x.^2 + exp(x);

en de le onder de naam \f.m" opslaat, dan heeft men binnen MATLAB de beschikking over de functie f . We maken enkele opmerkingen over bovenstaande le.

De eerste regel van de le moet het woord \function" bevatten. De gebruikte variabelen zijn locaal en MATLAB zelf hoeft ze niet te kennen.

Als x een array is, dan wordt y een array van functiewaarden. De puntkomma op het eind verhindert dat bij iedere functie-evaluatie er onnodige uitvoer op het scherm komt.

Test altijd de functie f uit. Oefening 2.23 Voer de de nitie van f zelf uit. Waar komt de le \f.m" te staan?

U kunt de de nitie van f veranderen m.b.v. het commando edit f.

Oefening 2.24 Laat in de de nitie van de functie f de puntkomma weg. Bereken enkele waarden. Wat is het eect?

x ; x0 0 ; x<0 Maak een eigen functie g in MATLAB. Gebruik hierbij de MATLABfunctie max. Controleer uw eigen functie met

Oefening 2.25 Beschouw de functie g(x) =

>> x = -1:0.2:1; >> g(x)

31

17. DE HELPFACILITEITEN

ans = Columns 1 through 7 0

0

0

0

0.8000

1.0000

0

0

0.2000


0.6000

Waarschuwing: De namen van functies en variabelen mogen niet hetzelfde zijn.

Oefening 2.26 Maak een eigen functie met de naam h bij de functie h(x) = x2 die bestand is tegen arrray-operaties. Voer de volgende opdrachten in de juiste volgorde uit. (a) >> h([1,2,3,4]) (b) >> h = 2.5 (c) >> h([1,2,3,4]) (d) >> h(1) (e) >> h(1.4) (f) >> clear h Wordt bij de onderdelen (c) t/m (e) h als een functie of als een variabele gezien? Waarom volgt er bij onderdeel (d) geen foutmelding? Kent MATLAB de functie h nog?

17 De helpfaciliteiten Dit is een van de belangrijkste paragrafen. In principe is alle informatie over MATLAB ook te vinden in de helpfaciliteiten van MATLAB. Deze handleiding is tot nu toe zo uitgebreid geweest om u een snelle start in MATLAB te geven. Omdat de commando's en functies in MATLAB vaak in meer algemene situaties gebruikt kunnen worden, is de informatie algemener dan voor u nodig is. U zult de voor u essentiele bestanddelen uit de aangeboden informatie moeten halen. We behandelen drie manieren om informatie te verkrijgen.

Directe informatie

Met het commando

32


>> help

of

>> help ``trefwoord''

Zonder speci catie krijgt u een algemene lijst van commando's, functies en `toolboxes' die binnen MATLAB ter beschikking staan. Met speci catie krijgt u gerichte informatie en worden ook vaak suggesties gedaan om bij andere trefwoorden informatie op te vragen.

Oefening 2.27 Vraag informatie op over de sinus en de functie load.

Informatie via het \Help Window"

Via aanklikken van `HelpnHelp Window' of door de opdracht >> helpwin

Dan verschijnt het gevraagde `window' en kan men via klikken informatie opvragen.

Oefening 2.28 Open het \Help Window" en dubbelklik op de regel met \matlabnelfun". Er verschijnt een lijst met alle elementaire wiskundige functies in MATLAB. Dubbelklik op de regel met \atan" en bekijk de informatie. Probeer ook de mogelijkheden van het \See also"-menu uit. De mogelijkheid \More sin help(HTML)" hangt samen met de derde manier van informatie opvragen.

Informatie via de `hypertext'-documentatie

Via aanklikken van `HelpnHelp Desk(HTML)' of door de opdracht >> helpdesk

Dan wordt via een `Web browser' de \Help Desk" opgestart. Via klikken en menu's kan men nu gericht naar informatie zoeken. Het zoeken van de informatie kan alfabetisch of per onderwerp.

Oefening 2.29 Zoek via de \Help Desk" informatie over de sinus op.

Oefening 2.30 Klik in \Help Desk" achtereenvolgens `Getting Started' en `Command Line Editing' aan. Probeer zelf de pijltjestoetsen en de Èscape'-toets uit.

Via de \Help Desk" is het ook mogelijk de oÆciele MATLABdocumentatie te bekijken.

Oefening 2.31 Probeer de handleiding \Getting Started with MATLAB" te bekijken. Klik hiertoe `Getting Started' aan in de \Help Desk" en vervolgens het PDF-icoon van de genoemde handleiding aan. Als alles goed gaat, verschijnt er een `Window' met rechts een titelpagina en links een strook met driehoekjes en rechthoeken. Klik er een paar aan en zie wat er gebeurt.

Als het bovenstaande niet lukt, sla deze oefening dan verder over.

In de loop van de tijd zult u de helpfaciliteiten steeds gemakkelijker gaan gebruiken. De \Help Desk" geeft het gemakkelijkst toegang tot allerlei informatie.

18. OPGAVEN

33

18 Opgaven Oefening 2.32 Maak in MATLAB een rijvector met de vijfde machten van de eerste 40 natuurlijke getallen.

Oefening 2.33 Bereken in MATLAB het exacte product 1 2 : : : 24 25.

Oefening 2.34 Bereken in MATLAB de exacte som 1 + 23 + 33 + 43 + : : : + 1993 + 2003 . Oefening 2.35 Teken de gra ek van de functie f (x) = x=(1 + x) op het interval [0; 4].

Oefening 2.36 Teken de cirkel met straal 1 en middelpunt (0; 0) in het vierkant met 2 x 2 en 2 y 2. Zorg ervoor dat de cirkel er rond uitziet! Hints: Maak eerst een lijst van x-coordinaten en een lijst van y-coordinaten. Teken de cirkel. Pas de vorm van de guur aan door het commando axis in verschillende combinaties te gebruiken. Oefening 2.37 De periodieke functie f met periode 2 is gegeven door 8 <

t2 ; 0t2; f (t) = f (t 2) ; 2 t ; : f (t + 2) ; t < 0 :

De nieer deze functie in MATLAB. Maak gebruik van mod. Voor een getal r geeft mod(r,2) een getal s terug, zodanig dat 0 s < 2 en s en r een veelvoud van 2 van elkaar verschillen. Dus f (r) = f (s) = s2 . Zorg ervoor dat de eigen functie ook op arrays werkt.

34


Hoofdstuk 3

De numerieke gereedschapskist 1 Inleiding Het pakket MATLAB zelf is de numerieke `toolbox'. We maken enkele algemene opmerkingen en behandelen enkele belangrijke commando's. Bekijk van alle genoemde commando's de helpinformatie. Het aantal cijfers van getallen op het scherm kan veranderd worden met behulp van de commando's format long of format short. De berekeningen zelf veranderen niet. MATLAB rekent intern met 16 cijfers nauwkeurig. We gaan de functie f (x) = x3 x2 3 arctan(x) + 1 op het interval [ 2; 3] onderzoeken. Maak bij een te onderzoeken functie altijd een `M- le'.

Oefening 3.1 Maak een eigen functie f met f (x) = x3 x2 3 arctan(x) + 1.

2 De gra ek Het tekenen van de gra ek van f kan op de eerder beschreven manier met behulp van een array van functiewaarden. Een andere manier is gebaseerd op het gebruik van `strings'. >> ezplot('f(x)',[-2,3])

De expressie f(x) moet tussen enkele aanhalingstekens (') staan. Het interval wordt als een rijvector bestaande uit het linkereindpunt en rechtereindpunt gegeven. Het voordeel van deze opdracht is dat het interval gemakkelijk is te veranderen om de gra ek van f gedetailleerder te onderzoeken. 35

36

HOOFDSTUK 3. DE NUMERIEKE GEREEDSCHAPSKIST

3 De nulpunten Nulpunten van verreweg de meeste functies kunnen niet exact bepaald worden. Deze kunnen dan wel numeriek benaderd worden. Men spreekt in het algemeen toch van het numeriek bepalen van nulpunten. De functie f heeft drie nulpunten. Het eerste nulpunt ligt in het interval [ 2; 1] en de functiewaarden f ( 2) en f ( 1) zijn verschillend van teken. >> fzero('f',[-2,-1]) ans = -1.2780 >> fzero('f(x)',[-2,-1]) ans = -1.2780

Het is dus mogelijk om 'f' of 'f(x)' te gebruiken. In het laatste geval moet per se de variabele x gebruikt worden. >> f(ans) ans = -4.4409e-016

Aan de laatste uitvoer kunt u zien dat MATLAB een benadering voor het nulpunt geeft.

Oefening 3.2 Laat zien dat de benadering hooguit 10

4

van het nulpunt a igt.

Oefening 3.3 Bepaal de overige nulpunten van de functie f numeriek.

4 De extrema In het algemeen kunnen ook de extrema niet exact bepaald worden. In principe kunnen (de plaatsen van) de extrema numeriek worden benaderd. Men spreekt ook wel van het numeriek bepalen van de extrema. De plaatsen waar de extrema van f ongeveer liggen, kunnen uit de gra ek worden afgelezen. Er is een minimum in het interval [0; 2]. >> x1 = fmin('f',0,2)

4. DE EXTREMA

37

x1 = 1.0878 >> x1 = fmin('f(x)',0,2) x1 = 1.0878 >> f(x1) ans = -1.3784

De functie f heeft een minimum in x = 1:0878 met waarde 1:3784 . Het is wederom mogelijk om 'f' of 'f(x)' te gebruiken. In het laatste geval moet per se de variabele x gebruikt worden. Het interval, waarin het minimum ligt, moet als een rijvector worden ingevoerd. De functie f heeft een maximum op het interval [ 1; 0]. Helaas heeft MATLAB geen commando fmax.

Oefening 3.4 Teken de gra ek van de functie x 7! f (x) op het interval [ 2; 3]. Verander de le \f.m" niet!

De volgende opdrachten zullen nu duidelijk zijn. >> x2 = fmin('-f(x)',-1,0) x2 = -0.5902 >> f(x2) ans = 2.0456

De functie f heeft in x = 0:5902 een maximum ter waarde 2:0456 . Door fmin wordt de plaats van een minimum in het opgegeven interval benaderd. De waarde van het minimum moet u zelf uitrekenen. De nauwkeurigheden waarmee de plaatsen van de extrema worden uitgerekend, kunnen worden aangepast.

38


Oefening 3.5 Zorg ervoor dat de getallen met 15 cijfers op het scherm verschijnen. Voer de opdracht

>> fmin('f(x)',0,2,[0,10^n])

vier keer uit, waarbij u achtereenvolgens de exponent n vervangt door 4; 6; 8 en 10. Bekijk de gelijkblijvende cijfers van de verschillende uitkomsten.

Oefening 3.6 De functie f (x) = (x3

2)2 heeft precies een nulpunt. Bepaal het nulpunt numeriek en vergelijk de gevonden waarde met de werkelijke waarde. Hint: Teken de gra ek van f over het interval [0; 2].

Oefening 3.7 Van de functie f (x) = e2x 4ex +2 moeten alle nulpunten en extrema bepaald worden. Kies een geschikt interval om de functie met MATLAB te tekenen. Het interval is geschikt als alle nulpunten en extrema erin liggen en duidelijk te zien is waar deze liggen. Leg uit waarom het een geschikt interval is. Bepaal de nulpunten en extrema (aard en waarde).

5 Numerieke integratie Integralen van de vorm De integraal 115 12

Z

3 2

(x3

Z b a

x2

f (x)dx kunnen meestal niet exact worden bepaald.

3 arctan(x) + 1)dx is wel exact te berekenen. De uitkomst is

9 arctan(3) + 32 log(2) + 6 arctan(2) 6:0245344593535 .

Er zijn verschillende commando's om de integraal numeriek te bepalen (d.w.z. te benaderen). Zie ook de helpfaciliteiten. We gebruiken het commando quad en laten de getallen in 15 cijfers op het scherm verschijnen. >> quad('f',-2,3) ans = 6.02456548407172

Helaas is het niet mogelijk om 'f(x)' als eerste argument te gebruiken. De benadering wijkt al in het vijfde cijfer af. Ìn principe' probeert quad vier correcte cijfers te geven. De nauwkeurigheid van de benadering kan worden opgeschroefd. >> quad('f',-2,3,[10^-7 10^-7]) ans = 6.02453446058438

39

6. OPGAVEN Z

1 )dx numeriek. x2 1 Experimenteer met nauwkeurigheden. Vergelijk de resultaten met 883 , de exacte waarde van de integraal.

Oefening 3.8 Bepaal met behulp van quad de integraal

3

(x3 + x2 +

Zorg ervoor dat de getallen weer met vijf a zes cijfers op het scherm verschijnen.

6 Opgaven Oefening 3.9 Beschouw de functie f (x) = x2 exp( x).

Beantwoord de onderdelen (a) t/m (c) zonder gebruik te maken van de afgeleide f 0. (a) Waarom heeft de functie f een minimum in x = 0? (b) Waarom heeft f een maximum op het interval (0; 1)?. (c) Bepaal met fmin een maximum op (0; 1). (d) Hoever ligt de in onderdeel (c) gevonden benadering af van de werkelijke waarde?

Oefening 3.10 Beschouw de functie f (x) = x sin2 (x). (a) Teken de functie f op het interval [0; 2]. (b) Bepaal numeriek de plaats en de waarde van de maxima van f op het interval [0; 2]. (c) Bepaal numeriek de integraal

Z 2 0

f (x) dx.

40


Hoofdstuk 4

De symbolische gereedschapskist 1 Inleiding Dank zij de \(Extended) Symbolic Math Toolbox" is het mogelijk om binnen MATLAB aan formulemanipulatie te doen. In toekomstige versies van MATLAB zal deze `toolbox' verder in MATLAB worden ingebouwd en er zullen nog de nodige veranderingen optreden. Voer de volgende opdracht uit: >> syms x

Als er geen foutmeldingen verschijnen, dan heeft u de beschikking over de \(Extended) Symbolic Math Toolbox". De betekenis van deze opdracht wordt later uitgelegd. Met deze `toolbox' wordt het aantal functies uitgebreid. Als u in de \Help Desk" aan de rechterkant `Symbolic Math Toolbox Ref' aanklikt, dan verschijnt er een lijst van functies die of nieuw zijn of aangepast zijn vanwege deze `toolbox'. Voor directe informatie is het help commando geschikt. >> help atan ATAN Inverse tangent. ATAN(X) is the arctangent of the elements of X. See also ATAN2. Overloaded methods help sym/atan.m >> help sym/atan ATAN

Symbolic inverse tangent.

41

42

HOOFDSTUK 4. DE SYMBOLISCHE GEREEDSCHAPSKIST

Oefening 4.1 Verwijder alle variabelen uit de \Work Space".

2 Expressies met variabelen In MATLAB kan men met variabelen, die een waarde hebben, werken. Als men in een expressie een variabele gebruikt die geen waarde heeft, dan krijgt men een foutmelding. >> clear x >> exp(-x)+sin(x) ??? Undefined function or variable 'x'.

In MATLAB kunnen `strings' worden ingevoerd, zoals 'exp(-x) + sin(x)' of 'x^2 + 3 x'. Het zijn expressies, aan beide kanten voorzien van het '-teken. In feite hebben we dit gedaan bij ezplot, fzero en fmin. Beschouw het resultaat van >> ezplot('exp(-x) + sin(x)',[0,2*pi])

De functie f heeft in het interval [3; 4] een nulpunt en in het interval [4; 5] een minimum. >> fzero('exp(-x)+sin(x)',[3,4]) Zero found in the interval: [3, 4]. ans = 3.1831 >> fmin('exp(-x)+sin(x)',4,5) ans = 4.7213

Het is zonder meer mogelijk om aan een variabele een string toe te kennen. >> clear x >> y = 'exp(-x) + sin(x)' y = exp(-x) + sin(x)

De variabele y is nu gentroduceerd, omdat deze een `string' als waarde gekregen heeft, maar de variabele x is in MATLAB onbekend. Het gebruik van x aan de rechterkant van het is-teken geeft in dit geval geen foutmeldingen. Merk op dat aan de uitvoer niet te zien is of

3. SUBSTITUTIE

43

het een expressie met symbolen of een `string' betreft! Men kan MATLAB met symbolen laten werken als men de beschikking heeft over de \(Extended) Symbolic Math Toolbox" en variabelen tot symbool verklaart. Dit gaat als volgt. >> x=sym('x') x = x

of >> syms x

Het declareren van meerdere variabelen of symbolen tegelijkertijd: >> syms x y z

Men kan nu met deze variabelen naar hartelust manipuleren. >> y = x^3 + x^2 + 1 y = x^3+x^2+1 >> v=sin(x) v = sin(x) >> y*v ans = (x^3+x^2+1)*sin(x)

3 Substitutie >> syms a b >> y = 3*sin(a) + cos(b)

44


y = 3*sin(a)+cos(b) >> y =subs(y,a,2) y = 3*sin(2)+cos(b) >> subs(y,b,5) ans = 3.0116

4 Dierentieren en integreren We illustreren de opdrachten aan de hand van het voorbeeld: >> syms x y >> y = atan(x)

Een keer en drie keer dierentieren: y = atan(x) >> diff(y,x) ans = 1/(1+x^2) >> diff(y,x,3) ans = 8/(1+x^2)^3*x^2-2/(1+x^2)^2

De onbepaalde integraal >> int(y,x)

Z

arctan(x) dx:

5. NUMERIEKE WAARDEN

45

ans = x*atan(x)-1/2*log(x^2+1)

De bepaalde integraal

Z

7

2

arctan(x) dx:

>> int(y,x,2,7) ans = 7*atan(7)-1/2*log(2)-1/2*log(5)-2*atan(2)

Dit is de exacte waarde van de integraal. De numerieke waarde is dan >> double(ans) ans = 6.6367

De opdracht double maakt van een exact getal een decimaal getal. Verreweg de meeste integralen kunnen niet worden uitgerekend. Beschouw het voorbeeld Z 4

1

p

arctan(x) 1 + x3 dx.

>> int(exp(-x) * sqrt(1 + x^3), x, 1, 4) ans = int(exp(-x)*(1+x^3)^(1/2),x = 1 .. 4) >> double(ans) ans = 1.0017

De MATLABfunctie double zorgt voor een numerieke benadering van de integraal.

5 Numerieke waarden De functie sym zet numerieke getallen om naar exacte getallen en double doet het omgekeerde.

46


MATLAB

Resultaat

sym(sqrt(2)) sqrt(2) sym(2.3654) 11827/5000 double(sin(2)+log(3)) 2.0079

Oefening 4.2 Bereken in MATLAB de integraal Wat is de numerieke waarde van deze integraal?

Z

3 1

(cos(x

1 1) + ) dx. x

De manier om een idee te krijgen van de grootte van een exacte uitdrukking is het toepassen van het commando double.

6 Het manipuleren van expressies Voer eerst de opdracht >> syms a b c d x y

uit en experimenteer met de onderstaande opdrachten. MATLAB

Betekenis van opdracht expand((a+b)^3) Werk de haakjes weg. factor(a^3+3*a^2*b+3*a*b^2+b^3) Ontbind in factoren. [t n] = numden(a/b+c/d) Breng onder een noemer. Dan wordt t de teller en n de noemer. simplify((x^2+2*x+1)/(x+1)) Vereenvoudig de uitdrukking. Het schrijven van de uitdrukking x=(1 + x=(1 + x)) als teller gedeeld door noemer gaat als volgt: >> syms x >> y = x/(1+x/(1+x)); >> [t,n]=numden(y) t = x*(1+x) n = 1+2*x

Teller en noemer zijn niet verder te vereenvoudigen.

7. SYMBOLEN, STRINGS EN GETALLEN

47

7 Symbolen, strings en getallen De functie f (x) = 2xe x sin(x) wordt met MATLAB onderzocht op het interval [0; ]. Voer de volgende commando's uit. >> syms x >> y = 2*x*exp(-x)-sin(x)

Nu zijn de variabelen x en y in MATLAB als symbolen gede nieerd. Het tekenen van de functie f (x) in MATLAB gaat met het commando >> ezplot(y,[0,pi])

Aan het plaatje is meteen te zien dat f een maximum, een minimum en twee nulpunten heeft in het interval [0; ]. De plaatsen van de extrema en de nulpunten zijn niet exact te bepalen. Deze punten moeten numeriek benaderd worden. De MATLABfuncties fmin en fzero werken echter met strings, dat wil zeggen uitdrukkingen die tussen enkele aanhalingstekens (') staan. De strings mogen expressies in de variabele x zijn. Van een symbolische uitdrukking kan eenvoudig een string gemaakt worden met de functie char. >> fzero(char(y),[0.5,1]) Zero found in the interval: [0.5, 1]. ans = 0.8030 >> fmin(char(y),1.5,2) ans = 1.8409 >> fzero(char(y),[2.5,3]) Zero found in the interval: [2.5, 3]. ans = 2.7923

De plaats van het minimum en de nulpunten zijn gevonden. Ook de plaats van het maximum is gemakkelijk te bepalen. >> fmin(char(-y),0,0.5) ans =

48


0.3384

Oefening 4.3 Bepaal de extreme waarden van de functie f in het interval (0; ). De rand-

punten van het interval doen niet mee.

Het zal vaak voor komen dat u van een symbolische uitdrukking een string moet maken of omgekeerd. U kunt met class altijd achterhalen met wat voor expressie u bezig bent. >> class(y) ans = sym >> z = char(y); >> class(z) ans = char

Veronderstel dat u een array van functiewaarden van f wil maken zonder eerst een eigen functie in te voeren, bijvoorbeeld f ( 1); f ( 0:6); f ( 0:2); : : : ; f (1). Dat gaat als volgt. >> x = -1:0.4:1; >> z = vectorize(y) z = 2.*x.*exp(-x)-sin(x) >> eval(z) ans = -4.5951

-1.6219

Een controle is >> syms x >> subs(y,x,0.6) ans = 0.0939

-0.2899

0.1288

0.0939

-0.1057

49

8. OPGAVEN

De benvloeding van symbolen, strings en getallen onderling maakt het gebruik van MATLAB moeilijk. Door veel oefening krijgt u er op den duur toch gevoel voor. We geven een voor voorbeeld. >> syms x >> 2.3456 + x ans = 1466/625+x

De omschrijving van het numerieke waarde 2:3456 naar de exacte breuk 1466=625 doet vermoeden dat het laatste resultaat in zijn geheel een symbolische waarde geworden is, want in MATLAB zelf wordt iedere exacte breuk meteen omgezet naar een numerieke waarde. Een controle met class leert dat dit inderdaad zo is.

8 Opgaven Een wezenlijk onderdeel van het maken van de volgende opgaven is het op papier zetten van de uitwerkingen en de resultaten. Schrijf de essentiele stappen op en niet iedere toetsaanslag. Symboliseer eventueel de benodigde variabelen met syms.

Oefening 4.4 Teken de gra eken van de volgende functies op het interval [0; 2]. (a)

f (x ) =

arctan(x) . 1 + x2

(b)

f (x) =

1 + sin(x) . (1 + cos(x))2

Oefening 4.5 Voer de volgende opdrachten letterlijk in en verklaar het verschil. (a)

x ^ 1=3 .

(b)

x ^ (1=3) .

Oefening 4.6 Substitueer met MATLAB in de uitdrukking x2 y3 + sin(xy) voor x en y

respectievelijk 2 en 3.

Oefening 4.7 Bepaal met MATLAB de onbepaalde integraal

Z

Laat met MATLAB zien dat dit antwoord correct is. Hint: Gebruik de MATLABfunctie numden.

Oefening 4.8 Bepaal onderstaande integralen met MATLAB. (a)

R 0

x=2 sin(x) dx.

R (b) 01 e

2x

cos(3x) dx.

1 dx. 1 + x4

50 (c) (d)

HOOFDSTUK 4. DE SYMBOLISCHE GEREEDSCHAPSKIST R1 0

R4

log10 x dx.

2)3 j dx. Vraag de numerieke waarde van het resultaat op. Bereken deze integraal door het integratie-interval in twee stukken te splitsen en de absolute waarden weg te werken. Zijn de resultaten met elkaar in overeenstemming? 0

j(x

2

Oefening 4.9 Bepaal met MATLAB de volgende twee integralen en leg uit wat het resultaat betekent: p R2 a. sin(t2 ) 1 + t3 dt. 0

b.

R

e

t2 dt.

Oefening 4.10 Teken met MATLAB de gra ek van f (x) = 3+sin(x) op het interval [ 1; 2]. Zorg ervoor dat de oorsprong (0; 0) zichtbaar is. Oefening 4.11 Bekijk het resultaat van de MATLABopdrachten >> syms x >> ezplot(asin(sqrt(x)),[-5,5])

Waarom is de gra ek alleen voor 0 x 1 getekend? Maak met MATLAB een betere tekening van de gra ek.

Oefening 4.12 Bereken na uitvoering van de opdracht >> syms x n

met MATLAB de integraal

Z

xn dx .

(a) Voor welke n is dit antwoord incorrect? (b) Geef aan n de waarde uit onderdeel (a) en bereken de integraal opnieuw. Is het antwoord nu correct?

Oefening 4.13 Verklaar het resultaat van de opdracht sqrt(x^2).

Oefening 4.14 Beschouw de functie f (x) = (log(1 + x) + x2 =2)=(1 + x)2 . Voer de volgende MATLABopdrachten uit

>> syms x y >> y = (log(1+x)+x^2/2)/(1+x)^2

Beantwoord de volgende vragen:

51

8. OPGAVEN

(a) Bereken f 0(x). (b) Bereken f 0(1). (c) Teken de gra eken van f (x) en f 0 (x) in een plaatje.

Oefening 4.15 Gegeven is de functie f (x) = e2x 5ex + 6. (a) Teken met MATLAB de gra ek van f op een geschikt interval. (b) Bepaal numeriek met fmin de plaats van het minimum van de functie f . (c) Bereken de functiewaarde en de waarde van de afgeleide in het door u gevonden punt.

(d) Hoeveel nulpunten heeft f ? p

Oefening 4.16 Beschouw de functie f (x) = x 2 + 8 3 (x 1)2 . (a) Teken de gra ek (eerst even uit de hand) van f . (b) Bepaal de nulpunten en de extrema van f . (c) Laat zien dat u alle nulpunten en extrema gevonden heeft.

Oefening 4.17 Beschouw de functie f (x) = x sin(x). (a) Bepaal de extrema van de functie f , die binnen het interval (0; 2) liggen. (b) Teken de gra eken van f (x) en sin(x) voor 0 x 2 in een plaatje. (c) Verklaar de ligging van de extrema van f ten opzichte van die van de sinus.

Oefening 4.18 Beschouw de vergelijking

px sin

1 1 = ; x>0: x 4

p

(a) Schets de functie x sin( x1 ) en de constante functie

1 4

in een plaatje.

(b) Waarom heeft de vergelijking voor x dicht bij 0 geen oplossingen? (c) Waarom heeft de vergelijking voor grote x geen oplossingen? (d) Hoeveel oplossingen zijn er? (e) Benader op z'n minst een oplossing.

Oefening 4.19 De vergelijking ex = ax(1 x) in de onbekende x heeft voor een positieve a

precies een oplossing tussen 0 en 1. Geef met behulp van plaatjes een schatting voor deze a. Opmerking: Het getal a is exact te berekenen.

52


Oefening 4.20 De functie f (x) = e

x2 (x +

px) heeft een maximum in de buurt van x = 0:6.

(a) Bepaal (een benadering van) de plaats x0 van het maximum. (b) Bereken f 0(x0 ).

Oefening 4.21 Bepaal met MATLAB f (10) (x), x 6= 1, van de functie f (x) = sin( Het antwoord past op een regel! Verklaar het resultaat.

x2 + 2x 3 ). x 1

Hoofdstuk 5

Lineaire Algebra 1 Inleiding MATLAB is een programmapakket voor het uitvoeren van matrixberekeningen. Een matrix is in dit verband niets anders dan een rechthoekig schema van getallen dat we ook kunnen opvatten als een array. De naam MATLAB is een afkorting van matrix laboratory. Er bestaan verschillende versies van MATLAB. Bij deze handleiding zal worden uitgegaan van MATLAB voor Windows versie zoals genstalleerd op uw notebook, inclusief de 'Extended Symbolic Math Toolbox'. Zoals uit de naam MATLAB blijkt, is de basisrekengrootheid een matrix. MATLAB kent voor het rekenen met matrices twee soorten bewerkingen; de matrixbewerkingen, zoals deze een rol spelen in de lineaire algebra, en de array-bewerkingen, waarmee u al heeft kennis gemaakt in de training MATLAB in het eerste trimester. De array-bewerkingen onderscheiden zich van de in de lineaire algebra gebruikelijke matrixbewerkingen doordat er een punt voor het bewerkingsteken wordt geplaatst. MATLAB kent ook matrixfuncties. Matrixfuncties zijn veelal opdrachten die, voor een gegeven matrix, binnen de lineaire algebra gebruikelijke berekeningen uitvoeren. Een aantal van deze functies zal binnen de cursus Lineaire Algebra worden behandeld.

2 Matrices 2.1 Wat is een matrix? Zij m; n 2 IN . Een m n matrix is een rechthoekig schema van getallen, gerangschikt in m rijen en n kolommen: 0 B B B @

a11 a12 : : : a1n a21 a22 : : : a2n .. .. .. .. . . . . am1 am2 : : : amn

1 C C C A

De getallen aij (1 i m, 1 j n) heten de elementen van de matrix. In deze notatie heet i de rij-index, j de kolom-index. We schrijven i.h.a. bovenstaande matrix kort op als de 53

54

HOOFDSTUK 5. LINEAIRE ALGEBRA

matrix A, of ook wel (aij ). MATLAB kent in wezen maar een rekengrootheid, namelijk een m n-matrix (een matrix met m rijen en n kolommen met m 1 en n 1). Zo wordt een scalar in MATLAB opgevat als een 1 1-matrix.

2.2 Het invoeren van matrices Een matrix kan op verschillende manieren ingevoerd worden, namelijk: 1. tussen [ en ], waarbij de rijen gescheiden worden door ; en de elementen gescheiden worden door spaties of komma's. Bijvoorbeeld >> A = [1 2 3;4 5 6;7 8 9]

of >> A = [1,2,3;4,5,6;7,8,9]

2. tussen [ en ], waarbij elke rij op een nieuwe regel begint, b.v.: >> A = [1 2 3 4 5 6 7 8 9]

3. via externe les, waarin de waarden van de variabelen zijn opgeslagen of waarin deze berekend worden. Een matrix element mag uit alle MATLAB uitdrukkingen bestaan. Zo is b.v.: >> a = [-1.3

sqrt(3)

(1+2+3) * 4/5]

een toegelaten opdracht. We kunnen ook matrices samenvoegen tot een nieuwe matrix. Zo geeft C=[A;a]

de matrix A met als vierde rij daaraan toegevoegd de matrix a. De formaten van de matrices moeten hiervoor wel bij elkaar aansluiten.

Oefening 5.1 Voer de volgende 4 bij 4 matrix in. 0 B

A = B @

1 7 3 0

2 2 2 1

1 0 1 4

3 1 1 8

1 C C A

55

2. MATRICES

Oefening 5.2 Voer de volgende 2 bij 3 matrix in.

B =

1 1:5 2 7 p 3 2 0:625

Oefening 5.3 Voer de volgende matrices in A=

1 2 3 4

; B=

5 6 7 8

0 B

; C=B @

1 1 1 1

2 2 2 2

1 C C A

Maak nu een matrix D met in de eerste twee kolommen de matrices A en B boven elkaar, en als 3e en 4e kolommen de matrx C .

2.3 Matrices met symbolische elementen Het commando `sym' construeert symbolische getallen, variabelen en objecten. Lees wat "help sym"je over dit commando te vertellen heeft. Met het commando `syms' kunnen meerdere symbolische objecten gelijkijdig geconstrueerd worden. De namen voor de objecten moeten in dit geval met een letter beginnen. Symbolische objecten kunnen als matrixelementen optreden. >> syms a b >> A=[ 1 2 a 3 b 4 1 2 3] A = [ 1, 2, a] [ 3, b, 4] [ 1, 2, 3]

Oefening 5.4 Voer de volgende 2 bij 3 matrix in. A =

a11 a12 a13 a21 a22 a23

2.4 Het invoeren van vectoren Een rijvector is in feite een matrix met een rij, en een kolomvector een matrix met een kolom. Vectoren worden dus net zo ingevoerd als matrices. Als het verschil tussen de opeenvolgende elementen van de rij steeds hetzelfde is, kan de rij ook gevormd worden door `variabele'=`beginwaarde' : `stapgrootte' : èindwaarde'

56


Hierbij is `beginwaarde' het eerste element van de rij, èindwaarde' het laatse element van de rij en `stapgrootte' het verschil tussen de opeenvolgende elementen. Bijvoorbeeld >> x=1:-.2:0

geeft x = 1.0000

0.8000

0.6000

0.4000

0.2000

0

Als het verschli tussen de opeenvolgende elementen van de rij een is kan de `stapgrootte' worden weggelaten. Bijvoorbeeld >> x=1:5

geeft als resultaat x = 1

2

3

4

5

Voor een kolomvector is het vaak het makkelijkst om eerst een rijvector te maken en deze vervolgens te transponeren (rijen als kolommen schrijven, zie x2.7). ' geeft in MATLAB de getransponeerde van een matrix aan. Bijvoorbeeld >> x=[1:5]'

geeft als resultaat de getransponeerde van de rij [1 2 3 4 5]. D.w.z. x = 1 2 3 4 5

Oefening 5.5 Maak de kolomvector met beginwaarde -1, eindwaarde 9 en stapgrootte 0.5. Oefening 5.6 Maak de kolomvector met beginwaarde 9, eindwaarde 0 en stapgrootte 1.

2. MATRICES

57

Oefening 5.7 Maak de kolomvector met beginwaarde a, eindwaarde a +20 en stapgrootte 2.

Merk op dat bij het transponeren het programma er van uit gaat dat de symbolen complexe waarden kunnen aannemen. (zie de opmerkingen bij het transponeren van matrices blz.58) Men kan dit oplossen door te transponeren met het commando x:0 (arraybewerking). (Een andere mogelijkheid is a als 'real' te declareren. Zoek op hoe dit moet met de help-functie (`help sym' of `help syms')).

2.5 Bijzondere matrices Behalve door direkte invoer, kunnen enkele bijzonder matrices ook rechtstreeks gevormd worden met MATLAB-opdrachten. Deze matrices zijn: { Een m n nulmatrix, d.w.z. een matrix waarin alle elementen de waarde nul hebben, wordt gevormd met de opdracht `zeros(m; n)'. { Een m n matrix, d.w.z. een matrix waarin alle elementen de waarde 1 hebben, wordt gevormd met de opdracht ònes(m; n)'. { Een mn matrix met willekeurige getallen wordt gevormd met de opdracht `rand(m; n)'.

{ Een n n diagonaalmatrix, d.w.z. een matrix met alleen op de diagonaal (de elementen met gelijke rijen kolomindex) elementen ongelijk aan nul, wordt gevormd met de opdracht `diag(v)', waarbij v een vector is met de n diagonaalelementen. { Een m n eenheidsmatrix, d.w.z. een diagonaalmatrix met enen op de diagonaal, wordt gevormd met de opdracht èye(m; n)' Een vierkante m n-matrix wordt gevormd als alleen m als index gebruikt wordt. Bijvoorbeeld een 3 3-eenheidsmatrix wordt gevormd met de opdracht: >> eye(3)

Wil men een matrix met dezelfde dimensies als een gegeven matrix A dan kan dit door gebruik te maken van het commando `size'. Het commando >> zeros(size(A))

vormt een nulmatrix met dezelfde dimensies als A. Van deze commando's accepteert alleen het commando `diag' een matrix als argument. Wanneer bij het commando `diag' een matrix A wordt gebruikt i.p.v. de vector v krijgt men als antwoord een kolomvector met de diagonaalelementen van A.

Oefening 5.8 Voer met behulp van bovenstaande commando's de volgende matrices in: a. De 6 bij 6 eenheidsmatrix. b. De 5 bij 10 nulmatrix.

58


c. De 5 bij 15 matrix waarin alle elementen de waarde 1 hebben. d. Een willekeurige 5 bij 5 matrix. e. Een diagonaalmatrix met op de diagonaal de getallen 1 t/m 5.

2.6 Indices Vaak is het gewenst om met onderdelen van matrices of met afzonderlijke elementen te werken. Ook dit is eenvoudig te realiseren in MATLAB. Elk element van een matrix is te bereiken met behulp van indices. De opdracht hiervoor heeft de algemene vorm A(m; n). Hierbij is A de matrix, waarvan het element met rij-index m en kolomindex n wordt gespeci ceerd. Als A een vector is, is een n index voldoende. Bij de matrix A uit oefening 5.1 levert A(3; 2) de waarde van dit element, namelijk 2. Hiermee is het o.a. mogelijk om een afzonderlijk element van de matrix te veranderen. De opdracht: >> A(3,2)=3

verandert bij de matrix A de oorspronkelijke waarde elementen van A blijven hetzelfde.

2 in de nieuwe waarde 3. Alle andere

Men kan in plaats een element van A ook meerdere elementen van A tegelijk aangeven. Dit doet men door in plaats van een losse index een vector van indices mee te geven. >> A(3,[2 4])

Geeft de elementen a32 en a34 uit A gerangschikt als in A, d.w.z. een 1 2 matrix. Het commando >> A([1 2 3],[2 4])

geeft dus de 3 2 matrix gevormd door de elementen van A in de eerste drie rijen en in de kolommen twee en vier. Gezien de opmerkingen uit x2.4 kunnen we dit laatse ook bewerkstelligen met het commando >> A(1:3,[2 4])

Wanneer alleen `:' wordt gebruikt, worden alle rijen of kolommen bedoeld. Dit kan in de volgende vormen voorkomen: { A(:; j )

is de j-de kolom van A.

{ A(i; :)

is de i-de rij van A.

{ A(:; :)

is hetzelfde als A.

59

2. MATRICES

Met name wanneer je met hele rijen of kolommen operaties uit wil voeren is deze mogelijkheid goed te gebruiken. Zo levert, bij de oorspronkelijk matrix A het commando: >> A(2,:)-7*A(1,:)

als resultaat de rijvektor op, die ontstaat door zeven maal de eerste rij van de tweede rij af te halen. Met het commando: >> A(2,:)=A(2,:)-7*A(1,:)

wordt hetzelfde resultaat berekend maar dit wordt nu toegekend aan de tweede rij van de matrix A. De matrix A is dus door dit commando veranderd.

Oefening 5.9 Voer bovenstaande commando's uit met de matrix A gegeven in oefening 5.1. Oefening 5.10 Ga na wat het eect is van het commando >> A([1 3],[1 3])=10*ones(2)

op de matrix A gegeven in oefening 5.1.

2.7 Matrixbewerkingen in MATLAB Binnen MATLAB is het mogelijk om rekenkundige bewerkingen en opdrachten uit te voeren met matrices als geheel. Hiertoe kent MATLAB een groot aantal commando's. Een aantal van deze commando's zullen in de loop van de cursus Lineaire Algebra (zie [2]) nog besproken en gebruikt worden. Meer over de wiskundige context vindt men in [3]. Voor dit moment zijn de volgende operaties binnen MATLAB van belang:

vermenigvuldigen:

Het vermenigvuldigen van twee matrices gebeurt op speciale wijze. Voor een m n matrix A = (aij ) en een n p matrix B = (bjk ) de ni eren we het product AB P van A met B als de m p matrix (cik ) waar cik = nj=1 aij bjk . Dus: 0 B B @

.. .. .. . . . ai1 ai2 : : : ain .. .. .. . . .

10

1

0 : : : b1k : : : .. . C B : : : b : : : CB 2 k C B CB = : : : c B C .. ik : : : A@ A @ . .. . : : : bnk : : :

1 C C A

Het product van de matrices A en B kan alleen gevormd worden als de formaten (dimensies) van A en B passend zijn: het aantal kolommen van A moet gelijk zijn aan het aantal rijen van B . Als AB bestaat hoeft BA niet te bestaan. Als zowel AB als BA bestaan geldt in het algemeen niet AB = BA. Het matrixproduct van A en B is binnen MATLAB altijd gede nieerd als A of B een getal is. Het vermenigvuldigen van een matrix met een getal (scalaire vermenigvuldiging) komt overeen met het vermenigvuldigen van alle matrixelementen met dat getal. Voor matrixvermenigvuldiging gebruikt men in MATLAB het symbool `*' dus

60


>> A*B

optellen en aftrekken: Het optellen >> A+B

en aftrekken >> A-B

van twee matrices A en B komt overeen met het optellen en aftrekken van de afzonderlijke elementen van de twee matrices. Deze bewerkingen zijn gede nieerd als de dimensies van de matrices gelijk zijn of als een van de twee matrices een getal is. In dat laatste geval wordt het getal bij elk element van de matrix opgeteld of afgetrokken.

machtsverheen: Het commando >> A^p

verheft de matrix A tot de p-de macht. Als p een positief geheel getal is, dan wordt A ^ p berekent door herhaald vermenigvuldigen van A met zichzelf.

transponeren: Zij A = (aij ) een m n matrix. dan is de getransponeerde AT van A de n m matrix (bij ) (1 i n, 1 j m) met bij = aji. Dus: 0

1T

0

1

a11 : : : a1n a11 : : : am1 B .. C B . .. C T .. A = @ ... A =@ . . A am1 : : : amn a1n : : : amn

In MATLAB is A0 de getransponeerde van de matrix A. Bijvoorbeeld: >> [1 2 3]'

vormt de kolommatrix: >> 1 2 3

Als: >> x=[1 2 3]

wordt dezelfde matrix verkregen met de opdracht: >> x'

Merk op dat indien een matrix A complexe, niet reele, elementen bevat het commando A0 de matrix niet alleen spiegelt maar ook van ieder element de complex geconjugeerde neemt. Wil men alleen de gespiegelde dan kan dit met het commando A:0 Ter illustratie een tweetal voorbeelden.

61

2. MATRICES

>> A=[2+3*i 4+5*i 2 3] A = 2.0000 + 3.0000i 2.0000

4.0000 + 5.0000i 3.0000

>> A' ans = 2.0000 - 3.0000i 4.0000 - 5.0000i

2.0000 3.0000

>> A.' ans = 2.0000 + 3.0000i 4.0000 + 5.0000i

2.0000 3.0000

Hiermee dient men bij het gebruik van symbolische matrixelementen rekening te houden. >> A=[a b 1 2] A = [ a, b] [ 1, 2] >> A' ans = [ conj(a), [ conj(b),

1] 2]

>> A.' ans = [ a, 1] [ b, 2]

De comando's `n' en `=': De `deling' is bij matrices niet gede nieerd. In MATLAB wordt echter bij matrices

62


wel de deelstreepnotatie gebruikt. We zullen spreken van een symbolische deling (binnen MATLAB ook aangeduid met `mldivide' en `mrdivide' ). De betekenis is alsvolgt >> X=A\B

geeft de oplossing X van de matrixvergelijking A X = B , en >> X=B/A

geeft de oplossing X van de matrixvergelijking X A = B . Omdat AX in het algemeen ongelijk is aan XA moet er onderscheid gemaakt worden tussen een symbolische deling aan de `rechterkant' en een symbolische deling aan de `linkerkant'. Wanneer de op te lossen vergelijking meerdere oplossingen heeft geeft MATLAB er slechts een. Wanneer de vergelijking geen oplossingen heeft geeft MATLAB de `kleinste kwadraten oplossing'. Het verdient dus aanbeveling het antwoord van MATLAB te controleren. We komen hierop terug in x3 Meerinformatie kunt u vinden in bijv. de \helpdesk"door onder \index"te kijken bij \arithmetic operators".

Oefening 5.11 Bepaal van de eerder ingevoerde matrix A (Oefening 5.1) de getransponeerde matrix AT d.m.v. het commando: B=A'. Bereken AAT via: C = A B . Ga na dat C symmetrisch is. Bereken AT A via: D = B A. Vergelijk het resultaat met C = A B . Bereken 3A + 5B 3 via: E = 3 A + 5 B ^ 3.

a. Voer de vector v = (1 2 3 4) in. Bereken vervolgens v A en A v0 , met A de matrix uit de voorgaande oefening. Wat gebeurt er door deze vermenigvuldigingen met de rijen en kolommen van A?

Oefening 5.12

b. Maak een diagonaalmatrix D met op de diagonaal de elementen van de vector v. Bereken D A en A D. Wat gebeurt er door deze vermenigvuldigingen met de rijen en kolommen van A?

Oefening 5.13 Bereken met de vector v uit de vorige oefening v v0 en v0 v. Ga na dat er sprake is van 'gewone' matrixvermenigvuldiging.

Oefening 5.14 Herhaal de opdrachten uit oefening 5.11 voor de matrix 0

A=@

1 b 3

2 a 3 2 1 c

1 A

:

2.8 Matrixfuncties in MATLAB In deze paragraaf noemen we een aantal commando's die een matrix als argument hebben (Denk b.v. ook aan `zeros(A)', blz. 55). Bij de training hoeft u de commando's in deze

3. OPLOSSEN VAN STELSELS LINEAIRE VERGELIJKINGEN

63

paragraaf niet uit te voeren. De meeste van deze commando's zullen worden gebruikt in het college Lineaire Algebra (zie [2]) waar ook duidelijk zal worden wat de genoemde begrippen precies betekenen. Op dit moment gaan we hier niet verder op in. Het commando >> inv(A)

berekent de inverse van een vierkante matrix A als A inverteerbaar is, d.w.z. berekent een matrix B zodat A B = B A = I . Als er geen inverse bestaat (A singulier) of als het numerieke resultaat niet betrouwbaar is zal MATLAB dit melden. ([2, week 2]). Het commando >> det(A)

berekent de determinant van A. A moet een vierkante matrix zijn. ([2, week 2]). Het commando >> rank(A)

geeft de rang van de matrix A. ([2, week 3]). Het commando >> eig(A)

geeft een vector met daarin de eigenwaarden van de A. De matrix A moet vierkant zijn. Dit wordt vaak gebruikt in de vorm >> [S,D]=eig(A)

Dit geeft als resultaat een matrix S en een matrix D, zodat A S = S D. Iedere kolom van S een eigenvector en D een diagonaalmatrix met op de diagonaal de eigenwaarden van A behorend bij de eigenvectoren in de overeeenkomstige kolommen van S . ([2, week 6]).

Opmerking: Bovenstaande commando's werken ook met symbolische matrices. I.h.a. moet

men er op bedacht zijn dat matrixfuncties mogelijkerwijs niet werken voor symbolische matrices. Omgekeerd zijn er ook commando's die alleen werken voor symbolische matrices.

3 Oplossen van stelsels lineaire vergelijkingen 3.1 Het commando A n b Binnen MATLAB is het mogelijk om rechtstreeks de oplossing van een stelsel vergelijkingen te bepalen. Het commando om Ax = b rechtstreeks op de lossen luidt:

64


>> x=A\b

(Binnen MATLAB kan men meer informatie vinden via het trefwoord "mldivide".) Soms geeft MATLAB, na de opdracht x = Anb, de volgende waarschuwing. >> Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = .... .

We gaan op de precieze betekenis van deze waarschuwing hier niet in. De waarde van x die op het scherm verschijnt, moeten we niet vertrouwen. Een andere mogelijke waarschuwing is >> Warning: Matrix is singular to working precision.

Dus: Verschijnt er een van bovenstaande waarschuwingen, dan is de waarde van x onbetrouwbaar. Indien deze meldingen niet worden gedaan, wil het nog niet zeggen dat A x = b. I.h.a. zal dit commando het stelsel oplossen met een Gauss-algoritme. Echter, dit commando geeft soms ook een antwoord x wanneer een stelsel vergelijkingen geen (exacte) oplossing heeft. In dat geval wordt een ander algoritme gebruikt. De betekenis van deze 'kleinstekwadratenoplossing' x komt in de cursus Lineaire algebra nog aan bod ([2, week 4]). Het is daarom belangrijk een antwoord dat je met dit commando bepaald hebt expliciet te controleren door A x uit te rekenen en te vergelijken met de gewenste uitkomst b.

Geldt niet A x = b, dan is het stelsel overbepaald. Geldt wel A x = b dan is het stelsel regulier of onderbepaald. MATLAB levert namelijk hoogstens een oplossing. Bij onderbepaalde stelsels, d.w.z. stelsels met meerdere oplossingen, moet de algemene oplossing dan anders bepaald worden, bijv. door vegen van de matrix. Bovenstaande methode werkt ook indien A en b symbolische matrices zijn. Men krijgt een waarschuwing als de oplossing niet bestaat.

3.2 De rijgereduceerde trapvorm MATLAB kent een commando om een matrix te vegen met het Gauss-Jordan algoritme, d.w.z. om de rijgereduceerde trapvorm van een matrix te bepalen. Binnen MATLAB kunnen we de rijgereduceerde trapvorm (reduced row echelon form) van een (symbolische) matrix A bepalen met het commando >> rref(A)

Passen we dit commando toe op de uitgebreide coeÆcientenmatrix van een stelsel vergelijkingen dan kunnen we uit het resultaat de oplossing bepalen.

3. OPLOSSEN VAN STELSELS LINEAIRE VERGELIJKINGEN

65

Oefening 5.15 Beschouw de matrix 0

A=@

1 2 3

2 3 1

3 2 1

1

0

A

en de vector b = @

6 14 2

1 A

en het stelsel Ax = b. Los dit stelsel op met het commando A\b en met het rref commando.

Oefening 5.16 Beschouw de matrix A=

1 2 2 1

3 3

en de vector b =

4 4

en het stelsel Ax = b. Los dit stelsel op met commando A\b. Is de oplossing uniek? Bepaal alle oplossingen door gebruik te maken van rref([A,b]). Maak symbolische matrices AA = sym(A) en bb = sym(b) en los het stelsel ook symbolisch op met het commando AA\bb.

Oefening 5.17 Beschouw de matrix 0

A=@

1 0 1 2 1 4

1

0

A

en de vector b = @

10 16 17

1 A

en het stelsel Ax = b. Ga met rref([A,b]) na dat dit stelsel strijdig is. Wat geeft het commando A\b? Is dit een oplossing? Maak symbolische matrices AA = sym(A) en bb = sym(b) en kijk wat het resultaat is van het commando AA\bb.

3.3 Het oplossen van stelsels met de Symbolic Toolbox De symbolic toolbox van MATLAB maakt het mogelijk ook vergelijkingen op te lossen waarin de coeÆcienten symbolische uitdrukkingen zijn. Het commando waarmee symbolische vergelijkingen (niet noodzakelijk lineair) kunnen worden opgelost is het commando 'solve'. Dit commando wordt alsvolgt gebruikt solve('expressie1','expressie2','onbekende1','onbekende2')

Hierin zijn expressie1 en expressie2 de vergelijkingen, en onbekende1 en onbekende2 de variabelen waarnaar men wil oplossen. Voorbeeld: >> syms a b c d p q >> [x1,x2]=solve('a*x1 +b*x2=p','c*x1 +d*x2=q','x1','x2') x1 =

66


(-b*q+p*d)/(-b*c+d*a) x2 = (-c*p+q*a)/(-b*c+d*a)

Meer detailinformatie over dit commando is te verkrijgen via de online-help van het programma.

Oefening 5.18 Beschouw het stelsel vergelijkingen uit oefening 5.16 en los dit stelsel op met behulp van het commando `solve'. Hoe ziet de algemene oplossing eruit in vectornotatie?

Oefening 5.19 Maak BK opgave 1.5.11 door de oplossing via symbolische matrices te bepalen in termen van de parameter a. Is de oplossing voor alle a gede nieerd? Beschouw de waarden van a waarvoor deze oplossing niet is gede nieerd apart, bijv. door gebruik te maken van het commando voor substitutie `subs'.

Oefening 5.20 Herhaal de opdracht uit oefening 5.19 voor BK opgave 1.5.13. Onderzoek p of er oplossingen zijn voor a = 6 via `rref' en via Ànb'.

4 Numerieke aspecten bij gebruik van MATLAB Wanneer we van een pakket als MATLAB gebruik maken dan moeten we er rekening mee houden dat het een numeriek pakket is. Bij numerieke algoritmen spelen de aspecten van eÆcientie, afrondfouten en datafouten een rol. We zullen ons beperken tot een korte illustratie. Deze onderwerpen komen in latere colleges in meer detail aan bod. Opmerking. Wanneer men binnen de symbolic toolbox numeriek rekent maakt men geen gebruik van de numerieke algoritmen van MATLAB maar van de algoritmen van Maple.

4.1 Getalrepresentatie en afrondfouten, een voorbeeld Alle waarden waarmee gerekend wordt worden gerepresenteerd met een beperkt aantal cijfers. De computer rond alle reele getallen af tot een getal dat past in de gebruikte representatie. Voorbeeld (standaard outputrepresentatie in MATLAB): 1 p3 = 0.3333 2 = 1.4142 e 10 = 4.5400e-005 e10 = 2.2026e+004 1 = 0.1000 10 1 = 1.0000e-006 106

4. NUMERIEKE ASPECTEN BIJ GEBRUIK VAN MATLAB

67

Intern rekent MATLAB met een grotere precisie dan op het scherm wordt weergegeven. Dat de outputrepresentatie tot onverwachte uitkomsten aanleiding kan geven wordt duidelijk in de volgende oefening.


9 0 0 1

Bereken met de hand A2 en A3 . Kunt u een formule geven voor An ? Bereken nu met MATLAB A10 . Is dit antwoord exact? Wanneer we een antwoord met meer cijfers willen representeren kan dit m.b.v. het commando >> format long

Geef het commando 'format long' en daarna opnieuw de opdracht A10 . vergelijk de uitkomst met het eerder verkregen antwoord. Geef vervolgens het commando >> format

om weer terug te keren naar een representatie met 5 cijfers. We zien dus dat 'afrondfouten' gevolgen kunnen hebben voor het resultaat. Gebruik 'help format' voor meer informatie over de mogelijke representaties van matlab output. De rol die afrondfouten in berekeningen spelen kan vaak beperkt worden door het slim inrichten van de algoritmen.

4.2 Datafouten, een voorbeeld Andere fouten waar je mee te maken kunt krijgen zijn 'datafouten'. In de volgende oefening zullen we illustreren dat, bij het oplossen van een stelsel Ax = b, een kleine verandering in de data grote gevolgen kan hebben voor de oplossing. Hierbij spelen intrinsieke eigenschappen van A een rol.

Oefening 5.22 (ML) Beschouw twee stelsels vergelijkingen, nl. 8 > > < > > :

en

8 > > < > > :

10x1 + 7x2 + 8x3 + 7x4 = 32 7x1 + 5x2 + 6x3 + 5x4 = 23 8x1 + 6x2 + 10x3 + 9x4 = 33 7x1 + 5x2 + 9x3 + 10x4 = 31 2x1 + x2 + 5x3 + x4 = x1 + x2 3x3 4x4 = 3x1 + 6x2 2x3 + x4 = 2x1 + 2x2 + 2x3 3x4 =

9 5 8 3

68


a. Bepaal van beide stelsels de oplossingen. Door bij deze stelsels de rechterleden iets te veranderen en te kijken naar het eect daarvan op de oplossing is het mogelijk om een indruk te krijgen hoe gevoelig de matrices zijn voor kleine veranderingen in de data. b. Los het eerste stelsel ook op met rechterleden (32:1; 22:9; 32:9; 31:1) en (32:01; 22:99; 32:99; 31:01). Los het tweede stelsel ook op met rechterleden (9:1; 5:1; 7:9; 3:1) en (9:01; 5:01; 7:99; 3:01). Ga na wat het eect van deze veranderingen op de oplossing is. Door de coeÆcienten van de vergelijkingen iets te veranderen en te kijken naar het eect daarvan op de oplossing is het ook mogelijk om een indruk te krijgen van de gevoeligheid voor kleine veranderingen in de data. c. Verander de elementen van de oorspronkelijke coeÆcienten{matrices door er de matrix 0:1 rand(4) bij op te tellen (in dat geval wordt bij elk element een willekeurig getal tussen 0 en 0.1 opgeteld). Bepaal nu de oplossingen van deze aangepaste stelsels bij de oorspronkelijke rechterleden, respectievelijk (32; 23; 33; 31) en (9; 5; 8; 3). En ga na wat het eect van deze veranderingen op de oplossing is. d. Bij stelsels die gevoelig zijn voor kleine veranderingen in de data noemen we de coeÆcientenmatrix slecht-geconditioneerd. Welk van beide stelsels vergelijkingen heeft de best geconditioneerde coeÆcientenmatrix?

Oefening 5.23 Met de MATLAB opdracht A = hilb(n) wordt A de n bij n Hilbert matrix, gede nieerd door aij = 1=(i + j

1).

a. Bepaal de 5 bij 5 Hilbert matrix. b. Los op Ax1 = b1 = (2:2833; 1:4500; 1:0929; 0:8845; 0:7456)T . c. Los op Ax2 = b2 = (2:2834; 1:4501; 1:0928; 0:8844; 0:7457)T . d. Vergelijk b1 met b2 en x1 met x2 . De matrix A is blijkbaar slecht geconditioneerd.

4.3 EÆcientie Het inrichten van algoritmen dient natuurlijk zo eÆcient mogelijk te gebeuren. Als maat voor de eÆcientie van een algoritme wordt wel gebruikt het aantal bewerkingen (optellingen, vermenigvuldigingen etc.) dat voor het uitvoeren van het algoritme nodig is. Bij MATLAB heet een enkele bewerking een ' op' ( oating point operation). Na a oop van iedere sessie geeft MATLAB het aantal gebruikte ops. De op-teller kan op nul gezet worden met het commando

5. MATLAB EN LINEAIRE ALGEBRA

69

>> flops(0)

Het aantal tot op dat moment gebruikte ops krijgt men met het commando >> flops

Oefening 5.24 Hoeveel ops zijn er nodig voor het vermenigvuldigen van twee n n-

matrices? Controleer dit door bijv. twee random 5 5- en 10 10-matrices met elkaar te vermenigvuldigen. Zet voor iedere vermenigvuldiging de teller op nulen vraag na de vermenigvuldiging het aantal ops op.

Oefening 5.25 Doe hetzelfde als in de voorgaande oefening maar nu voor het matrixproduct van een m n matrix met een n k-matrix.


A11

;

;

A22

met A11 en A22 willekeurige 40 bij 40 matrices en ; een 40 bij 40 nulmatrix. Het berekenen van A A kan op twee manieren, namelijk rechtstreeks en gebruik makend van het feit dat A11 A11 ; AA=A= ; A22 A22 In het laatste geval is vermenigvuldiging van een nulmatrix met een andere matrix niet nodig. Bepaal voor beide manieren, m.b.v. MATLAB, het aantal benodigde ops en vergelijk de resultaten

5 MATLAB en Lineaire Algebra We besluiten deze training met een aantal opdrachten die door toepassing van matrixrekening kunnen worden opgelost. MATLAB is hierbij het aangewezen gereedschap. Hoewel de opdrachten ook met de hand kunnen worden opgelost (eventueel ter controle) dient men bij iedere opdracht na te gaan hoe zoveel mogelijk van het rekenwerk op een eÆciente manier met MATLAB uitgevoerd kan worden. Vaak zal de oplossing verkregen worden door het oplossen van een stelsel vergelijkingen in dat geval moeten alle mogelijke oplossingen bepaald worden.

5.1 Opdrachten Oefening 5.27 Bepaal alle 2 2 matrices die commuteren met

1 0 0 3

:

70


Oefening 5.28 Bepaal alle 2 2 matrices die commuteren met

1 3 0 2

:

Oefening 5.29 Beschouw n deeltjes in het vlak met plaatsvectoren r1 ; r2 ; ; rn en massa's m1 ; m2 ; ; mn . Voor de plaatsvector rzp van het zwaartepunt van dit stelsel van n massa's geldt dan rzp =

1 (r m + r2 m2 + + rn mn ) ; M 1 1

met M = m1 + m2 + + mn . Beschouw een plaat in het vlak met hoekpunten (1,2), (5,1) en (4,4) waarvan de massa geheel geconcentreerd is in de hoekpunten. Hoe kunnen we een massa van 1 kg verdelen over de hoekpunten zodanig dat het zwaartepunt in het punt (3,3) ligt? Beschouw nu een plaat met vier hoekpunten waaronder in ieder geval de reeds gegeven drie punten (kies zelf een vierde punt). Wat zijn nu de mogelijke massaverdelingen?

Oefening 5.30 De impulsvector P van een stelsel van n deeltjes met massa's m1 ; m2 ; ; mn en snelheden v1 ; v2 ; ; vn is gegeven door P = m1 v1 + m2 v2 + + mn vn :

Beschouw twee deeltjes met snelheden v 1 = (2; 2; 2)T en v2 = (5; 7; 9)T . De deeltjes botsen en hebben na botsing de snelheden w1 = (4; 6; 4)T en w2 = (4; 5; 8)T . Wat kunnen we, onder de veronderstelling van behoud van impuls, zeggen over de massa's van de twee deeltjes?

Oefening 5.31 Teken in het gebied 6 x 6, 6 y 6 de driehoek met hoekpunten p1 =

1 2

; p2 =

5 1

; p3 =

4 3

Zij D de matrix D=

0 1

1 0

Teken de driehoek met hoekpunten bi = Dpi . Teken ook de driehoeken met hoekpunten D2 pi en D3 pi . Wat is het eect van de vermenigvuldiging met D? Met welke matrix moeten we vermenigvuldigen om te spiegelen in de x-as?

Hoofdstuk 6

Dierentiaalvergelijkingen 1 M- les 1.1 M- les maken Tot nu toe is MATLAB hoofdzakelijk gebruikt in een commando-gestuurde mode; de commando's worden regel voor regel ingevoerd, waarbij de commando's op een regel onmiddelijk worden uitgevoerd. MATLAB kan echter ook een reeks van commando's, die opgeslagen zijn in een le, uitvoeren. Dergelijke les heten M- les naar de extensie `.m'. Een M- le bestaat uit een serie opdrachten in de MATLAB taal. Andere programma's in die taal kunnen een M- le aanroepen. Een M- le kan ook zichzelf aanroepen en zo een recursief proces vormen. Een M- le is een gewone ASCII- le die geschreven kan worden met iedere gewone teksteditor. Het ligt echter het meest voor de hand om de in MATLAB ingebouwde faciliteit te gebruiken. Deze wordt opgeroepen via het menu 'File/New/M- le'. Hiermee wordt de 'MATLAB Editor/Debugger' geopend. Een M- le is een le met de extensie '.m'. De le bevat een reeks matlabopdrachten. Ook alle standaard matlabcommando's liggen vast in een le 'commando.m'. De opdrachten in de le commando.m worden uitgevoerd door binnen MATLAB de opdracht >> commando

te geven.

1.2 Script- en function les Er bestaan twee soorten M- les, de zogenaamde script les en de function les (zie ook x16). In de informatica is een script een commandoreeks die genterpreteerd wordt in plaats van gecompileerd. Een script le bestaat dus uit een aantal commandoregels die uitgevooerd 71

72

HOOFDSTUK 6. DIFFERENTIAALVERGELIJKINGEN

worden als commando's die in de commando-gestuurde mode worden ingevoerd. De variabelen kunnen dus reeds eerder gede nieerde variabelen zijn en nieuwe variabelen kunnen door de script le worden gentroducerd. De variabelen zijn zogenaamde globale variabelen. Bijvoorbeeld in de le "rijom.m", een le die de volgorde van de elementen in een rij omdraait, kunnen de volgende opdrachten voorkomen [r k]=size(x); y=x(k:-1:1)

Deze script le wordt aangeroepen met de opdracht `rijom'. Het programma neemt aan dat de variabele x bestaat. Na a oop van het programma staan de variabelen x, r, k en y in het geheugen. Een function le begint met een opdracht die de in- en uitvoer variabelen van de opdracht de nieert. Een voorbeeld van een function le, die hetzelfde doet als bovenbeschreven script le, is: function y=wissel(x) % De functionfile , % draait de volgorde van de elementen van de rij x om % en bergt het resultaat op in matrix y. [r k]=size(x); y=x(k:-1:1)

De eerste regel begint met `function' om aan te geven dat het niet om een `script' gaat. De naam van de functie is `wissel' en de naam van de le is wissel.m . De functie kent 1 parameter x. Deze variabelenaam is alleen van belang binnen de functie le en bestaat niet buiten de functie le. Het is een zogenaamde locale variabele. Het resultaat van de functie heet y. Ook y is een locale variabele. De variabelen die je in een functie le gebruikt zijn dus afgeschermd van de andere variabelen en benvloeden de werking van de variabelen buiten de functie niet. Het %-teken geeft aan dat er een commentaarregel volgt. Het aaneengesloten commentaarblok na de eerste regel van een function le bevat doorgaans een beschrijving van de functie. Het commando `help wissel' geeft de deze beschrijving op het scherm. De aanroep à=wissel(b)' zal nu de kolommen van een vooraf gede nieerde matrix b in volgorde omwisselen en het resultaat opbergen in matrix a. Er is bij een functie le sprake van locale variabelen. Wil men variabelen in een functie le ook buiten de functie le betekenis geven dan kan dit met het commando `global' (zie ook help global).

Oefening 6.1 Schrijf de bovenbeschreven script le rijom.m Maak het MATLAB geheugen leeg met het commando `clear' en de nieer de rij x=[1 3 5 7 9 11]

Voer vervolgens de script le uit. Wat is de verwachte uitvoer van deze script le ? Bekijk met het commando `who' welke variabelen in het geheugen staan.

2. DIFFERENTIAALVERGELIJKINGEN

73

Oefening 6.2 Schrijf de bovenbeschreven function le wissel.m. Maak het MATLAB geheugen leeg met het commando `clear' en de nieer de rij a=[1 3 5 7 9 11].

Pas de function le toe op a en bekijk daarna met het commando `who' welke variabelen in het geheugen staan.

2 Dierentiaalvergelijkingen 2.1 Het numeriek oplossen van dierentiaalvergelijkingen Voor het berekenen (eigenlijk 'benaderen') van de oplossing van de dierentiaalvergelijking wordt gebruik gemaakt van algoritmen die verwant zijn aan het Euler algoritme dat in het college Calculus is behandeld. Het Euler algoritme komt, in het kort, neer op het volgende. Hoewel men de oplossing van de dierentiaalvergelijking niet expliciet kent, kent men wel de afgeleide van de oplossing in ieder punt (deze wordt gegeven door de di.vgl.), m.a.w. men kent het richtingsveld. Dit richtingsveld wordt nu gebruikt om de oplossing te benaderen. Vanuit een gegeven beginwaarde wordt het volgende punt van de oplossing berekend door vanuit de beginwaarde een stap te zetten in de richting van het richtingsveld. Vanuit het zo verkregen punt doet men weer hetzelfde, etc. De benadering van de oplossing hangt dus af van de beginwaarde, het richtingsveld en de stapgrootte. In MATLAB kan men dierentiaalvergelijkingen numeriek oplossen met behulp van de commando's òde23' of òde45'. De onderliggende algoritmen maken gebruik van Runge-KuttaFehlberg integratie met variabele stapgrootte, d.w.z. dat het algoritme de stapgrootte groter maakt wanneer de oplossing minder varieert. òde23' maakt gebruik van 2e en 3e orde formules, òde45' van 4e en 5e orde formules. Met òde45' wordt een hogere nauwkeurigheid bereikt. Beide commando's zijn toepasbaar voor stelsels eerste orde dierentiaalvergelijkingen van de vorm dy = f (t; y) dt

Hierin is y de toestandsvector en f de vectorwaardige functie die de afgeleide van de toestand geeft als functie van t en y. De vorm waarin òde23' gebruikt wordt is: >> [t,x]=ode23('filename',[t0,tm],y0)

Hierin is filename de naam van de .m le waar de dierentiaalvergelijking is gede nieerd, [t0,tm] het interval waarover gentegreerd wordt en y0 de kolomvector met beginwaarden. Het algoritme in het 'programma' ode23 bepaalt eerst in een punt de afgeleide (x_ van de functie gede nieerd in de functie le ` lename'. Hiermee wordt dan het volgende punt bepaald, etc. het algoritme bepaald zelf met welke stappen het interval [t0; tm] doorlopen wordt.

74


Omdat òde45' met hogere orde formules werkt zijn hierbij veelal minder integratiestappen nodig. Als gevolg hiervan geeft òde45' dus soms minder `gladde' plaatjes dan òde23'. De output van òde' geeft een vector t waarin de tijdstippen waarop de oplossing x is uitgerekend (de stappen waarmee het interval [t0; tm] doorlopen is) en een rij(vector) x waarin in de bijbehorende functiewaarden zijn uitgerekend. We zullen e.e.a. illustreren met een voorbeeld waarin òde23' wordt gebruikt. Voor òde45' gaat alles analoog.

Voorbeeld 6.1 Beschouw de dierentiaalvergelijking x_ = 1

tx

Om te beginnen maken we de le divgl.m met als inhoud function xdot=diffvgl(t,x) xdot=1-t*x;

De punt-komma voorkomt dat de output op het scherm komt. De opdracht >> [t,x]=ode23('diffvgl',[0,1],3)

berekent nu voor het tijdsinterval [0; 1] waarden van x(t) en met beginwaarde x(0) = 3. Men kan de oplossingen plotten met >> plot(t,x)

In bovenstaande wordt de output van de ode-solver toegekend aan de vextoren t en x. Er bestaat echter ook de eigenschap ÒutputFcn' waarmee men de ode-solver kan aangeven wat er met de output moet gebeuren. De standaard (default) OutputFcn is òdeplot'. D.w.z. dat de berekende waarden worden geplot zodra zij zijn berekend, de plot wordt voortdurend aangepast. Het commando >> ode23('diffvgl',[0,1],3)

levert dus uiteindelijke hetzelfde resultaat als 'plot(t,x)'.

Oefening 6.3 Bereken met òde23' de oplossing van de dierentiaalvergelijking x_ = 1

2tx

met beginwaarde x(0)=1, op het interval [0,3]. Plot het resultaat. Doe dit op de manier van voorbeeld 6.1. Wat is het resultaat indien men de toekenning [t,x]= weglaat?

2. DIFFERENTIAALVERGELIJKINGEN

75

Oefening 6.4 Beschouw het beginwaardeprobleem x_ =

tx ; x(0) = 3 : (1 + x2 )

a. Bereken met ode23 x(t) op het interval [0,2]. Hoeveel stappen zijn er gebruikt? Plot de gra ek. b. Wat is de benadering voor x(2) in onderdeel a. ?

2.2 Stelsels dierentiaalvergelijkingen De opdrachten òde23' en òde45' kunnen ook in veel algemenere situaties worden toegepast.

Voorbeeld 6.2 Beschouw het stelsel dierentiaalvergelijkingen x_ 1 = 5x1 x_ 2 = 7x1

2x2 4x2

door de beginwaarden x1 (0) = 2 en x2 (0) = 8 wordt precies een oplossing vastgelegd. Met de opdrachten òde23' en òde45' kunnen dergelijke dierentiaalvergelijkingen worden aangepakt. Om te beginnen maken we de le divgl.m met als inhoud function xdot=diffvgl(t,x) xdot=[5*x(1)-2*x(2);7*x(1)-4*x(2)];

Merk op dat xdot een kolomvector moet zijn. De opdracht >> [t,x]=ode23('diffvgl',[0,1],[2,8])

berekent nu voor het tijdsinterval [0; 1]de waarden van x1 (t) en x2 (t) met beginwaarden x1 (0) = 2 en x2 (0) = 8. De variabele x krijgt als waarde een array met twee kolommen, de eerste kolom bevat de waarden van x1 (t) en de tweede de waarden van x2 (t). Het plotten gebeurt door >> plot(t,x)

Omdat x uit twee kolommen bestaat, worden er twee gra eken getekend. Indien men alleen het commando De opdracht >> ode23('diffvgl',[0,1],[2,8])

76


geeft worden ook de gra eken van x1 (t) en x2 (t) getekend omdat per default de ÒutptFcn' odeplot wordt gebruikt. Merk op dat men de function le bijvoorbeeld ook alsvolgt kan programmeren function xdot=diffvgl(t,x) A=[5,-2;7,-4]; x=[x(1);x(2)]; xdot=A*x;

2.3 Het richtingsveld De oplossingen van de dierentiaalvergelijking in voorbeeld 6.2 zijn krommen in de (t; x1 ; x2 )ruimte gegeven door (t; x1 (t); x2 (t). Het richtingsveld wordt gegeven door de raakvectoren aan de oplossingskrommen, dus door de vectoren (1; x_ 1 (t); x_ 2 (t)). Merk op dat men het richtingsveld kan a eiden met behulp van de dierentiaalvergelijkingen, de oplossingen hoeft men niet expliciet te weten. In het algemeen tekent men de oplossingen en het richtingsveld alleen in het (x1 ; x2 )-vlak. Men tekent feitelijk de projectie op het (x1 ; x2 )-vlak. Men noemt dit vlak ook wel het fasevlak of toestandsruimte (de toestand van het systeem op tijdstip t wordt gegeven door (x1 (t); x2 (t))). De geprojecteerde oplossing wordt een fase-pad of ook wel (deel van) een integraalkromme genoemd. In een gegeven punt (x1 (t0 ); x2 (t0 ) van de toestandsruimte kan men de rictingsvector uitrekenen m.b.v. de gegeven dierentiaalvergelijkingen. In het voorbeeld 6.2 volgt voor het punt (x1 (t0 ); x2 (t0 ) = (1; 2) de richtingsvector (x_ 1 (t0 ); x_ 2 (t0 ) = (5 1 2 2; 7 1 4 2) = (1; 1). Het richtingsveld krijgt men door in ieder punt de bijbehorende richtingsvector te bepalen. Om een richtingsveld te kunnen tekenen moet men eerst in het (x1 ; x2 )-vlak de punten aangeven waar men de richtingsvectoren wil tekenen. Hiervoor kan men het commando `meshgrid' gebruiken. Het commando >> [X,Y]=meshgrid(-1:.1:1 ,-2:.2:2)

maakt een matrix X waarvan iedere rij gelijk is aan de vector -1:.1:1 en een matrix Y waarvan iedere kolom gelijk is aan de vector -2:.2:2. Het elementsgewijs combineren van deze twee matrices geeft de roosterpunten (x; y), waarbij x de vector -1:.1:1 doorloopt en y de vector -2:.2:2 doorloopt, d.w.z. de punten (Xij ; Yij ). De vectoren X en Y kunnen nu gebruikt worden om berekeningen in de roosterpunten uit te voeren. Bijvoorbeeld U = 5 X 2 Y en V = 7 X 4 Y berekent voor het voorbeeld 6.2 in ieder roosterpunt de componenten (U; V ) van de richtingsvectoren. Om een richtingsveld in een serie roosterpunten te tekenen kent MATLAB het commando `quiver'. Het commando >> quiver(X,Y,U,V)

3. HET SYMBOLISCH OPLOSSEN VAN DIFFERENTIAALVERGELIJKINGEN

77

tekent de vectoren met componenten (Uij ; Vij ) in de punten (Xij ; Yij ). De vectoren worden automatisch geschaald. Hierdoor wordt het richtingsveld soms wat onduidelijk. Bij de dierentiaalvergelijking uit voorbeeld 6.2 worden de vectoren U en V verkregen door x_ 1 en x_ 2 te berekenen in de punten van het rooster, d.w.z. U = 5 X 2 Y en V = 7 X 4 Y .

Oefening 6.5 Teken het richtingsveld voor de dierentiaalvergelijking uit voorbeeld 6.2 op het rooster gegeven door de vectoren x = [ 1 : :1 : 1] en y = [ 1 : :1 : 1].

2.4 Plotten van integraalkrommen De output van ode23 voor een stelsel van twee vergelijkingen als in voorbeeld 6.2 kan ook gebruikt worden voor het tekenen van fase-paden. Hiertoe moet men de oplossingen x1 (t) en x2 (t) plotten in het fasevlak, d.w.z. het (x1 ; x2 )-vlak. Dit kan door na ode23 de kolommen van de matrix x tegen elkaar uit te zetten met >> plot(x(:,1),x(:,2))

Dit geeft een guur met een fasepad. Wil men verschillende fasepaden met verschillende plot-opdrachten in dezelfde guur laten tekenen dan moet men na de eerste plot-opdracht het commando `hold on' geven. (zie x14). Tekent men fasepaden en richtingsveld in een guur dan zullen de richtingsvectoren raken aan de fasepaden. Men kan ook een fasepad laten tekenen door de ode-solver als eigenschap 'OutputFcn' mee te geven 'odephas2'. Het commando >> ode23('diffvgl',[0,1],[2,8],odeset('OutpuFcn','odephas2'))

resulteert in een 2D fasepad vanuit de gegeven beginvoorwaarde.

3 Het symbolisch oplossen van dierentiaalvergelijkingen 3.1 Eerste-orde dierentiaalvergelijkingen Voor het symbolisch oplossen van dierentiaalvergelijkingen kent de symbolic toolbox van MATLAB het commando 'dsolve'. >> dsolve('Dx1=a*x1') ans = exp(a*t)*C1

men kan hierbij eventueel een beginvoorwaarde toevoegen.

78


>> dsolve('Dx1=a*x1','x1(0)=2') ans = 2*exp(a*t)

Merk op dat in de meeste gevallen de oplossing van een dierentiaalvergelijking niet op te schrijven is als een nette formulefunctie. In deze gevallen heeft symbolisch oplossen geen zin en is men aangewezen op numerieke technieken.

Oefening 6.6 Ga na of de dierentiaalvergelijkingen uit voorbeeld 6.1 en oefening 6.4 symbolisch op te lossen zijn. Merk op dat nieuwe functies worden gentroduceerd. Ga met `help' na wat dit voor functies zijn.

3.2 Stelsels eerste-orde dierentiaalvergelijkingen Er kunnen zonodig meerdere vergelijkingen en beginvoorwaarden worden ingevoerd. In dit geval moet men aangeven wat de output variabelen zijn. >> [x1,x2]=dsolve('Dx1=5*x1 - 2*x2','Dx2=7*x1 - 4*x2','x1(0)=0','x2(0)=5') x1 = 2*exp(-2*t)-2*exp(3*t) x2 = -2*exp(3*t)+7*exp(-2*t)

Indien men de functies x1 en x2 wil plotten kan dit met 'ezplot(x1)' resp. 'ezplot(x2)'. Wil men de toestand (x1; x2) plotten dan dient men eerst de oplossingen x1 en x2 voor een aantal waarden van t te berekenen. Bijvoorbeeld: >> >> >> >>

t=[0:.1:1]; X=subs(x1,t); Y=subs(x2,t); plot(X,Y)

Oefening 6.7 Werk voorbeeld 6.2 uit. Plot in het fasevlak (in een guur) de oplossingskrom-

men met beginwaarden (x1 (0); x2 (0)) resp. (2,8), (2,7), (2,6), (0.1,0.1), (-2,-8), (-2,-7), (-2,-6), (-0.1,-0.1) voor t 2 [0; 1]. (Het verdient aanbeveling hierbij met àxis' de assen vast te kiezen, bijv. met x1 en x2 tussen -10 en 10.) Doe dit zowel numeriek als symbolisch.

4. MATLAB EN DIFFERENTIAALVERGELIJKINGEN - OPDRACHTEN

79

3.3 Hogere-orde dierentiaalvergelijkingen Men heeft binnen matlab ook de mogelijkheid om hogere orde dierentiaalvergelijkingen symbolisch op te lossen zonder deze eerst om te schrijven tot een stelsel eerste orde dierentiaalvergelijkingen. Bijvoorbeeld de direntiaalvergelijking my00 (t) + cy0 (t) + ky(t) = 0

met beginvoorwaarden y(0) = 1, y0 (0) = 0 kan men oplossen met het commando >> dsolve('m*D2y+c*Dy+k*y=0','y(0)=1','Dy(0)=0')

Oefening 6.8 Beschouw de dierentiaalvergelijking y00 (t) + 2ky0 (t) + 25y(t) = 0 ; y(0) = 1 ; y0 (0) = 1 :

a. Bepaal de oplossing van deze vergelijking voor k=0. De oplossingsfunctie beschrijft een trilling. Je kunt dit visualiseren door de oplossing te plotten. Is de trilling gedempt? Is de trilling periodiek? Wat is de frequentie van deze oplossing? b. Bepaal de oplossing van deze vergelijking voor k=1. Wat is de eigenfrequentie van het systeem? De oplossingsfunctie beschrijft een trilling. Je kunt dit visualiseren door de oplossing te plotten. Is de trilling gedempt? Is de trilling periodiek? Wat is de frequentie van deze oplossing? Plot de oplossing en de functie e kt in dezelfde guur. c. Bepaal de oplossing van de vergelijking y00 (t) + 2y0 (t) + 25y(t) = cos(5t) y(0) = 1 ; y0 (0) = 1 :

De oplossingsfunctie beschrijft een trilling. Je kunt dit visualiseren door de oplossing te plotten. Is de trilling gedempt? Is de trilling periodiek? Wat is (op den duur) de frequentie en ampitude van deze oplossing?

4 MATLAB en dierentiaalvergelijkingen - Opdrachten We besluiten deze training met een aantal opdrachten waarin de in voorgaande paragrafen behandelde matlabcommando's een rol spelen. U dient zelf de te volgen strategie te bepalen en te beslissen of u numeriek dan wel symbolisch wilt werken.

4.1 Opdrachten Oefening 6.9 Beschouw de dierentiaalvergelijking van een harmonische oscillator gegeven door een gedempt massa-veer systeem (vgl. [4, hfdst 3]). my00 (t) + cy0 (t) + ky(t) = 0

80


We gaan er van uit dat we een veer hebben die door een gewicht van 98 Newton 20 cm wordt uitgerekt. Dus k =490[N/m] en m =10[kg]. De demping c is zo gekozen dat c2 < 4mk, d.w.z. onderkritische demping. 10y00 (t) + 50y0 (t) + 490y(t) = 0 Deze dierentiaalvergelijking is desgwenst te schrijven als een stelsel eerste orde dierentiaalvergelijkingen. Stel x1 (t) = y(t) en x2 (t) = y0 (t). Dan krijgen we de vergelijkingen: x_ 1 = x2 x_ 2 = 5x2

49x1

Plot de oplossing y(t) met beginwaarde y(t) = 0:2, y0 (t) = 0 voor een voldoende groot tinterval en plot deze oplossing ook in het fasevlak. Wat gebeurt er bij andere beginwaarden? Plot ook voor bovenkritische demping, c > 140, de oplossing met dezelfde beginvoorwaarden. Kijk of je bij kritische demping, c = 140, beginvoorwaarden kunt vinden zodat de oplossing tenminste een nulpunt heeft (tenminste een doorgang door de evenwichtsstand).

Oefening 6.10 Beschouw de Lotka-Volterra vergelijkingen voor een prooi-roofdier model (voor meer informatie zie [5]) dp = (a1 b1 r)p dt dr = ( a2 + b2 p)r dt

Met a1 ; a2 ; b1 ; b2 > 0. Hierin is p(t) het aantal prooidieren en r(t) het aantal roofdieren. De prooidieren hebben een onuitputtelijke voedselbron en er is geen onderlinge concurrentie. Wanneer er geen roofdieren zijn ontwikkelt de populatie zich exponentieel. De prooidieren vormen voor de roofdieren de enige voedselbron. Wanneer er geen prooidieren zijn sterft de roofdierpopulatie uit. Men kan aantonen dat p(t) en r(t) aan de volgende relatie moeten voldoen a2 log p

b2 p + a1 log r

b1 r = constant

Deze vergelijking geeft voor iedere constante een gesloten kromme in het (p; r)-vlak. De populaties ontwikkelen zich dus volgens een periodiek patroon. a. Neem de waarden a1 = 1:00, b1 = 0:10, a2 = 0:50 en b2 = 0:02. Plot de oplossingen p(t), r(t) en de bijbehorende kromme in het fasevlak voor enige zelfgekozen beginwaarden en een zelf gekozen tijdsinterval. Merk op dat het punt p = ab11 = 25, r = ab22 = 10 een stationair punt is in het fasevlak (constante oplossing). Opmerking: Vanwege numerieke oorzaken kan het zijn dat de getekende fasekromme niet een gesloten kromme is maar een spiraalkromme. Dit komt alleen naar voren bij plotten met grote t-intervallen. òde45' geeft hier het beste resultaat. b. De parameters a1 en a2 bepalen de ontwikkeling wanneer er geen interactie is. De parameters b1 en b2 bepalen de interactie tussen roofdier en prooi. Bekijk voor een aantal verschillende parametercombinaties hoe de populaties zich ontwikkelen. Merk op dat alleen het eerste kwadrant van het (p; r)-valk van belang is.

4. MATLAB EN DIFFERENTIAALVERGELIJKINGEN - OPDRACHTEN

81

Oefening 6.11 Beschouw een aangedreven veer-demper systeem gegeven door de vergelijking

my + cy_ + ky = cx_ + kx

Hierin is x de aandrijving of input en y de responsie of output. m=870 kg, k=70000 Nm 1 , c=5000 Nsm 1 (men kan denken aan de vering van een auto). We veronderstellen voor de input x(t) = sin(t). Bereken (kies zelf beginvoorwaarden, niet (0,0)) en teken de output (kies een interval waarop het periodiek gedrag zichtbaar is).

Oefening 6.12 Beschouw als in de voorgaande opdracht een aangedreven veer-demper systeem gegeven door de vergelijking

my + cy_ + ky = k sign(sin(2t)) m=870 kg, k=70000 Nm 1 , c=5000 Nsm 1 . Ga na wat de functie in het rechterlid precies voorstelt. Merk op dat de matlab-functie 'sign' geen symbolisch argument accepteert. Bereken (kies zelf beginvoorwaarden, niet (0,0)) en teken de output (kies een interval waarop het periodiek gedrag zichtbaar is).

82


Bibliogra e [1] The Student Edition of MATLAB, Version 4, User's Guide / The Mathworks Inc., Prentice Hall, Englewood Clis, N.J., 1995. [2] Studiehandleiding Lineaire Algebra voor W (2Y650), diktaat nr. 2532, Technische Universiteit Eindhoven, 1997. [3] B. Kolman, Introductory Linear Algebra with Applications, 6th ed.,Prentice-Hall internatinal Inc., Upper Sadle River, 1997. [4] Syllabus Calculus voor W (2Y130), Technische Universiteit Eindhoven, 1998. [5] E.C. Pielou, Mathematical Ecology, John Wiley & Sons, New York, 1977. [6] E. Part-Enander et al, The MATLAB handbook, Addison-Wesley, Harlow, 1996.

83

Deel 2 Handleiding Matlab voor W en BMT trimester 1.3, dagdeel 5 H.A. van Essen Faculteit Werktuigbouwkunde 2000/2001

Programmeren in Matlab In dit deel worden eerst enkele algemene opmerkingen over het programmeren in MATLAB vermeld. Daarna worden de twee opgaven van middag 5 gepresenteerd. Je kunt beginnen met het doornemen van deze handleiding en het uit proberen van de voorbeelden op je eigen computer. Probeer zoveel mogelijk de tekst van de MATLAB M-files over te nemen uit deze handleiding, (selecteer hiervoor in Acrobar Reader 4 het “T” icoontje en in versie 3 het “ABC” icoontje) dat bespaart veel overtypwerk. Natuurlijk mag je ook andere voorbeelden en variaties uitproberen. De opgaven die je na het doorwerken van het eerste gedeelte moet maken bouwen voort op de gegeven voorbeelden: probeer deze voorbeelden dus eerst te begrijpen.

1 Command window en directories Wanneer je MATLAB opstart, krijg je het command window te zien. Hierin staat een prompt >> waar achter je commando’s kunt intypen die je wilt uitvoeren. In deze training gaan we onze eigen MATLAB files schrijven en nieuwe commando toevoegen aan MATLAB. Het is zeer onverstandig om je eigen files in de standaard directory van MATLAB te plaatsen. De standaard directory van MATLAB is vanaf Realease 11 C:\matlab\work of C:\ProgramFiles\matlab\work. Deze directory staat tussen de MATLAB programmadirectories. Je kunt beter een eigen directory aanmaken voor je nieuwe files. Je kunt zien in welke directory je staat met het commando pwd (van print working directory). Met het commando ls (van list) kun je zien welke files in deze directory staan. Het is mogelijk dat u in een eerdere training al een eigen directory heeft aangemaakt, bijvoorbeeld D:\Matlabtr, natuurlijk kunt u die gebruiken. MAAK ANDERS EEN NIEUWE DIRECTORY AAN VOOR JE EIGEN MATLAB FILES. Dit kun je het beste doen in de Windows Explorer of Verkenner, dus buiten MATLAB zelf. Als je eigen files in een eigen directory bewaart, moet Matlab die natuurlijk wel kunnen vinden. Daartoe kijkt MATLAB altijd eerst in de huidige directory (herinner het commando pwd). Om MATLAB duidelijk te maken in welke directory jouw eigen files gevonden kunnen worden, kun je twee dingen doen: 1. Je kunt naar je eigen directory toelopen met het commando cd (van change directory), zodat dit de huidige directory wordt. Bijvoorbeeld cd C:\piet\matlab\training3. Je kunt ook stapsgewijs van de ene naar de andere directory lopen. Vanuit C:\matlab\bin naar C:\piet\matlab\training type je dan achtereenvolgens: >> cd .. (dit betekent een stap omhoog, let op: cd spatie punt punt) >> cd .. >> cd piet >> cd matlab >> cd training3

Als je "lange namen" of namen met spaties kiest, dan moet je het path met apostroffen opgegeven, anders herkent MATLAB het niet. Bijvoorbeeld: cd 'C:\My Matlab Directory'.

2

PROGRAMMEREN IN MATLAB

2. Je kunt ook je eigen directory toe aan het MATLAB-path, d.w.z. de verzameling van directories waarin MATLAB automatisch zoekt naar files. Je kunt je eigen directory toevoegen aan het standaard path m.b.v. het commando addpath of m.b.v. de pathbrowser. Deze kun je opstarten met een button op de werkbalk van het command window. (Probeer sowieso de vier meest rechtse buttons eens uit: workspace browser, path browser, opstarten SIMULINK, MATLAB help). Als je de path browser voor het eerst opent, moet je hem vaak wat groter maken. Klik vervolgens op de browse button en selecteer in het venster rechts de directory waar je je eigen files wilt gaan opslaan. Als je op OK drukt, dan verschijnt de gekozen directory in het vakje current directory linksboven. Klik vervolgens op de tweede button van links "add to path" van de werkbalk of selecteer "add to path " onder het menu path bovenaan. Je kunt een directory eenmalig of permanent toevoegen aan het path door het path al dan niet te saven. PROBEER DIT EENS UIT.

2 MATLAB editor en workspace Je eigen MATLAB files kun je schrijven in de M-file editor. Deze kun je starten door in de command window te typen >>edit

of >>edit filename

(alleen voor bestaande files als je al in de goede directory zit)

of met de muis op de werkbalk: file, new M-file. In de editor kun je je eigen file ook saven. LET ER OP DAT JE DE GOEDE DIRECTORY SELECTEERT. De MATLAB workspace is het computergeheugen waarin alle variabelen staan opgeslagen die op dat moment gedefinieerd zijn. Met het commando whos kun je kijken welke variabelen op dat moment in het geheugen staan en welke afmeting (of dimensies) ze hebben. De werkelijke inhoud of waarde van de variabelen kun je inzien door de naam van de variabele in te typen in de MATLAB command window. Je kunt het geheugen schoon (leeg) maken m.b.v. het commando clear all.

3 MATLAB help MATLAB kent uitgebreide ingebouwde help faciliteiten. Wanneer je weet welk commando je moet gebruiken, maar niet meer precies weet wat de syntax of wat de volgorde van de argumenten is, dan kun je (vrijwel) altijd >> help commandonaam

intypen achter de prompt van het command window. Wanneer je niet precies weet welk commando je nodig hebt, of als je nieuwsgierig bent naar de vele ingebouwde commando’s die MATLAB kent, dan kun je met >> helpwin

2


3

browsen door de helpfiles van MATLAB. Als je ook de HTML of PDF helpfiles op je computer hebt geïnstalleerd, dan kun je ook de HelpDesk gebruiken. Je kunt deze opstarten vanuit het Help menu op de werkbalk van het command window.

4 Enkele opmerkingen over variabelen • Variabelen in MATLAB worden opgeslagen als arrays of matrices, je kunt de inhoud elements-gewijs adresseren. Een dubbele punt (:) geeft een hele rij of kolom aan. >>A=[1 2; 3 4]; >>A(2,1) ans=3 >>A(:,2) ans= 2 4 >>A(1,:) ans= 1 2

•

•

Variabelen - zijn case-sensitive (dus Items niet gelijk aan ITEMS) - zijn maximaal 31 karakters lang - beginnen met een letter (dus dit mag: X51483 en dit mag niet: 1abcd) - met de namen: ans, pi, eps, flops, inf, NaN (nan), i, j,nargin, nargout, realmin, realmax hebben een speciale betekenis of waarde in MATLAB. Vermijd deze namen dus! Display formaat >>x=1/3 x=0.3333 >>v = x-0.3333 v = 3.3333e-005 >>format long >>x x=0.33333333333333 >>v v=0.00003333333333

•

Je mag je eigen files ook geen willekeurige namen geven. Matlab commando’s zoals for of end, maar ook prod en eig mag je NIET kiezen. Soms krijg je geen echte foutmelding maar alleen onbegrijpelijke resultaten.LET HIEROP.

5 Programmeren Programmeren is het “vertalen van een technische vraag naar computer-instructies” of te wel het schrijven van een computer-programma. Een standaard procedure voor het oplossen van een technisch probleem met behulp van programmeren in MATLAB is als volgt 1. 2. 3. 4. 5.

probleem analyseren, probeer de oplossystematiek te bepalen (op papier) formulering uitwerken (op papier) M-file schrijven in de MATLAB editor testen en debuggen probleem oplossen

In deze training gaan we vooral in op de punten 3 en 4 met betrekking tot MATLAB. Daarnaast gaan we in op enkele belangrijke elementen en constructies uit programmeertalen die je nodig hebt om een goed programma te schrijven.

3

4


6 Twee soorten M-files In MATLAB kun je twee soorten programma’s maken. Omdat beide types de extensie .m hebben, noemen we ze M-files. We kennen de script M-file, en de function M-file. In het algemeen zijn M-files • niet te begrijpen voor anderen… • niet meer begrijpbaar voor jezelf na ongeveer twee weken… …tenzij uitvoerig van commentaar voorzien (met een % teken voor de regel): % dit is een commentaar regel in een M-file % en wordt dus niet uitgevoerd als commando % helpt alleen om te begrijpen wat er staat. 6.1 script M-file Een script M-file voert van te voren ingevoerde programmaregels (die ook achter elkaar achter de prompt in het command window kunnen worden ingetypt) in een keer uit. Handig als je hetzelfde programmaatje meerdere malen wilt uitvoeren (voor verschillende parameterwaarden bijvoorbeeld) of wanneer je het programma wilt bewaren. Een script M-file is een op zichzelf staand programma dat je steeds opnieuw kunt aanroepen door de naam van de M-file achter de prompt in de command window in te typen Voorbeeld 1: Als je een M-file schrijft met daarin cd cd cd cd cd

.. .. piet matlab training

en die opslaat onder de naam piet.m, dan kun je door het intypen van het commando >> piet

in één keer van de opstart directory C:\matlab\bin naar C:\piet\matlab\training. Voorbeeld 2: plotten van de functie y = sin (x ) + cos(x ) in één commando •

schrijf een M-file met als naam plotsplusc.m waarin staat % maak het geheugen schoon clear all % maak een rij met x-waarden x = [-10:0.1:10]; % bereken de bijhorende functie waarden y = sin(x) + cos(x); % teken y als functie van x plot(x,y)

4

PROGRAMMEREN IN MATLAB •

5

Sla deze file op in je eigen directory. Gebruik deze M-file (d.w.z. uitvoeren van de commando's die in de file staan)) door gewoon >> plotsplusc

in te typen. Dit noemen we het "aanroepen" van de file. 6.2 function M-file Met behulp van een function M-file kun je zelf nieuwe MATLAB commando's schrijven. Function M-files zijn te vergelijken met een function of subroutine in andere programmeertalen. Op deze manier kun je nieuwe commando's maken, door het combineren van bestaande. Een function M-file staat niet op zichzelf, maar communiceert via in- en uitgangsvariabelen: dit noemen we de argumenten van de functie. Voorbeeld 3: y = sin (x ) + cos(x ) •

Schrijf een M-file met als naam sinpluscos.m waarin staat

function y=sinpluscos(x) % aanroepen met y = sinpluscos(x) % deze functie berekent de som van de sinus % en de cosinus van x en stopt deze waarde in y. a=sin(x); b=cos(x); y=a+b;

waarin − x is de ingangsvariabele (argument) − y is de uitgangsvariabele (argument) − a en b zijn lokale en daarmee tijdelijke variabelen: deze komen niet in de MATLAB workspace te staan, ze hebben achteraf dus geen waarde. •

sinpluscos is nu een nieuw MATLAB commando geworden. Je kunt dit commando nu gewoon gebruiken. Let op dat je wel argumenten mee moet opgeven! Bijvoorbeeld:

>> sinpluscos(pi/4) ans = 1.4241

of >> piet = sinpluscos(pi/4) piet = 1.4241

of x = [-10:0.1:10]; y = sinpluscos(x);

% maak een rij met x-waarden % bereken de bijhorende functie waarden

5

6

PROGRAMMEREN IN MATLAB plot(x,y);

% teken y als functie van x

Let op: - de functienaam moet gelijk zijn aan de bijbehorende filenaam (behalve .m extensie). - voor het aanroepen van de function (het nieuwe MATLAB commando) hoef je niet dezelfde namen van de argumenten (in dit geval x en y) als in de file te gebruiken, maar het mag wel. - de tweede regel van de functie, de commentaarregel, is automatisch de help-tekst voor de functie. Als je deze regel invoert, kun je dus >> help functienaam

deze tekst op je scherm krijgen. Vaak is het handig op deze plaats iets te schrijven over wat de functie precies doet en hoe de argumenten gekozen moeten worden. In dit geval geeft dit: >> help sinpluscos aanroepen met y = sinpluscos(x) deze functie berekent de som van de sinus en de cosinus van x en stopt deze waarde in y.

7 Programmeertaal constructies In deze training maken we kennis met twee speciale constructies die veel gebruikt worden in computerprogramma's: de repeterende constructie en de voorwaardelijke constructie. 7.1 repeterende constructies: FOR loop In MATLAB kennen we twee repeterende constructies. De eerste is de FOR loop en de tweede is de WHILE loop. Nu bespreken we eerst de FOR loop. Deze constructie biedt de mogelijkheid bepaalde stukjes programma meerdere malen achter elkaar uit te voeren waarbij de waarde van bepaalde variabelen steeds verandert. Een zogenaamde loopvariabele houdt bij hoe vaak de programmaregels moeten worden herhaald. De algemene vorm van een FOR loop is FOR loopvariable = start : step : end programmaregels END Wanneer step wordt weggelaten wordt automatisch stapgrootte 1 aangenomen. Voorbeeld 4 : een FOR loop: % definieren van constanten n n = 10; % aanmaken van een matrix met n rijen % en 1 kolom met allemaal nullen gevuld A=zeros(n,1);

6


7

% vullen van deze kolom met de waarden 1 tot % en met n in een FOR loop for k=1:n A(k,1) = k; % hier staat: het element van A op de k-de % rij en de 1-ste kolom krijgt de waarde k end A % afdrukken resultaat op het scherm.

De werking van deze loop is als volgt: na het for commando krijgt de loopvariabele k de waarde 1 en worden de programmaregels tot het end commando uitgevoerd. Dat wil zeggen dat A(1,1) gelijk wordt aan 1. Na het end commando springt MATLAB terug naar de for regel en krijgt de loopvariabele k de waarde 2. Vervolgens worden de programmaregels weer uitgevoerd: A(2,1) wordt 2. Dit herhaalt zich tot de laatste stap, waarin de loopvariabele k de waarde van n dus 10 krijgt. Hierna gaat het programma verder met de regels na end. Dan wordt A dus geprint op het scherm. k = 1 : n betekent “voor k is 1 tot en met n in stapjes van 1”. Je kunt ook in grotere stapjes van 1 naar n gaan of beginnen bij een andere waarde: for k=2:2:10 betekent “voor k is 2 tot en met 10 in stapjes van 2” (dus 2,4,6,8,10). Natuurlijk kun je ook andere variabelen kiezen als i en n. Voorbeeld 5 % initieren van parameter p p = 0; for k=1:10 p = p + 1; end p

In deze loop wordt telkens de waarde van p met één opgehoogd. De nieuwe waarde van p gelijk wordt aan de oude waarde van p (dat was dus 0) plus 1, dus p wordt 1. DIT IS EEN ASSIGNMENT OPDRACHT EN DUS GEEN VERGELIJKING. Zo'n assignment constructie kun je heel goed gebruiken om een teller of een sommatie binnen een loopconstructie bij te houden. Het is ook mogelijk om meerdere for loops te nesten. Op zo’n manier kun je makkelijk een nbij-m matrix vullen door element voor element de matrix af te lopen.: Voorbeeld 6: Een matrix vullen met oplopende getallen. % definieren van dimensies van een matrix nr = 4; % aantal rijen nk = 3; % aantal kolommen p=1;

% start getal

% aanmaken van een matrix A met nr rijen % en nk kolommen met allemaal nullen gevuld % in dit geval is het is niet perse noodzakelijk de matrix % reeds van te voren aan te maken met allemaal % nullen omdat je toch alle elementen een waarde geeft, % maar het is wel netjes, het voorkomt fouten, % en soms kun je nuttig gebruik van maken. A=zeros(nr,nk);

7

8

PROGRAMMEREN IN MATLAB % vullen van deze matrix met random waarden van 0 tot 1 % in een geneste dubbele FOR loop for r=1:nr for k=1:nk A(r,k) = p; % hier staat: het element van A op % de r-de rij en de k-de kolom krijgt de % krijgt de waarde p p+p+1; end end A

Ga stap voor stap na hoe dit voorbeeld werkt; Kunt u het eindresultaat A begrijpen? Probeer te begrijpen hoe de matrix rij voor rij gevuld wordt. 7.2 voorwaardelijke constructies Een van de meest gebruikte programmeertaal constructies is de if constructie. Hiermee is het mogelijk om aan de hand van een relationele test op logische variabelen bepaalde programmaregels wel of juist niet uit te voeren. De algemene vorm van deze constructie is IF logische expressie programmaregels ELSEIF logische expressie programmaregels ELSE programmaregels END De ELSEIF en ELSE statements are optioneel, die mag je dus ook weglaten. Een logische expressie is waar of niet-waar, true of false, 1 of 0. Bij het evalueren van een logische expressie maken we gebruik van relationele en logische operatoren:

Relationeel (vergelijkend) < kleiner dan <= kleiner dan of gelijk aan > groter dan >= groter dan of gelijk aan == gelijk aan ~= ongelijk aan

Logisch & | ~

en of niet

Voorbeelden: 1 < 2 WAAR, 1 == 1 WAAR, 1 == 2 NIET WAAR, 1 ~= 2 WAAR etc. Pas op: 1.0000 == 1 soms WAAR, meestal NIET WAAR !!!!!! t<0  f (x ) = 0 als tekent op Voorbeeld 7 : maak een MATLAB programma dat de functie  t≥0  f (x ) = 1 het domein [-10,10]. Schrijf een script M-file met als inhoud: % rij met tijdstappen aanmaken

8


9

t = -10:0.1:10 % bepaal het aantal tijdstappen in array t met het commando size n=max(size(t)) % initialiseer de functiewaarden f op 0 f = zeros(1,n) % bepaal de functie waarden in een FOR loop for i=1:n if t(i) < 0 % let op: t(i) betekent het i-de element % van de rij t. i is een index dwz een % adres, terwijl t(i) de inhoud of % waarde van dat adres bevat. t(i) is in % dit geval hetzelfde als t(1,i). f(i) = 0; else f(i) = 1; end end % teken de functie plot(t,f)

Ga na dat, omdat we f van te voren op de waarde nul gezet hebben, bovenstaande IF END loop eenvoudiger gemaakt kan worden. Voorbeeld 8 : het gebruik van een for loop en conditionele tests voor het berekenen van het inproduct van twee vectoren. Hiervoor schrijven we een nieuwe function M-file, dat wil zeggen, we creëren ons eigen MATLAB commando voor het berekenen van het inproduct. We noemen ons nieuwe commando inprod . Om dit te realiseren schrijven we een function M-file inprod.m. We noemen de twee vectoren die we willen vermenigvuldigen in deze file a en b. Het uitgangsargument noemen we ipr. function ipr = inprod(a,b) % test of a en b vectoren zijn met gelijke lengte %( anders is het niet mogelijk het inproduct te bepalen. [ra ka]=size(a); % het commando size geeft nu de afmetingen van % matrix a door aan de variabelen ra en ka, dit % zijn respectievelijk het aantal rijen en het % aantal kolommen van a [rb kb]=size(b); if ka~=1 error('eerste argument is geen vector') end if kb~=1 error('tweede argument is geen vector') end if ra~=rb error('vectoren kunnen niet vermenigvuldigd worden') end % initialiseren van resultaat (een getal) ipr = 0; for p=1:ra ipr = ipr + a(p,1)*b(p,1); end

9

10


7.3 repeterende constructies: WHILE loop Een tweede repeterende constructie in MATLAB is de WHILE loop. Bij de FOR loop moet je van tevoren aangeven hoeveel keer je de loop wilt uitvoeren. De WHILE loop biedt de mogelijkheid om een loop te combineren met een voorwaardelijke constructie: De loop wordt uitgevoerd zolang een bepaalde expressie waar of niet-waar is. De algemene vorm van een WHILE loop is WHILE logische expressie programmaregels END Voorbeeld 9: Een "eerlijke" loting. Hierin wordt een WHILE loop gebruikt Matlab kiest random een getal tussen 1 en 10. Vervolgens moeten verschillende spelers om de beurt raden welk getal dit is. % laat MATLAB een willekeurig getal tussen 1 en 10 bepalen n = round(rand(1)*10+0.5); % initieren stopvariabele stop = 0; while stop == 0 % % % %

hier staat: zolang de variabele stop gelijk is aan 0, worden de programmaregels tussen while en de bijbehoerende end steeds opnieuw herhaalt. Om te stoppen moet dus BINNEN deze lus de variabele stop een waarde ONGELIJK aan 0 gegeven worden. keuze = input('Geef een geheel getal van 1 tot en met 10 en ENTER:'); % met het commando input kun de gebruiker om input vanaf het toetsenbord % vragen, sluit af met enter

if keuze == n % het juiste getal is gekozen, we kunnen stoppen. stop = 1; disp('*** ***') disp('U heeft gewonnen !!!!') disp('*** ***') else % een fout getal is gekozen, we moeten doorgaan disp('Helaas onjuist: volgende speler') end end

8 Debuggen Wanneer je eenmaal een programma hebt geschreven, kun je het uitvoeren. Meestal doe je dit door de naam van de (script) M-file in te typen in de command window. Soms geeft MATLAB dan een foutmelding of gebeurt er iets anders wat je niet verwacht. Een foutmelding van MATLAB is meestal duidelijk genoeg om je programma te verbeteren. Dit is met name het geval bij syntax fouten, een vergeten END commando of het bewerken

10


11

van matrices die niet de juiste dimensies hebben (bijv. vergeten te transponeren). Wanneer je de fout echter niet herkent, of als je programma wel loopt, maar je onverwachte of rare resultaten krijgt, moet je je programma debuggen. De eenvoudigste manier van debuggen is het bekijken van tussenresultaten. Deze kun je op je scherm laten afdrukken door het weghalen van de ; achter de betreffende regels. Op deze manier is het bijvoorbeeld mogelijk om bij te houden wat er met een bepaalde variabele in een loop gebeurd. Andersom kun je bij een programma dat wel goed werkt overal ; achter zetten om onnodige tussenresultaten te vermijden. Iets ingewikkelder is het gebruik van de MATLAB debugger. Deze is ingebouwd in de editor en staat je toe om stap voor stap door je programma te lopen, zodat je kunt zien waar de fout optreedt. Je kunt ook zogenaamde breakpoints aangeven waarnaar het programma toeloopt en dan pauzeert. Je kunt vervolgens stap voor stap door je programma heen op zoek naar een eventuele fout. De waarde van alle variabelen in je programma kun je tussendoor bekijken door er met de muis op te gaan staan. PROBEER HET GEBRUIK VAN DE DEBUGGER ZELF EENS UIT. Je kunt een breakpoint aanmaken in de editor met een icoontje. Vervolgens start je het programma gewoon op in Matlab. Probeer dit eens met het vorige voorbeeld uit. MAAK NU DE OPGAVEN

Opgaven Opgave 1 MATLAB kent standaard een matrixvermenigvuldiging. De matrices A en B , die we zelf aanmaken in de MATLAB workspace kunnen vermenigvuldigd worden door middel van prod = A*B. Het resultaat is dan de matrix prod. Wanneer de dimensies van de matrices niet overeenkomen, dan genereert MATLAB een foutmelding. In deze opgave gaan we zelf een nieuw MATLAB commando schrijven voor het uitvoeren van een matrix vermenigvuldiging. Een nieuw MATLAB commando kun je maken m.b.v. een M-file function. De volgende functie—een file met als naam matverm.m , opgeslagen in je eigen directory-- zou kunnen voldoen: function resultaat = matverm(mat,kol) resultaat = mat*kol

Deze functie kan worden aangeroepen d.m.v. prod = matverm(A,B).Herinner dat mat en kol de ingangsargumenten van de functie zijn en dat resultaat het uitgangsargument is. De namen in de aanroep hoeven dus niet overeen te komen met de namen in de functie file, alleen de plaats in de aanroep is belangrijk. In de bovenstaande functie gebruiken we nog steeds de ingebouwde matrix vermenigvuldiging. De opdracht is om nu zelf een Matlab functie te schrijven die twee matrices van willekeurige dimensies met elkaar vermenigvuldigd waarbij je alleen gebruikmaakt van de scalaire vermenigvuldigingsoperator * . Met andere woorden, je moet

11

12


de matrixvermenigvuldiging zodanig uitschrijven dat je alleen (vele malen achter elkaar!) twee gewone getallen vermenigvuldigt. • Ga uit van het bovenstaande voorbeeld. Noem je m-file bijvoorbeeld matverm.m. Gebruik als (lokale) namen voor de ingangsargumenten mat en kol of A en B. Noem het uitgangsargument bijvoorbeeld prod of resultaat. • Analyseer het probleem: Ga na hoe je element voor element het product van een matrix met een kolom bepaalt. Let goed op de volgorde van de vermenigvuldiging: A*B ongelijk aan B*A. Ga in eerste instantie uit van twee 2bij2 of 2bij3 matrices om te ontdekken hoe dit ook al weer zat. Het programma moet dit straks voor allerlei verschillende matrix dimensies kunnen, probeer dus de systematiek te ontdekken. Dit is een typisch voorbeeld van programmeren: analyseer het probleem aan de hand van een voorbeeld dat je nog makkelijk kunt bevatten en generaliseer dat vervolgens tot een computerprogramma dat de taak automatiseert voor veel ingewikkelder problemen. • De functie moet achtereenvolgens het volgende doen: − Bepaal de dimensies van de ingangsargumenten mat en kol. Gebruik hiervoor het commando size. − Ga na of deze dimensies geschikt zijn voor vermenigvuldiging, genereer anders een melding op het scherm. Gebruik hiervoor een test met if end, een (fout)melding op het scherm kan geprint worden met error('text') of disp('text'). Ga (met help) na wat deze functies doen en wat de verschillen zijn. − Bepaal de dimensie van de resultaat-matrix op basis van de dimensies van de ingangsargumenten en initialiseer deze matrix met allemaal nullen als elementen. (Hiervoor kun je het commando zeros gebruiken.) − Vul, element voor element, in een for loop, de resultaat matrix. De waarde op elk element is een sommatie van verschillende producttermen(een inproduct). Hiervoor moet je (dus binnen de eerste for loop) nog een for loop gebruiken samen met de assignment constructie som = som + verandering. (de nieuwe waarde van som wordt gelijk aan de vorige waarde van som plus de waarde van de verandering). Zie hiervoor het voorbeeld van het berekenen van het inproduct. − Zorg ervoor dat de resultaat-matrix het uitgangsargument van de functie is. • Probeer de functie uit voor allerlei matrices A en kolommen B die je zelf aanmaakt. Controleer je programma door het resultaat te vergelijken met de ingebouwde matrixvermenigvuldiging van Matlab. Controleer ook de foutmeldingen. •

VOOR DE LIEFHEBBERS, breid je programma uit naar een vermenigvuldiging van twee matrices.

Opgave 2 In deze opgave moet je in een MATLAB M-file een menu aanmaken. Een menu is in vele gevallen handig om afhankelijk van de (waarde van de) input van de gebruiker een bepaalde actie te ondernemen. In Matlab zijn vele mogelijkheden om een dergelijke menu-structuur te programmeren. Vandaag doen we dit recht toe recht aan met het IF ELSEIF ELSE END commando. Voor het printen van tekst op het scherm kun je het commando disp gebruiken. Voor het vragen naar input van de gebruiker vanaf jet toetsenbord kun je het commando input gebruiken. Zoek het juiste gebruik van de commando's op in de ingebouwde help of kijk nog even in de gegeven

12


13

voorbeelden. Maak je programma zo flexibel mogelijk, dan kun je het later nog eens voor serieuzere doelen gebruiken! Als eenvoudig voorbeeld kiezen we nu het plotten van een functie in verschillende kleuren en of lijn soorten. Het gaat nu niet om wat we plotten, maar hoe we het plotten. We kijken naar plot-opties zoals kleuren en lijnen, assen e.d. en naar de mogelijkheid om meerdere plots in een figuur of op een pagina te krijgen. Eerst wat informatie over de verschillende plotmogelijkheden van MATLAB. %voorbeeld voor het plotten van rode rondjes x = [-2:0.1:2] plot(x,x.^2, 'ro') %voorbeeld voor het plotten van groene streep stippellijn plot(x,x.^2, 'g-.')

Als je het bovenstaande commando uitprobeert, zie je dat MATLAB standaard het vorige plot-figuur overschrijft waneer je een nieuwe plot opdracht opgeeft. De standaard naam voor het plot-window is figure(1). Wanneer je de vorige plot wilt bewaren, bijvoorbeeld om twee versies te vergelijken, dan kun je een nieuw plot-window openen door figure(2) in te typen. % voorbeeld voor het plotten van rode rondjes % en groene streep stippellijn in aparte windows x = [-2:0.1:2] figure(1) plot(x,x.^2, 'ro') figure(2) plot(x,x.^2, 'g-.')

Meerdere lijnen in één figuur kun je krijgen door plot opdrachten te combineren, of door gebruik van het hold commando. Hold houdt het huidige figuur open voor nieuwe plotopdrachten. Vergeet niet het figuur af te sluiten met hetzelfde commando hold. % voorbeeld voor het plotten van rode rondjes % en groene streep stippellijn in een window x = [-2:0.1:2] % combineren van plotopdrachten plot(x,x.^2, 'ro', x,x.^2,'g-.') % gebruik van het hold commando plot(x,x.^2, 'ro') hold plot(x,x.^2,'g-.') hold

Ga met help plot zelf na welke andere plotkleuren en markeringen in Matlab gebruikt kunnen worden. Grootheden op de assen kunnen met de commando's xlabel en ylabel worden opgegeven. Een titel kan met title worden opgegeven. Eventueel kun je met het commando legend een legenda opgeven. Verschillende grafieken op één pagina / in één window kun je maken met het commando subplot. subplot(m,n,p) deelt het huidige figuur window op in een m-bij-n matrix van deelfiguren en selecteert het p-de deelfiguur om iets in te plotten. De figuren worden rij voor rij geteld van linksboven naar rechtsonder.

13

14


% voorbeeld voor het plotten van rode rondjes % en groene streep stippellijn in aparte subplots x = [-2:0.1:2] subplot(2,1,1); plot(x,x.^2, 'ro') subplot(2,1,2); plot(x,x.^2,'g-.')

Opdracht menu file: Schrijf een menu-file. Onze menu-file is een script M-file die zodra die aangeroepen wordt (naam intypen en “enter” geven achter de Matlab prompt) − een aantal mogelijkheden op het scherm weergeeft samen met een vraag naar welke mogelijkheid je kiest, − vervolgens wacht tot de gebruiker een keuze heeft opgegeven en via enter bevestigd heeft. − de input van de gebruiker controleert en de gewenste opdracht uitvoert. − een foutmelding geeft indien de input van de gebruiker niet overeenkomt met de mogelijkheden en de vraag daarna opnieuw stelt. − de vraag opnieuw stelt na het uitvoeren van de opdracht om eventueel een andere keuze uit te voeren. Het is daarom handig ook een “stop criterium” in te bouwen, bijvoorbeeld door een stop = 0 while stop ~= 1 … end loop gebruiken, waarin stop een parameter is die in het begin de waarde 0 heeft en die de waarde 1 krijgt als de gebruiker het stopcriterium kiest. Zolang stop gelijk aan nul is worden de programma regels tussen while en end continu herhaalt. Kies zelf een functie om te plotten en een aantal leuke plotopties. Denk naast verschillende kleuren en lijnen ook eens aan subplots of meerdere figuren. Voorbeeld: In welke kleur of markering wilt u de functie plotten? toets toets toets toets

1 2 3 4

voor rode rondjes voor blauwe sterretjes voor gele lijnen om te stoppen

Wat is uw keuze?

Extra uitbreidingsmogelijkheden voor deze menufile In Matlab bestaan handige functies voor het snel analyseren van wiskundige functies. De commando's fplot, fmin en fzero worden gebruikt om functies te plotten, minima en nulpunten te bepalen. Probeer onderstaande commando’s eens uit en probeer daarna deze functies toe te voegen aan je eigen menu-file. Er wordt dan bijvoorbeeld gevraagd om de functie die je wilt analyseren zelf in te geven en te kiezen of je een bepaald bereik wilt plotten, nulpunten of minima wilt bepalen. % functie definitie, let op gedefinieerd als een % TEXT string tussen apostrofes! HIEROM KUN JE OOK VRAGEN IN EEN MENU !!! f='x^4-100*x^2+3/10*x-6'; % plotten van deze functie tussen x=-15 en x=+15. % de y-schaal wordt automatisch aangemaakt. xl=-15; xh=15; fplot(f, [xl xh])

14


15

% bepalen minimum tussen x=-15 en x=+15. xmin = fmin(f, xl, xh) % in de plot is te zien dat dit minimum niet uniek is. % om de andere minima te bepalen moet je een ander x-bereik kiezen. % om de bijbehorende functiewaarde (y-as waarde) te definieren kun je % recht toe recht aan substitueren, x=xmin ymin=1/100*x^4-x^2+3/10*x-6

15

Deel 3 Handleiding Matlab voor W en BMT trimester 1.3, dagdeel 6, W H. Haitjema Faculteit Werktuigbouwkunde 2000/2001

Verwerken van meetgegevens In onderstaande opgaven wordt een aantal mogelijkheden om met Matlab meetgegevens te verwerken doorgenomen. Door de tekst goed te lezen en de voorbeeld-commandoregels uit te voeren moet men in staat zijn enkele ‘functions’ volgens de opdrachten te maken. Deze zijn later van nut wanneer ‘echte’ metingen moeten worden verwerkt.

1

Beschrijving van waarnemingen met één variabele

Iemand wil de valversnelling meten door een kogel in vacuum van 1 m hoogte te laten vallen en de tijd te meten met een stopwatch. Hij herhaalt deze meting 100 maal. Als we in werkelijkheid deze metingen doen kunnen we ze meteen intypen en aan de variabele x toekennen bijvoorbeeld voor 7 metingen: >>x=[2.15 2.12 2.13 2.23 2.28 2.12 2.29]

Een andere mogelijkheid is om de metingen eerst in een file op te slaan en deze vervolgens met commando’s zoals load of fscanf in te lezen. Dit wordt hier niet verder behandeld. Om een aantal verwerkingsmogelijkheden te illustreren zullen we geen werkelijke metingen gebruiken maar een honderdtal waarden zoals deze gemeten zouden kunnen zijn. Deze metingen xi simuleren we met behulp van de Matlabfunctie randn. Ga na wat het verschil is tussen deze functie en de functie rand met behulp van de help functie >>help rand

of >>help randn

De gesimuleerde metingen plaatsen we als volgt in array ‘x’: >>x=2.21+0.1*randn(1,100);

Bekijk het resultaat met >>plot(x,’.’);

De ‘.’ zorgt ervoor dat alleen de meetpunten wordt getekend en dat er geen lijn wordt getrokken. We zullen nu nagaan op welke manieren we deze metingen met Matlab kunnen beschrijven.

2

VERWERKEN VAN MEETGEGEVENS

2.8

meetwaarde

2.6 2.4 2.2 2 1.8 1.6 0

20

40 60 meting nummer

80

100

De statistische basisgrootheden zijn het gemiddelde, de variantie en de standaardafwijking. Het gemiddelde wordt gegeven door: x=

1 ⋅ ∑ xi n i

In Matlab wordt het gemiddelde gegeven door de functie mean >>xgem=mean(x)

De standaardafwijking wordt gegeven door

s=

n 1 ⋅ ∑ ( xi − x ) 2 = n − 1 i =1

2  n   ∑ xi 1  n 2  i =1   ∑ xi − n − 1  i =1 n  

      

Het linkerlid is de definitie; het rechterlid geeft weer hoe een computer (of een rekenmachine) deze berekent zonder dat alle afzonderlijke meetpunten in het geheugen hoeven worden opgeslagen. In Matlab wordt de standaardafwijking berekend met de functie std >>s=std(x)

Ga met de help- functie na wat het verschil is tussen >>std(x,0)


3

en >>std(x,1)

Het kwadraat van de standaardafwijking wordt variantie v genoemd: v = s2 In Matlab wordt deze berekend met de functie var: De standaardafwijking in het gemiddelde van n metingen xi wordt gegeven door sm = s/√n waarin s de standaardafwijking in een enkele meting is. OPDRACHT Maak nu een Matlab functie die van een aantal metingen n >1 het gemiddelde, de standaardafwijking en de standaardafwijking in het gemiddelde berekent. De functie heet ‘beschrijf’. De bijbehorende m-file heet dus ‘beschrijf.m’. De eerste regel van die functie moet er als volgt uitzien: function [xgem,sx,smx]=beschrijf(x)

Als invoergrootheid geldt de variabele x; uitvoergrootheden zijn ‘xgem’, ‘sx’ en ‘smx’. Bij een latere aanroep kunnen de variabelen elke naam hebben: bijvoorbeeld >>[gemidd,standafw,gstandafw]=beschrijf(cx)

zorgt ervoor dat het gemiddelde in variabele ‘gemidd’ staat etc. De variabele ‘cx’ moet wel bestaan en een passende dimensie hebben. bouw een veiligheid in voor als het aantal metingen 1 bedraagt en test de functie met >>[xg,sx,smx]=beschrijf([1,2,3])

Reken de uitkomst handmatig na Er zijn binnen Matlab en zijn statistische toolbox nog veel andere functies die eigenschappen van een reeks meetgegevens beschrijven. Het minimum van een aantal meetpunten in array x : >>min(x)

Het maximum van een aantal meettpunten in array x: >>max(x)

Het bereik van een aantal meetpunten in array x: >>range(x)

De mediaan van een aantal meetpunten in array x: >>median(x)

Het gemiddelde waarbij een percentage (hier 10%) van de extreme meetwaarden niet wordt meegerekend: >>trimmean(x,10)

De gemiddelde afwijking van punten van het gemiddelde (vergelijk met std; welke is altijd het grootst ?)

4


>>mad(x)

De meetwaarde die groter is dan 5% van de meetwaarden >>prctile(x,5)

De meetwaarde die groter is dan 95% van de meetwaarden >>prctile(x,95)

gebruik de help-functie om de betekenis van deze functies te achterhalen, bijvoorbeeld >>help mad

als je het echt niet meer weet kun je Matlab altijd nog naar het ‘waarom’ vragen >>why

Al deze functies zijn gedefinieerd voor een willekeurige verdeling van de metingen. Wanneer we meetresultaten willen koppelen aan een betrouwbaarheidsinterval is het van wezenlijk belang om iets over de verdeling van metingen te kunnen zeggen. Een snelle indruk van de verdeling van metingen wordt gegeven door een histogram; binnen Matlab wordt deze direct getekend met de functie hist, bijvoorbeeld: >>hist(x);

Dit histogram bevat 10 intervallen We kunnen ook aangeven op welke intervallen we het histogram we getekend willen zien: >>hist(x,1.5:0.05:2.5);

Hiermee wordt een histogram getekend, van x=1,5 tot 2,5 met stappen van 0,05

Een equivalent hiervoor is >>[n,y]=hist(x,1.5:0.05:2.5); >>bar(y,n);

De kans dat een meetwaarde in een bepaald interval terechtkomt wordt gevonden door de waarden in het histogram door het totaal aantal meetwaarden te delen. De lengte van een vector wordt in Matlab gegeven door de functie length: >>length(x)

Een kansfunctie in histogram-vorm wordt dus getekend met: >>bar(y,n/length(x));

Deze verdeling kunnen we vergelijken met een theoretisch kansdichtheidsfunctie. De meest toegepaste is de normale (Gaussische) verdeling. De kansdichtheidsfunctie voor een Gaussische verdeling in Matlab is normpdf(x,xgem,s) met xgem is het gemiddelde van de verdeling, s is de standaardafwijking van de verdeling en x de waarde waarvoor de kansdichtheid wordt berekend. De bekende Gauss-kromme kunnen we als volgt laten plotten: In vector ‘z’ plaatsen we 200 getallen lopend van 1,9 tot 2,5 met de functie linspace: >>z=linspace(1.9,2.5,200);


5

Voor elk element van ‘z’ bepalen we de kansdichtheidsfunctie ‘kdh’ voor een Gaussische verdelingsfunctie met een gemiddelde van 2.21 en standaardafwijking 0.1: >>kdh=normpdf(z,2.21,0.1);

De Gauss-kromme wordt nu geplot met >>plot(z,kdh)

De kansdichtheid bij meetwaarde ‘2’ wordt gegeven door >>normpdf(2,2.21,0.1)

De kans dat een meetwaarde (bijvoorbeeld) in het interval [1.98;2.02] ligt wordt nu (in benadering) gevonden door de kansdichtheid te vermenigvuldigen met de breedte van het interval. >>kans=0.04*normpdf(2,2.21,0.1)

De kans op het vinden van een waarde kleiner dan een gegeven waarde wordt gegeven door de zogeheten cumulatieve kansdichtheidsfunctie. In Matlab is deze te berekenen met de functie normcdf(x,xgem,s) waarin xgem is het gemiddelde van de verdeling, s is de standaardafwijking van de verdeling en x de waarde waarvoor de zogeheten ‘overschrijdingskans’ wordt berekend, uitgaande van een normale verdeling. De kans om een kleinere waarde dan ‘2’ te vinden bij een gemiddelde van 2,21 en een standaardafwijking van 0,1 wordt dus gegeven door: >>kans=normcdf(2,2.21,0.1)

De inverse functie hiervan is in Matlab norminv(p,xgem,s) met p de overschrijdingskans, xgem het gemiddelde en s de standaardafwijking. De waarde die met 97,5 % kans wordt overschreden wordt dus gegeven door >>norminv(0.975,2.21,0.1)

De kans om een meting in een bepaald (klein) interval aan te treffen is de kansdichtheidsfunctie maal dit interval. We gebruiken dit om de metingen te vergelijken met een normale verdeling op de volgende manier: Een histogram wordt berekend; bijvoorbeeld met 15 intervallen >>[n,y]=hist(x,15);

De kansdichtheid volgens de normale verdeling ter plaatse van y vermenigvuldigen we met de breedte van het interval en met het aantal metingen; zodat we in het histogram vergelijkbare waarden krijgen. >>z=(y(2)-y(1))*normpdf(y,mean(x),std(x))*length(x);

een histogram wordt getekend >>bar(y,n);

de grafiek wordt vastgehouden zodat een volgende grafiek binnen dezelfde assen wordt getekend. >>hold on

De Gauss-kromme wordt in zwart getekend (optie k):

6


>>plot(y,z,’k’); >>hold off

We zien nu visueel dat de metingen aardig volgens de theoretische normale manier zijn verdeeld. Eventueel is een vloeiender kromme te tekenen door de Gauss-curve bij meer punten op de x-as te berekenen. OPDRACHT Met behulp van bovenstaande gegevens breiden we nu de gemaakte functie ‘beschrijf’ uit. De eerste regel aanroep blijft function [xgem,sx,smx]=beschrijf(x)

maar nu wordt er een histogram geplot samen met de theoretische Gaussische verdeling. Het aantal te nemen intervallen N bedraagt volgens NEN 1047 voor een aantal metingen n: n < 40 N=6 40 # n # 400 N ≈ √n n > 400 N=20 Gebruik if, elseif en end om binnen de functie dit aantal intervallen te bepalen. Gebruik de functie round om het aantal intervallen op een geheel getal af te ronden. vergelijk hiermee de gesimuleerde metingen >>x=2.21+0.1*randn(1,100); >>beschrijf(x);

en >>x=2.21+0.1*(2*rand(1,100)-1.0)*sqrt(3); >>beschrijf(x);

gecombineerde metingen verdelen zich al snel als een Gauss-verdeling. Probeer het gemiddelde van twee rechthoekige verdelingen: >>x=2.21+0.1*(rand(1,100)+rand(1,100)-1.0)*sqrt(6); >>beschrijf(x);

Bij een beperkt aantal metingen moet, om tot een betrouwbaarheidsinterval voor xi of voor het gemiddelde van x te komen, nog worden gecorrigeerd voor het feit dat behalve het gemiddelde van x ook de standaardafwijking niet de ‘werkelijke’ is maar slechts een schatting daarvan. Om een bepaald betrouwbaarheidsinterval te verkrijgen moeten we de standaardafwijking vermenig vuldigen met een bepaalde factor (de ‘dekkingsfactor’, ook wel t-factor of k-factor genoemd) die afhangt van het aantal metingen en de gevraagde betrouwbaarheid. Voor een 95% betrouwbaarheid


7

(tweezijdig) en genoeg metingen is deze factor ongeveer 2; bijvoorbeeld voor de gesimuleerde metingen xi: >>k_factor=tinv(0.975, length(x)-1)

Hiervan maakt ook de functie normfit gebruik. Daarnaast geeft normfit nog een onzekerheidsinterval voor de standaardafwijking. Probeer >>[xgem,s,intervalxgem,intervals]=normfit(x,0.05)

en interpreteer het resultaat. Ga na hoeveel metingen er nodig zijn om de standaardafwijking met een relatieve onzekerheid van 5% te kunnen schatten. Hoe goed is de standaardafwijking bepaald bij 5 metingen ? OPDRACHT Gebruik tinv of normfit voor het verder aanpassen van de function beschrijf: function [xgem,sx,smx,tsx,tsmx]=beschrijf(x)

invoer x meetdata uitvoer sx standaardafwijking in x smx standaardafwijking in gemiddelde van x tsx vergrote standaardafwijking voor 95% betrouwbaarheidsinterval van x tsmx idem voor gemiddelde van x er wordt weer een histogram geplot samen met de theoretische Gaussische verdeling. De theoretische Gaussische verdeling wordt vloeiend geplot door er altijd 200 punten voor te nemen. 20

aantal waarnemingen

18 16 14 12 10 8 6 4 2 0 1.8

1.9

2

2.1

2.2 interval

2.3

2.4

2.5

2.6

8

2


Verwerken van waarnemingen met twee variabelen

Iemand meet nu de valversnelling door een kogel in vacuum van verschillende hoogten tussen de 1 en 10 m te laten vallen en de tijd te meten met een stopwatch. Er worden 100 verschillende hoogten genomen. Deze metingen xi simuleren we als volgt: Eerst genereren we 100 hoogten tussen de 1 en 10 meter met de functie linspace in de vector ‘h’: >>h=linspace(1,10,100);

Met de ‘workspace browser’ kunnen we nagaan of de vector ‘h’ de waarden bevat die we willen: klik op de knoppenbalk op de ‘workspace browser’. Je ziet nu een overzicht van de variabelen die in gebruik zijn. Dubbelklik op ‘h’: je ziet dan de inhoud in een soort spreadsheet-vorm. Je kunt de waarden zelfs wijzigen door erop te dubbelklikken. Definieer de valversnelling >>g=9.81;

De valtijd in seconden wordt gegeven door t =

2⋅h : g

>>t=sqrt(2*h/g);

In de gemeten tijd zit echter een standaardafwijking van 0,2 seconde: >>t=t+0.2*randn(1,100);

De ‘metingen’ die we hebben gedaan kunnen worden geplot: >>plot(h,t,'.')

met weer de ‘.’ optie om punten i.p.v. een lijn te krijgen. De mate waarin twee grootheden x en y samenhangen wordt gegeven door de covariantie die is gedefinieerd als n

cov( x, y ) = ∑ i =1

( xi − x) ⋅ ( yi − y ) n −1

De covariantie is positief als bij metingen (xi,yi) zowel xi als yi tegelijk groter of kleiner dan het gemiddelde zijn. De covariantie is negatief als bij metingen (xi,yi) xi groter dan het gemiddelde is terwijl yi kleiner is dan het gemiddelde en omgekeerd. Als bij een meetserie het eerste ongeveer even vaak voorkomt als het tweede gaat de covariantie naar 0; aan de sommatie zijn dan evenveel positieve als negatieve bijdragen. De Matlab functie cov geeft een matrix met op de diagonaal de varianties en daarnaast de covarianties: >>cov(h,t)

of >>cov(t,h)


9

De covariantie heeft de dimensie van x· y. Daarom is een waarde vaak weinig inzichtelijk en wordt er vaker gerekend met de correlatiecoëfficiënt r die is gedefinieerd als: cov( x, y ) sx ⋅ s y In Matlab wordt een matrix van coëfficiënten gegeven door corrcoef: rx , y =

>>cormat=corrcoef(h,t)

Het niet-diagonaal element hiervan is de hierboven gedefinieerde correlatiecoëfficiënt: >>rxy=cormat(1,2)

Deze correlatie-coefficiënt geeft aan in hoeverre twee variabelen iets met elkaar te maken hebben: bij een absolute waarde kleiner dan ±0,5 zullen de twee grootheden nauwelijks een (lineair) verband hebben. Bij veel metingen en experimenten waarbij een oorzakelijk verband het uitgangspunt is (bv. het meten van de veerconstante: kracht tegen uitrekking) zal de correlatiecoëfficiënt in de buurt van de 1 liggen en ze zegt dan weinig meer. Veel vaker gebruikt wordt het aanpassen van een aantal meetpunten aan een rechte lijn; de zogeheten lineaire kleinste-kwadraten methode. Bij deze methode wordt door de punten (x,y) een rechte lijn getrokken waarvoor geldt: y = a +b⋅ x a en b zijn de coëfficiënten die dusdanig worden berekend dat de som van de kwadraten van de afstand van elk punt tot deze lijn minimaal is. In Matlab worden deze coefficiënten berekend met de functie polyfit(x,y,n); n is de orde van de berekening: 1 is lineair, 2 is kwadratisch, etc. De uitvoer is een array met de coëfficiënten. >>coeff=polyfit(h,t,1)

De hoogste-orde coëfficiënt staat vooraan in de array ‘coeff’. Dus >>a=coeff(2) >>b=coeff(1)

De ‘y’-waarde van de berekende rechte lijn bij een gegeven ‘x’ wordt gegeven door de functie polyval(coeff,x), waarbij coeff en x de variabelen zijn. Dit kunnen we gebruiken om de meetpunten en de berekende lijn in een grafiek te plotten: De berekende tijd voor de hoogte h volgens de kleinste-kwadratenlijn noemen we tber, dus: >>tber=polyval(coeff,h);

Nu kunnen we zowel de metingen als de berekende kleinste-kwadratenlijn in een grafiek plotten; de metingen als punten en de kleinste-kwadratenlijn als een lijn. >>plot(h,t,’.’,h,tber)

10


Het verschil tussen de gemeten en de berekende punten noemen we de residuen: >>resid=t-tber;

Een maat voor hoe ‘goed’ de lijn de metingen volgt wordt gegeven door de standaardafwijking in het verschil tussen de berekende lijn en de gemeten punten die wordt gecorrigeerd voor het feit dat door twee meetpunten een rechte lijn met standaardafwijking 0 gaat: >>stda=std(resid)*sqrt((length(resid)-1)/(length(resid)-2)

We kunnen het resultaat ook bekijken met de eerder gemaakte functie ‘beschrijf’: >>beschrijf(resid)

We krijgen nu niet alle uitvoer-variabelen, maar de grafiek wordt getekend zoals we in de functie ‘beschrijf’ hebben geprogrammeerd. De coëfficiënten hebben elk ook een standaardafwijking en een correlatie, maar omdat Matlab hiervoor geen eenvoudige standaard-voorzieningen heeft worden deze hier niet verder beschouwd.

OPDRACHT maak een functie met als m-filenaam ‘lkk.m’. De eerste regel moet zijn: function [a,b,s]=lkk(x,y)

uitvoer moet zijn

a: constante in kleinste-kwadraten-aanpassing b: richtingscoëfficiënt in kleinste-kwadraten-aanpassing s: de standaardafwijking die aangeeft hoe goed de aanpassing is verder moet er een grafiek worden getekend met de meetpunten samen met de kleinste-kwadraten lijn. 8 7

tijd in seconden

6 5 4 3 2 1 0

2

4 6 valhoogte in meter

8

10

Het is ook mogelijk om hogere-orde polynomen aan de metingen aan te passen, bijvoorbeeld als 3e orde:


11

Noem de orde van het polynoom k: >>k=3;

De coëfficiënten berekenen we met polyfit >>coeff=polyfit(h,t,k);

We berekenen punten met polyval >>tber=polyval(coeff,h);

plotten de meetpunten met de berekende aanpassing >>plot(h,t,’.’,h,tber);

berekenen residuen >>resid=t-tber;

berekenen de standaardafwijking >>stda=std(resid)*sqrt((length(resid)-1)/(length(resid)-k-1)

Bij weinig of niet-uniform verdeelde meetpunten kunnen we ook standaard bijvoorbeeld 200 punten nemen om de berekende kromme te tekenen: >>hh=linspace(min(h),max(h),200); >>ttber=polyval(coeff,hh); >>plot(h,t,’.’,hh,ttber);

OPDRACHT maak een functie met als m-filenaam ‘multifit.m’. De eerste regel moet zijn: function [s,coeff]=multifit(x,y,k)

invoer

x,y: meetdata y als functie van x k: de orde van de aanpassing (>1) uitvoer moet zijn s: de standaardafwijking die aangeeft hoe goed de aanpassing is coeff: de coefficienten verder moet er een grafiek worden getekend met de meetpunten samen met de uit de coefficienten berekende lijn. Voor de berekende lijn worden altijd 100 equidistante punten genomen. Daarnaast wordt met de functie ‘beschrijf’ een grafiek getekend van de verdeling van de residuen, samen met de theoretische Gauss-verdeling. Ga nu na tot welke orde van de aanpassing de standaardafwijking duidelijk kleiner wordt. Uiteindelijk heeft het geen zin om hogere ordes te nemen als de standaardafwijking ongeveer 0,2 is. Doe nu hetzelfde voor de variabelen omgekeerd: de hoogte als functie van de tijd met een 2e orde aanpassing: >>[s,coeff]=multifit(t,h,2)

12


12

valhoogte in meter

10

8

6

4

2

0

1

2

3

4 5 tijd in seconden

6

7

8

bereken uit de waarde van coeff(1) de valversnelling. In Matlab zit een demonstratieprogramma voor polynoomfitting: polytool. Probeer: >>polytool(t,h)

Er zijn diverse mogelijkheden: veranderen van de orde van de aanpassing, uitvoer van diverse variabelen met een betrouwbaarheidsinterval, met de muis een assenkruis in de grafiek bewegen. Door met de muis in de getoonde grafiek te bewegen worden verschillende berekende waarden getoond. Het is ook mogelijk variabelen van een willekeurige functie aan metingen aan te passen. Binnen Matlab kan dit met de functie nlinfit. De gebruiker moet zelf de functie definiëren en startwaarden voor de variabelen geven. Uitgebreide behandeling valt buiten het kader van deze oefenmiddag.

Deel 4 Handleiding Matlab voor W en BMT trimester 1.3, dagdeel 6, BMT P.P.J. van den Bosch Faculteit Elektrotechniek 2000/2001

#JCARPGA?J #LEGLCCPGLE +C?QSPCKCLR ?LB !MLRPMJ %PMSN

?QGA SQ?EC MD 1'+3*',) UGRF +2*

2FC 1'+3*',) *G@P?PW @PMUQCP 3N BQD@SD @ MDV 2(,4+(-* LNCDK EHQRS NODM SGD 2(,4+(-* KHAQ@QX AQNVRDQ AD SXOHMF HM SGD

,

3+

! BNLL@MC VHMCNV

RHLTKHMJ 3GHR VHKK NODM @ VHMCNV @R RGNVM ADKNV

3GD SNNKA@Q NE SGD 2HLTKHMJ KHAQ@QX AQNVRDQ NEEDQR SGD ENKKNVHMF ED@STQDR EQNL KDES SN QHFGS

•

$REATE A NEW MODEL

.ODMR @ MDV DLOSX VNQJRGDDS VGDQD @ RHLTKHMJ LNCDK B@M AD

BNMRSQTBSDC

•

0PEN A MODEL

•

4TAY ON TOP

•

.ODMR @ CH@KNF ANW VGDQD XNT B@M RODBHEX @M DWHRSHMF LNCDK SN AD NODMDC

%NQBDR SGD KHAQ@QX AQNVRDQ VHMCNV SN @KV@XR RS@X NM SNO NM XNTQ CDRJSNO

'IND A BLOCK IN THE LIBRARY BROWSER 6GDM XNT RODBHEX @ AKNBJ M@LD HM SGD DMSQX EHDKC MDWS SN SGD HBNM BKHBJHMF SGD HBNM VHKK AQHMF XNT SN SGD KNB@SHNM NE SGHR AKNBJ

3GD 2(,4+(-* KHAQ@QX HR CHUHCDC HMSN RDUDQ@K B@SDFNQHDR .M SNO VD @KV@XR G@UD SGD 2HLTKHMJ

RS@MC@QC KHAQ@QX BNMS@HMHMF LNRS NE SGD RS@MC@QC AKNBJR 3GD B@SDFNQHDR KHRSDC ADKNV SGD

2

SIMULINK

2HLTKHMJ KHAQ@QX CDODMC NM SGD KNB@KKX HMRS@KKDC ED@STQDR (LONQS@MS ENQ "NMSQNK 3GDNQX HR NE BNTQRD "NMSQNK 2XRSDL 3NNKANW KHAQ@QX B@SDFNQX B@M AD DWO@MCDC AX BKHBJHMF SGD RHFM VGHBG VHKK SGDM RGNV SGD TMCDQKXHMF GHDQ@QBGX

S SGD KNVDRS KDUDK XNT EHMC SGD AKNBJR SG@S B@M AD CQ@FFDC SN XNTQ VNQJRGDDS

6GDM XNT RDKDBS @ AKNBJ KDES LNTRD BKHBJ HM SGD KHAQ@QX AQNVRDQ XNT B@M QD@C @ RGNQS CDRBQHOSHNM @MC UHDV HSR @OOD@Q@MBD HM XNTQ VNQJRGDDS @S SGD ANSSNL NE SGD KHAQ@QX AQNVRDQ VHMCNV RDD EHFTQD ADKNV

!MLQRPSARGLE ? 1GKSJGLI KMBCJ #LOCKS !KNBJR EQNL SGD 2HLTKHMJ KHAQ@QX B@M AD CQ@FFDC KDES LNTRD CNVM @MC CQ@F SN XNTQ VNQJRGDDS

%NQ DW@LOKD @ESDQ CQ@FFHMF SGD AKNBJ ?2HLTKHMJ"NMSHMTNTR3Q@MREDQ %BM SN XNTQ VNQJRGDDS HS VHKK KNNJ KHJD SGD EHFTQD ADKNV

SIMULINK

3

6GDM XNT CN @ ?QHFGS LNTRD BKHBJ NUDQ SGD AKNBJ HM XNTQ VNQJRGDDS @ LDMT ENQ RDSSHMF U@QHNTR

AKNBJ OQNODQSHDR VHKK ONO TO 3@JD @ KNNJ @MC SQX @ EDV NE SGDL 3GD OQNODQSHDR SG@S BNMSQNK SGD @OOD@Q@MBD NE SGD AKNBJ B@M @KRN AD ENTMC HM SGD %NQL@S LDMT !KNBJR ONRHSHNMDC NM SGD VNQJRGDDS B@M @KV@XR AD QDONRHSHNMDC AX OTSSHMF SGD BTQRNQ NUDQ SGD AKNBJ SGDM OQDRR CNVM SGD KDES LNTRD ATSSNM @MC CQ@F SGD AKNBJ 6GDM XNT LNUD SGD BTQRNQ NUDQ @ AKNBJ @MC GNKC SGDQD ENQ @ LNLDMS MNS OQDRRHMF @MX LNTRD ATSSNM RNLD HMENQL@SHNM @ANTS SG@S AKNBJ VHKK ONO TO 3GD @LNTMS NE HMENQL@SHNM SG@S XNT RDD B@M AD RDS AX SGD HSDLR HM SGD RTA LDMT

4GCU JMAI "?R? 2GNQ EQNL SGD

L@HM LDMT A@Q

#LOCK PARAMETERS

KLNRS DUDQX AKNBJ G@R HSR NVM RODBHEHB @ RDS O@Q@LDSDQR SG@S LTRS AD RODBHEHDC ADENQD XNT B@M QTM SGD RHLTK@SHNM 3N RODBHEX SGD !KNBJ O@Q@LDSDQR NODM SGD !KNBJ O@Q@LDSDQR CH@KNF ANW AX KDES BKHBJHMF NM SGD AKNBJ 3GD O@Q@LDSDQ RODBHEHB@SHNM B@M AD CNMD HM SVN V@XR

#HQDBSKX EHKK HM SGD MTLDQHB U@KTDR NE SGD O@Q@LDSDQR RDD OHBSTQD ADKNV

2ODBHEX SGD M@LDR NE U@QH@AKDR SG@S NQ VHKK AD RODBHEHDC HM SGD ,

3+

! VNQJRO@BD ADENQD XNT

QTM SGD RHLTK@SHNM RDD OHBSTQD ADKNV

2ODBHEHB GDKO NM SGD LD@MHMF NE SGD AKNBJ O@Q@LDSDQR B@M AD ENTMC AX OQDRRHMF SGD 'DKO ATSSNM HM SGD !KNBJ /@Q@LDSDQR CH@KNF ANW

$OPYING BLOCKS AND BLOCK NAMES .MBD @ AKNBJ G@R ADDM CQ@FFDC SN XNTQ VNQJRGDDS SGD D@RX V@X SN L@JD @ BNOX NE SG@S AKNBJ HR SN GNKC CNVM SGD QHFGS LNTRD ATSSNM NUDQ SGD AKNBJ @MC CQ@F HS SN SGD OK@BD VGDQD XNT V@MS SGD BNOX SN AD -NV XNT G@UD @ AKNBJ VHSG SGD R@LD AKNBJ O@Q@LDSDQR @MC OQNODQSHDR @R SGD NQHFHM@K NMD SGD NMKX SGHMF SG@S G@R BG@MFDC HR SGD AKNBJ M@LD VGHBG HR CHROK@XDC @ANUD NQ ADMD@SG DUDQX SGD

4

SIMULINK

AKNBJ 3GD AKNBJ M@LD NE SGD BNOX HR SGD M@LD NE SGD NQHFHM@K AKNBJ OKTR @ MTLADQ 8NT B@M BG@MFD SGD AKNBJ M@LD AX KDES BKHBJHMF HS (S VHKK SGDM CHROK@X @ BTQRNQ HM HS @MC HR DCHS@AKD

$ONNECTING BLOCKS &DMDQ@KKX @ RHLTKHMJ LNCDK BNMS@HMR LNQD AKNBJR CQ@FFDC EQNL SGD 2HLTKHMJ KHAQ@QX AQNVRDQ NQ BNOHDC @R DWOK@HMDC @ANUD $@BG AKNBJ B@M G@UD DHSGDQ @M HMOTS @M NTSOTS NQ ANSG (MOTSR NE AKNBJR B@M AD BNMMDBSDC SN NTSOTSR NE NSGDQ AKNBJR AX CQ@FFHMF @QQNVR /TS SGD BTQRNQ NM SGD AKNBJR HMOTS ONQS RHFM ONHMSHMF SNV@QCR SGD AKNBJ NQ SGD NTSOTS ONQS RHFM ONHMSHMF @V@X EQNL SGD AKNBJ SGD BTQRNQ VHKK BG@MFD HMSN @ BQNRR G@HQ RG@OD -NV GNKC CNVM SGD KDES LNTRD ATSSNM @MC CQ@F SNV@QCR SGD ONQS XNT V@MS SN BNMMDBS SN

8NT B@M BNMMDBS NMD NTSOTS SN RDUDQ@K HMOTSR RDUDQ@K NTSOTSR SN NMD HMOTS HR MNS ONRRHAKD 3N CN SGHR ITRS ?S@O @ BNMMDBSHNM KHMD @KQD@CX FNHMF EQNL SGD NTSOTS AX LNUHMF SGD BTQRNQ NUDQ SGD KHMD SGD BTQRNQ RG@OD VHKK MNS BG@MFD MNV Õ -NV GNKC CNVM SGD QHFGS LNTRD JDX @MC CQ@F SN SGD HMOTS ONQS XNT V@MS SN BNMMDBS SN

.MBD @ BNMMDBSHNM @QQNV G@R ADDM CQ@VM HSR QNTSHMF B@M AD BG@MFDC 3N CN SGHR KDES BKHBJ SGD BNMMDBSHNM @QQNV HS VHKK MNV AD ROKHS HMSN KHMD RDFLDMSR L@QJDC VHSG RL@KK AK@BJ RPT@QDR 8NT B@M QDONRHSHNM D@BG KHMD RDFLDMS NE SGD BNMMDBSHNM @QQNV AD LNUHMF SGD BTQRNQ NUDQ SGD KHMD RDFLDMS SGDM GNKC CNVM SGD KDES LNTRD ATSSNM @MC CQ@F SGD KHMD RDFLDMS (M BNLOKDW CH@FQ@LR SGD CDE@TKS QNTSHMF NE RHLTKHMJ LTRS NESDM AD BG@MFDC HM SGHR V@X SN FDS @ BKD@QDQ CH@FQ@L

6GDM XNT QDONRHSHNM BNMMDBSDC AKNBJR SGD BNMMDBSHNM @QQNVR VHKK ?ENKKNV SGD AKNBJ @R XNT CQ@F HS @MC RS@X @SS@BGDC

SIMULINK

5

1?TGLE ?LB MNCLGLE 1'+3*',) KMBCJQ MDV RHLTKHMJ LNCDK B@M AD R@UDC NQ @M DWHRSHMF NMD NODMDC TRHMF SGD HSDLR HM SGD

HM SGD MNQL@K VHMCNVR E@RGHNM

$GJC LDMT NE SGD L@HM LDMT A@Q NE SGD LNCDK VNQJRGDDS VHMCNV .E BNTQRD

XNT B@M @KRN TRD SGD BNQQDRONMCHMF HBNMR NM SGD SNNKA@Q 2(,4+(-* LNCDK M@LDR G@UD SGD DWSDMRHNM

LCK

6GDM XNT RS@QS , 3+ ! XNT B@M @KRN NODM @M DWHRSHMF 2(,4+(-* LNCDK AX SXOHMF SGD LNCDK EHKD M@LD VHSGNTS DWSDMRHNM HM SGD , SGD ,

3+

! O@SG NQ ,

3+

3+

! BNLL@MC VHMCNV L@JD RTQD SG@S SGD EHKD KNB@SHNM HR HM

!R VNQJHMF ENKCDQ HR SGD R@LD @R VGDQD XNTQ LNCDK EHKD HR

0SLLGLE ? QGKSJ?RGML .MBD @ 2(,4+(-* CH@FQ@L G@R ADDM ATHKC @KK AKNBJR BNMMDBSDC @MC AKNBJ O@Q@LDSDQR RODBHEHDC XNT B@M QTM SGD RHLTK@SHNM NE SGD LNCDK 'NVDUDQ ADENQD CNHMF RN XNT EHQRS G@UD SN RODBHEX SGD RHLTK@SHNM O@Q@LDSDQR

5HE 4IMULATION 1ARAMETERS

.ODM SGD 1GKSJ?RGML .?P?KCRCPQ BG?JME @MV SGQNTFG SGD 1GKSJ?RGML .?P?KCRCPQ LDMT NM SGD L@HM LDMT A@Q NE XNTQ LNCDK VNQJRGDDS VHMCNV (S VHKK KNNJ KHJD SGD EHFTQD ADKNV

%NQ MNQL@K RHLTK@SHNM SGD

1MJTCP N?EC HR SGD NMKX NMD XNT MDDC SN BNMRHCDQ

3GD 2NKUDQ O@FD @KKNVR XNT SN

• •

2DS SGD RHLTK@SHNM RS@QS @MC RSNO SHLDR "GNNRD SGD RNKUDQ @MC RODBHEX HSR O@Q@LDSDQR

6

SIMULINK •

2DKDBS NTSOTS NOSHNMR

CJMU WMS DGLB ? KMPC BCR?GJCB BCQAPGNRGML MD RFC GRCKQ ML RFC QMJTCP N?EC RFC RCVR AMKCQ DPMK +?RJ?@gQ FCJN N?ECQ $MP QGKNJC KMBCJQ HSQR SQC RFC QCRRGLEQ ?Q GL RFC DGESPC ?@MTC

5HE 4IMULATION TIME

8NT B@M BG@MFD SGD RS@QS SHLD @MC RSNO SHLD ENQ SGD RHLTK@SHNM AX DMSDQHMF MDV U@KTDR HM SGD 1R?PR RGKC @MC 1RMN RGKC EHDKCR 3GD CDE@TKS RS@QS SHLD HR RDBNMCR @MC SGD CDE@TKS RSNO SHLD HR RDBNMCR

2HLTK@SHNM SHLD @MC @BST@K BKNBJ SHLD @QD MNS SGD R@LD %NQ DW@LOKD QTMMHMF @ RHLTK@SHNM ENQ RDBNMCR VHKK TRT@KKX MNS S@JD RDBNMCR 3GD @LNTMS NE SHLD HS S@JDR SN QTM @ RHLTK@SHNM CDODMCR NM L@MX E@BSNQR HMBKTCHMF SGD LNCDKfR BNLOKDWHSX SGD RNKUDQfR RSDO RHYDR @MC SGD BNLOTSDQfR BKNBJ RODDC

4OLVERS OPTIONS 2HLTK@SHNM NE 2HLTKHMJ LNCDKR HMUNKUDR SGD MTLDQHB@K HMSDFQ@SHNM NE RDSR NE NQCHM@QX CHEEDQDMSH@K DPT@SHNMR .#$R 2HLTKHMJ OQNUHCDR @ MTLADQ NE RNKUDQR ENQ SGD RHLTK@SHNM NE RTBG DPT@SHNMR

!DB@TRD NE SGD CHUDQRHSX NE CXM@LHB RXRSDL ADG@UHNQ RNLD RNKUDQR L@X AD LNQD DEEHBHDMS SG@M

NSGDQR @S RNKUHMF @ O@QSHBTK@Q OQNAKDL 3N NAS@HM @BBTQ@SD @MC E@RS QDRTKSR S@JD B@QD VGDM BGNNRHMF SGD RNKUDQ @MC RDSSHMF O@Q@LDSDQR

8NT B@M BGNNRD ADSVDDM U@QH@AKD RSDO @MC EHWDC RSDO RNKUDQR 4?PG?@JCQRCN

QMJTCPQ

B@M LNCHEX SGDHQ

RSDO RHYDR CTQHMF SGD RHLTK@SHNM 3GDX OQNUHCD DQQNQ BNMSQNK @MC YDQN BQNRRHMF CDSDBSHNM $GVCBQRCN

QMJTCPQ S@JD SGD R@LD RSDO RHYD CTQHMF SGD RHLTK@SHNM 3GDX OQNUHCD MN DQQNQ BNMSQNK @MC CN MNS KNB@SD YDQN BQNRRHMFR %NQ @ SGNQNTFG CHRBTRRHNM NE RNKUDQR RDD

3QGLE +2*

(E XNT CN MNS BGNNRD @ RNKUDQ 2HLTKHMJ BGNNRDR NMD A@RDC NM VGDSGDQ XNTQ LNCDK G@R RS@SDR

3GD CDE@TKS RNKUDQ O@Q@LDSDQR OQNUHCD @BBTQ@SD @MC DEEHBHDMS QDRTKSR ENQ LNRS OQNAKDLR (M RNLD B@RDR GNVDUDQ STMHMF SGD O@Q@LDSDQR B@M HLOQNUD ODQENQL@MBD 8NT B@M STMD SGD RDKDBSDC RNKUDQ AX BG@MFHMF O@Q@LDSDQ U@KTDR NM SGD

1MJTCP O@MDK

4TEP 4IZES %NQ U@QH@AKD RSDO RNKUDQR XNT B@M RDS SGD L@WHLTL @MC RTFFDRSDC HMHSH@K RSDO RHYD O@Q@LDSDQR !X CDE@TKS SGDRD O@Q@LDSDQR @QD @TSNL@SHB@KKX CDSDQLHMDC HMCHB@SDC AX SGD U@KTD @TSN %NQ EHWDC RSDO RNKUDQR XNT B@M RDS SGD EHWDC RSDO RHYD 3GD CDE@TKS HR @KRN @TSN +?VGKSK QRCN QGXC 3GD +?V QRCN QGXC O@Q@LDSDQ BNMSQNKR SGD K@QFDRS SHLD RSDO SGD RNKUDQ B@M S@JD 3GD CDE@TKS HR CDSDQLHMDC EQNL SGD RS@QS @MC RSNO SHLDR

&DMDQ@KKX SGD CDE@TKS L@WHLTL RSDO RHYD HR RTEEHBHDMS (E XNT @QD BNMBDQMDC @ANTS SGD RNKUDQ

LHRRHMF RHFMHEHB@MS ADG@UHNQ BG@MFD SGD O@Q@LDSDQ SN OQDUDMS SGD RNKUDQ EQNL S@JHMF SNN K@QFD @ RSDO (E SGD SHLD RO@M NE SGD RHLTK@SHNM HR UDQX KNMF SGD CDE@TKS RSDO RHYD L@X AD SNN K@QFD ENQ SGD RNKUDQ SN EHMC SGD RNKTSHNM

KRN HE XNTQ LNCDK BNMS@HMR ODQHNCHB NQ MD@QKX ODQHNCHB ADG@UHNQ @MC

XNT JMNV SGD ODQHNC RDS SGD L@WHLTL RSDO RHYD SN RNLD EQ@BSHNM RTBG @R NE SG@S ODQHNC (M FDMDQ@K ENQ LNQD NTSOTS ONHMSR BG@MFD SGD QDEHMD E@BSNQ MNS SGD L@WHLTL RSDO RHYD 'LGRG?J QRCN QGXC !X CDE@TKS SGD RNKUDQR RDKDBS @M HMHSH@K RSDO RHYD AX DW@LHMHMF SGD CDQHU@SHUDR NE SGD RS@SDR @S SGD RS@QS SHLD (E SGD EHQRS RSDO RHYD HR SNN K@QFD SGD RNKUDQ L@X RSDO NUDQ HLONQS@MS ADG@UHNQ 3GD HMHSH@K RSDO RHYD O@Q@LDSDQ HR @

QSEECQRCB EHQRS RSDO RHYD 3GD RNKUDQ SQHDR SGHR RSDO RHYD

ATS QDCTBDR HS HE DQQNQ BQHSDQH@ @QD MNS R@SHREHDC

&RROR 5OLERANCES 3GD RNKUDQR TRD RS@MC@QC KNB@K DQQNQ BNMSQNK SDBGMHPTDR SN LNMHSNQ SGD DQQNQ @S D@BG SHLD RSDO #TQHMF D@BG SHLD RSDO SGD RNKUDQR BNLOTSD SGD RS@SD U@KTDR @S SGD DMC NE SGD RSDO @MC @KRN

SIMULINK

7

JMA?J CPPMP SGD DRSHL@SDC DQQNQ NE SGDRD RS@SD U@KTDR 3GDX SGDM BNLO@QD SGD KNB@K ?AACNR?@JC CPPMP VGHBG HR @ ETMBSHNM NE SGD QDK@SHUD SNKDQ@MBD PRMJ @MC @ARNKTSD SNKDQ@MBD ?RMJ (E SGD DQQNQ HR FQD@SDQ SG@M SGD @BBDOS@AKD DQQNQ ENQ ?LW RS@SD SGD RNKUDQ QDCTBDR SGD RSDO RHYD

CDSDQLHMD SGD DQQNQ SN SGD

@MC SQHDR @F@HM

•

0CJ?RGTC RMJCP?LAC LD@RTQDR SGD DQQNQ QDK@SHUD SN SGD RHYD NE D@BG RS@SD 3GD QDK@SHUD SNKDQ@MBD QDOQDRDMSR @ ODQBDMS@FD NE SGD RS@SDfR U@KTD 3GD CDE@TKS D LD@MR SG@S SGD BNLOTSDC

•

RS@SD VHKK AD @BBTQ@SD SN VHSGHM @QMJSRC RMJCP?LAC

HR @ SGQDRGNKC DQQNQ U@KTD 3GHR SNKDQ@MBD QDOQDRDMSR SGD @BBDOS@AKD DQQNQ

@R SGD U@KTD NE SGD LD@RTQDC RS@SD @OOQN@BGDR YDQN 3GD DQQNQ ENQ SGD HSG RS@SD

CH HR QDPTHQDC SN R@SHREX

(E XNT RODBHEX @TSN SGD CDE@TKS 2HLTKHMJ RDSR SGD @ARNKTSD SNKDQ@MBD ENQ D@BG RS@SD HMHSH@KKX SN D R SGD RHLTK@SHNM OQNFQDRRDR 2HLTKHMJ QDRDSR SGD @ARNKTSD SNKDQ@MBD ENQ D@BG RS@SD SN SGD L@WHLTL U@KTD SG@S SGD RS@SD G@R @RRTLDC SGTR E@Q SHLDR SGD QDK@SHUD SNKDQ@MBD ENQ SG@S RS@SD 3GTR HE @ RS@SD FNDR EQNL SN @MC QDKSNK HR D SGDM AX SGD DMC NE SGD RHLTK@SHNM SGD @ARSNK HR RDS SN D @KRN (E @ RS@SD FNDR EQNL SN SGDM SGD @ARSNK HR RDS SN (E SGD BNLOTSDC RDSSHMF HR MNS RTHS@AKD XNT B@M CDSDQLHMD @M @OOQNOQH@SD RDSSHMF XNTQRDKE 8NT LHFGS G@UD SN QTM @ RHLTK@SHNM LNQD SG@M NMBD SN CDSDQLHMD @M @OOQNOQH@SD U@KTD ENQ SGD @ARNKTSD SNKDQ@MBD (E SGD L@FMHSTCDR NE SGD RS@SDR U@QX VHCDKX HS LHFGS AD @OOQNOQH@SD SN RODBHEX CHEEDQDMS @ARNKTSD SNKDQ@MBD U@KTDR ENQ CHEEDQDMS RS@SDR 8NT B@M CN SGHR NM SGD (MSDFQ@SNQ AKNBJfR CH@KNF ANW

0UTPUT 0PTIONS 3GD -SRNSR MNRGMLQ @QD@ NE SGD CH@KNF ANW DM@AKDR XNT SN BNMSQNK GNV LTBG NTSOTS SGD RHLTK@SHNM FDMDQ@SDR 8NT B@M BGNNRD EQNL SGQDD ONOTO NOSHNMR

• • •

1DEHMD NTSOTS /QNCTBD @CCHSHNM@K NTSOTS /QNCTBD RODBHEHDC NTSOTS NMKX

0CDGLC MSRNSR

3GD 0CDGLC MSRNSR BGNHBD OQNUHCDR @CCHSHNM@K NTSOTS ONHMSR VGDM SGD RHLTK@SHNM

NTSOTS HR SNN BN@QRD 3GHR O@Q@LDSDQ OQNUHCDR @M HMSDFDQ MTLADQ NE NTSOTS ONHMSR ADSVDDM SHLD RSDOR ENQ DW@LOKD @ QDEHMD E@BSNQ NE OQNUHCDR NTSOTS LHCV@X ADSVDDM SGD SHLD RSDOR @R VDKK @R @S SGD RSDOR 3GD CDE@TKS QDEHMD E@BSNQ HR 3N FDS RLNNSGDQ NTSOTS HS HR LTBG E@RSDQ SN BG@MFD SGD QDEHMD E@BSNQ HMRSD@C NE QDCTBHMF SGD RSDO RHYD 6GDM SGD QDEHMD E@BSNQ HR BG@MFDC SGD RNKUDQR FDMDQ@SD @CCHSHNM@K ONHMSR AX DU@KT@SHMF @ BNMSHMTNTR DWSDMRHNM ENQLTK@ @S SGNRD ONHMSR "G@MFHMF SGD QDEHMD E@BSNQ CNDR MNS BG@MFD SGD RSDOR TRDC AX SGD RNKUDQ 3GD QDEHMD E@BSNQ @OOKHDR SN U@QH@AKD RSDO RNKUDQR @MC HR LNRS TRDETK VGDM TRHMF NCD 3GD NCD RNKUDQ HR B@O@AKD NE S@JHMF K@QFD RSDOR VGDM FQ@OGHMF RHLTK@SHNM NTSOTS XNT L@X EHMC SG@S NTSOTS EQNL SGHR RNKUDQ HR MNS RTEEHBHDMSKX RLNNSG (E SGHR HR SGD B@RD QTM SGD RHLTK@SHNM @F@HM VHSG @ K@QFDQ QDEHMD E@BSNQ

U@KTD NE RGNTKC OQNUHCD LTBG RLNNSGDQ QDRTKSR

.PMBSAC ?BBGRGML?J MSRNSR

3GD .PMBSAC ?BBGRGML?J MSRNSR BGNHBD DM@AKDR XNT SN RODBHEX CHQDBSKX

SGNRD @CCHSHNM@K SHLDR @S VGHBG SGD RNKUDQ FDMDQ@SDR NTSOTS 6GDM XNT RDKDBS SGHR NOSHNM 2HLTKHMJ CHROK@XR @M

-SNSR 2GKCQ EHDKC NM SGD 1MJTCP

O@FD $MSDQ @ ,

3+

! DWOQDRRHNM HM SGHR EHDKC SG@S

DU@KT@SDR SN @M @CCHSHNM@K SHLD NQ @ UDBSNQ NE @CCHSHNM@K SHLDR 3GD @CCHSHNM@K NTSOTS HR OQNCTBDC TRHMF @ BNMSHMTNTR DWSDMRHNM ENQLTK@ @S SGD @CCHSHNM@K SHLDR 4MKHJD SGD QDEHMD E@BSNQ SGHR NOSHNM BG@MFDR SGD RHLTK@SHNM RSDO RHYD RN SG@S SHLD RSDOR BNHMBHCD VHSG SGD SHLDR SG@S XNT G@UD RODBHEHDC ENQ @CCHSHNM@K NTSOTS .PMBSAC QNCAGDGCB MSRNSR MLJW

3GD .PMBSAC QNCAGDGCB MSRNSR MLJW BGNHBD OQNUHCDR RHLTK@SHNM

NTSOTS MLJW @S SGD RODBHEHDC NTSOTS SHLDR 3GHR NOSHNM BG@MFDR SGD RHLTK@SHNM RSDO RHYD RN SG@S SHLD RSDOR BNHMBHCD VHSG SGD SHLDR SG@S XNT G@UD RODBHEHDC ENQ OQNCTBHMF NTSOTS 3GHR BGNHBD HR TRDETK VGDM BNLO@QHMF CHEEDQDMS RHLTK@SHNMR SN DMRTQD SG@S SGD RHLTK@SHNMR OQNCTBD NTSOTS @S SGD R@LD SHLDR !MKN?PGLE -SRNSR MNRGMLQ

R@LOKD RHLTK@SHNM FDMDQ@SDR NTSOTS @S SGDRD SHLDR

"GNNRHMF

0CDGLC MSRNSR @MC RODBHEXHMF @ QDEHMD E@BSNQ NE FDMDQ@SDR NTSOTS @S SGDRD SHLDR

8

SIMULINK

"GNNRHMF SGD

.PMBSAC ?BBGRGML?J MSRNSR NOSHNM @MC RODBHEXHMF :< FDMDQ@SDR NTSOTS @S SGDRD

SHLDR @MC ODQG@OR @S @CCHSHNM@K SHLDR CDODMCHMF NM SGD RSDO RHYD BGNRDM AX SGD U@QH@AKD RSDO RNKUDQ "GNNRHMF SGD

.PMBSAC 1NCAGDGCB -SRNSR -LJW

NOSHNM @MC RODBHEXHMF :< FDMDQ@SDR NTSOTS @S

SGDRD SHLDR

4TARTING AND STOPPING THE SIMULATION

2S@QS @ RHLTK@SHNM SGQNTFG SGD LDMT HSDL

1GKSJ?RGML 1R?PR 2SNO @ RHLTK@SHNM SGQNTFG SGD LDMT

HSDL 1GKSJ?RGML 1RMN 8NT B@M @KRN TRD SGD BNQQDRONMCHMF SNNKA@Q HBNMR

% Instructie voorbeeld: Regeltheorie % % % % % % % % % % % % % % % % % % % % % % % % % % % % %

In 1958 is door Westheimer de beweging van het oog onderzocht. Hij kwam uiteindelijk % met een tweede orde differentiaal vergelijking die de oogbeweging beschrijft. Zijn model ging uit van een open loop systeem, er was dus geen terugkoppeling via de hersenen. De differentiaal vergelijking die hij opstelde was:

% % % % % %

Opdracht voor student: Voer nu de overdrachtsfunctie van bovenstaande diff. vgl. in Matlab in. Matlab heeft daarvoor een aantal mogelijkheden: met coefficienten van teller en noemer: sys=tf([1 3],[1 2 1]); met nulpunten, polen en versterking K: sys=zpk([-3],[-1 -1], 1); gebruik help(tf) en/of help(zpk) voor meer informatie.

.. . A2*theta + A1*theta + A0*theta = f(t) Hierbij is: A2= massatraagheid van de oogbal A1= viscose demping A0= veerconstante van de spieren. f(t)= kracht uitgeoefend op het systeem theta= ooghoek in het horizontale vlak Wanneer we bovenstaande diff. vgl. in een andere vorm schrijven, en we veronderstellen voor de uitgeoefende kracht f(t) een constante stapfunctie, krijgen we: .. . theta + 2*d*wn*theta + wn^2*theta = K Hierbij is: wn=sqrt(A0/A2) = natuurlijke resonantie frequentie d=A1/(2*sqrt(A0*A2) = dempingsconstante van het systeem K= stap functie voor de kracht op de oogbal. Uit metingen van Westheimer bleek: d=0.7 en wn=120 rad/s

% **************** aanpassing door student **************** % ********************************************************** % % % % % %

Veel eigenschappen van een dynamisch systeem zijn af te leiden uit responsie aan de uitgang van systeem op een stapvormig ingangssignaal. Deze zgn. stapresponsie kunnen we karakteriseren m.b.v. de volgende kenmerken: stijgtijd, settling tijd, doorschot en piektijd (zie Franklin blz 126) In matlab kan de stapresponsie van een systeem worden berekend m.b.v. de functie 'step' (zie matlab help).

% % % % % % % % %

Opdracht voor student: Bepaal de stapresponsie van het systeem met de overdrachtsfunctie zoals hierboven ingevuld. Analyseer deze stapresponsies als volgt: Door in het figuurvenster met de rechter muisknop te klikken verschijnen een aantal menu opties, kies hiervande optie 'Characteristics'. Binnen het sub-menu dat nu verschijnt kun je een of meerdere stapresponsie-eigenschappen aanvinken, deze worden dan in de plot gevisualiseerd door bolletjes. Door vervolgens de linkermuisknop op een bolletje ingedrukt te houden, verschijnen de numerieke waarden.

% **************** aanpassing door student ****************

% ******************************************************************************* % % % % % % %

Zoals bekend kunnen we ook veel zeggen over het gedrag van een dynamisch systeem wanneer we een overzicht hebben van de polen en nulpunten van de overdrachts-functie van het systeem. Bestudeer de Matlab help van de functie 'pzmap', en gebruik deze functie om van de hierboven gegeven overdrachtsfunctie van de oogbeweging een pool nulpuntenoverzicht te maken. Druk de waarden van de polen en nulpunten tevens af in Matlab's commandovenster.

% **************** aanpassing door student **************** % *******************************************************************************

2 % % % % % % % % % % % % % % %

INSTRUCTIE Om nu de relatie tussen ligging van de polen en het tijddomein gedrag van een dynamisch systeem te bekijken zullen we de polen gaan verplaatsten, en vervolgens van het nieuwe systeem een stapresponsie analyseren. Maak nu eerste van de oorspronkelijke systeempolen het reele deel 5x zo groot maak een stapresponsie en analyseer deze. Maak nu van de oorspronkelijke systeempolen het imaginaire deel 5x zo groot Analyseer ook dit systeem aan de hand van een de stapresponsie. Zorg er in beide gevallen voor dat de statische versterking van het systeem gelijk blijft zodat de stapresponsies alle drie naar dezelfde eindwaarde toegaan. Plot de drie stapresponsies in een grafiek. Plot het pool-nulpunten overzicht in een grafiek. Gebruik voor manipulatie van complexe getallen de functies 'real', 'imag', 'complex' en 'conj'.

% **************** aanpassing door student **************** % ******************************************************************************* % % % % %

Om de relatie tussen ligging van systeempolen en frequentiegedrag te bekijken kunnen we van elke overdrachtsfunctie zoals hierboven bepaald, de bodediagrammen door Matlab laten uit rekenen. Gebruik hiervoor het commando 'bode'. Laat de bode-% diagrammen van de drie systemen in een grafiek zien (gebruik hold on)

% **************** aanpassing door student **************** % *******************************************************************************

Handleiding Matlab voor W en BMT

Recommend Documents