Matlab∗, een mini-inleiding, P. de Groen
“Matlab” is een interactieve taal, ontworpen door Cleve Moler, die is ontstaan uit een demonstratieproject, waarmee studenten op een eenvoudige wijze zouden kunnen experimenteren met rekenmethoden in de lineaire algebra met gebruik van de standaardroutines uit LINPACK (pakket voor het oplossen van lineaire stelsels vergelijkingen) en EISPACK (pakket voor het eigenwaardeprobleem) zonder zelf de basisalgoritmen gedetailleerd te moeten implementeren.
Het ontwerp was zo succesvol, dat Moler een onderneming oprichtte (MathWorks) die het ont- werp heeft uitgebouwd tot een zeer krachtige programmeeromgeving met een uitgebreide verza- meling numerieke en grafische hulpmiddelen voor het oplossen (en simuleren) van problemen en het grafisch voorstellen van de oplossing.
De basis-datastructuur is de matrix. Met de instructie “p=5; q=7; A = rand(p,q)” wordt een re¨ele matrix gecre¨eerd met 5 rijen en 7 kolommen (in IR5×7 dus) bestaande uit uniform op [0 , 1] verdeelde random getallen. Een matrix met ´e´en kolom is een kolomvector, een matrix met ´e´en rij is een rijvector en een enkele “real” staat gelijk met een 1×1–matrix; de typen real en vector zijn dus geen aparte datatypen. Er wordt gerekend met de standaard IEEE 64-bits reals. Bewerkingen met matrices en vectoren volgen de standaard regels van de lineaire algebra.
Een matrix kun je altijd vermenigvuldigen met een scalair (= 1×1–matrix), twee matrices van gelijke afmetingen kun je optellen en aftrekken en als A∈ IRp×qen B∈ IRq×r, dan is het product de matrix AB ∈ IRp×r. Het accent wordt gebruikt voor transpositie, dus als A ∈ IRp×q dan A0 ∈ IRq×p. Het systeem rekent even gemakkelijk met complexe matrices volgens de geldende regels uit de lineaire algebra. Een complex getal wordt altijd voorgesteld in de vorm re¨eel deel + i× imaginair deel. Het accent geeft dan transpositie plus complex toevoegen aan.
Voorbeeld(>> is de matlab-prompt):
>> x=[1+i,1-i]
x =
1.0000 + 1.0000i 1.0000 - 1.0000i
>> x’
ans =
1.0000 - 1.0000i 1.0000 + 1.0000i
>> x’*x ans =
2.0000 0 - 2.0000i
0 + 2.0000i 2.0000
>> x*x’
ans = 4
>>
Er zijn allerhande mogelijkheden om deelmatrices te selecteren, b.v. als A ∈ Cp×q, dan is real(A)∈ IRp×q het re¨ele deel en imag(A)∈ IRp×q het imaginaire deel. A(:,k) is de k-de kolom (mits 1≤ k ≤ q) en A(1:3,2:2:q) is een matrix bestaande uit de elementen met even kolomindex uit de eerste 3 rijen van A. De opdracht x = A\b lost het stelsel vergelijkingen Ax = b op met de optimale methode, dus Gausseliminatie met rijverwisseling als A vierkant en goed-geconditioneerd is, en met een QR-ontbinding of een Singuliere-waardenontbinding als A slecht geconditioneerd of niet vierkant is. De dimensies van b en A moeten natuurlijk compatibel zijn. Het hele scala van standaard matrix- en vectorroutines, zoals FFT, QR, LU, Choleski, SVD en de berekening van eigenwaarden/eigenvectoren is beschikbaar.
Het systeem is op het eerste gezicht “commandline-oriented” en interpreterend, maar er zijn uitgebreide mogelijkheden tot het oproepen van routines (m-files zonder parametersubstitutie en
∗Datum: 5 okt 1998
1
2 zonder locale variabelen) en functies (met call by value invoerparameters en uitvoerparameters).
De standaard controlestructuren zijn (zie help if, help for, etc)
if <boolean> , <statements> else <statements> end for <variable> = <range> , <statements> end
while <boolean> , <statements> end
Een “ , ” (komma) en een “ ; ” (puntkomma) en een <eoln> (end-of-line) worden gebruikt om statements van elkaar te scheiden. Een opdracht in het command-window wordt ge¨ınterpreteerd en onmiddellijk uitgevoerd en het resultaat wordt in het commandwindow afgedrukt, tenzij de statement wordt afgesloten met een “ ; ”.
Omdat iedere programmaregel ge¨ınterpreteerd wordt, is de berekening van het inproduct for k = 1:100 , s = s + a(k) ∗ b(k) ; end
van de kolomvectoren a en b van lengte n, veel trager dan de opdracht s = a’ * b . Ook de statement [m,i]= max(abs(A(1:n,k))) is veel sneller dan het equivalente
m=0; for j = 1:n, if abs(A(j,k)) > m, i = j; m = abs(A(j,k)); end, end . Een declaratie van een matrix of vector is in principe niet nodig. Een datastruktuur wordt automatisch gecre¨eerd en/of uitgebreid, indien nodig. B.v. in
for k = 1:100, c(k) = a(k) ∗ b(k) ; end
wordt de vector c in iedere slag met ´e´en element verlengd (en dus worden alle reeds aanwezige elementen in iedere slag gecopieerd naar een nieuwe plaats in het geheugen). Veel sneller zal de volgende variant zijn, waarin vooraf de gehele vector ineens wordt gecre¨eerd:
c = zeros(100,1); for k = 1:100, c(k) = a (k) ∗ b(k) ; end Nog sneller kunnen we dit doen met een “elementsgewijze” vermenigvuldiging1 “ .∗ ”:
c = a .* b .
Let wel op, met “c = a .* b” wordt een nieuwe matrix van de gevraagde afmetingen gecre¨eerd, terwijl met “c(:,k) = a .* b” de k-de kolom van de bestaande matrix c wordt overschreven met nieuwe waarde. In dit laatste geval moet deze matrix c wel de goede kolomlengte hebben;
als er reeds k of meer kolommen zijn wordt de k-de overschreven en als er minder zijn, wordt de matrix uitgebreid met het benodigde aantal kolommen (met nullen).
Voor iedere functie maak je een aparte m-file met naam “<functienaam>.m” en eerste regel
“[a,b] = function <functienaam>(x,y,z)”, waarbij “x, y, z” de invoerparameters (call by value) en “a, b” de uitvoerparameters zijn. Het aantal actuele in- en uitvoerparameters mag kleiner zijn dan het aantal formele. Binnen de functie is het aantal actuele in- en uitvoerpara- meters op te vragen met “nargin” resp. “nargout”. De nietgespecificeerde invoerparameters kun je dan een default-waarde geven. De dimensies van een variabele zijn op te vragen met
“size”. Binnen een functie zijn alle variabelen locaal en kan er geen beroep worden gedaan op globaal aanwezige variabelen (behalve “pi” en “i”), terwijl een bij oproep van een m-file routine de letterlijke inhoud van de file ter plaatse wordt gesubstitueerd en uitgevoerd.
Het programma is voorzien van uitgebreide 2D en 3D grafische hulpmiddelen (met interac- tieve mogelijkheden, zie “help plot” en “help surf”) en goede documentatie via een “help
<commandnaam>” op de command line en via een help-window. Er is een uitgebreide verza- meling “toolboxen” voor ijle matrices, signaalverwerking, niet-lineaire vergelijkingen, spline- interpolatie, simulaties, wavelets, gewone en parti¨ele differentiaalvergelijkingen en er is zelfs een interface met Maple.
N.B. “i”, “eps” en “ans” zijn standaard variabelen, de imaginaire eenheid√
−1 , de standaard machineprecisie resp. het laatst berekende en in het command window afgedrukte anonyme antwoord. Je kan problemen verwachten als je deze herdefinieert.
1 De vermenigvuldigingsoperator “*”, de delingsoperator “/” en de machtsverheffing “^” volgen de gebruike- lijke rekenregels voor matrices. De elementsgewijze opratoren “.*”, “./” en “.^” voeren de bedoelde operatie uit op elementen met gelijke indices, mits beide operanden (matrices) gelijke afmetingen hebben. Als “A” een matrix is, dan is “A*A” of “A^2” het product van “A” met zichzelf en dit product is alleen gedefinieerd als “A”
vierkant is; “A.*A” of “A.^2” is een matrix van gelijke afmetingen als “A” met als elementen de kwadraten van de overeenkomstige elementen van “A”.
