Kapitel 4. Programmet MATLAB MATLAB (namnet h¨arlett ur MATrix LABoratory) ¨ar ett matematikprogram baserat p˚ a matrisalgebra, som blivit mycket anv¨ant f¨or fysikaliska och tekniska till¨ampningar. Den ursprungliga versionen, som skrevs p˚ a FORTRAN av Cleve Moler i University of New Mexico, saknade grafik, men bildade ¨and˚ a stommen f¨or senare versioner av MATLAB skrivna p˚ a C-spr˚ aket. Fo¨r att l¨ara sig anv¨anda MATLAB, kan man utom handb¨ockerna som levereras med programmet, ocks˚ a anv¨anda b¨ocker, som beskriver fysikaliska och tekniska till¨ampningar av MATLAB. En kort introduktion, MATLAB Primer som skrivits av Kermit Sigmon finns p˚ a n¨atet, se t.ex. http://www.internationaleconomics.net/documents/matlab.pdf, den nyaste versionen kan k¨opas i bokhandeln). En utfo¨rligare beskrivning av MATLAB finner man i The MATLAB Handbook skriven av Eva P¨art-Enander, Anders Sj¨oberg, Bo Melin och Pernilla Isaksson i Uppsala universitet. I Uppsala har ocks˚ a skrivits en svensk handledning (se http://www.it.uu.se/edu/bookshelf/). Det finns versioner av MATLAB f¨or m˚ anga slags datorer. De f¨oljande exemplen har testats med Windowsversionen, som ¨ar ganska l¨att att anv¨anda, t.ex. via mikron¨atet (TOUKO). MATLAB (v. 6.5) finns dessutom p˚ a flere av Unix-datorerna i HU (rock, soul, heavy, etc). Introduktion till vetenskapliga ber¨akningar I, Tom Sundius 2009
JJ J I II ×
1
Det finns ocks˚ a en studentversion av MATLAB som man kan k¨opa f¨or ca 70 =C . MATLAB skiljer sig fr˚ an de symboliska programmen d¨arigenom, att den i allm¨anhet ger numeriska approximationer f¨or de r¨akneoperationer man ber den utf¨ora. I de nyaste versionerna av MATLAB (6 och h¨ogre) ing˚ ar dock ¨aven symboliska verktyg. MATLAB startas p˚ a de flesta datorer genom att man skriver matlab (f¨or att undvika att den ¨oppnar m˚ anga grafiska f¨onster p˚ a en unix–dator, vilket tar tid, kan man skriva matlab -nojvm). Programmet startar d˚ a upp och ger prompten >>. Om man sj¨alv vill ¨ova sig i att anv¨anda MATLAB, kan man pr¨ova kommandot intro, eller demo, som demonstrerar vilka slags problem som MATLAB klarar. Om man vill ha hj¨alp med n˚ agot kommando i MATLAB, skriver man help kommando. F¨or att avsluta en session med MATLAB anv¨ander man kommandot quit eller exit. Genom att trycka p˚ a upp˚ at-pilen f˚ ar man p˚ a nytt ett tidigare kommando, med v¨ansterpil flyttar man sig till v¨anster ett steg, och delete (backspace) tar bort ett tecken. Om man tycker att n˚ agon r¨akneoperation tar f¨or l˚ ang tid, kan man avbryta den med Ctrl-C.
Introduktion till vetenskapliga ber¨akningar I, Tom Sundius 2009
JJ J I II ×
2
4.1. Grunderna till MATLAB
MATLAB har en inbyggd hj¨alpfunktion, som ger information om olika kommandon, som t.ex. >> help pi PI 3.1415926535897.... PI = 4*atan(1) = imag(log(-1)) = 3.1415926535897.... Med kommandot lookfor kan man so¨ka upp s˚ adana kommandofiler, som inneh˚ aller en given textstr¨ang: >> lookfor root matlabro.m: % MATLABROOT Root directory of MATLAB installation. SQRT Square root. SQRTM Matrix square root. ROOTS Find polynomial roots. CPLXROOT Riemann surface for the n-th root.
Introduktion till vetenskapliga ber¨akningar I, Tom Sundius 2009
JJ J I II ×
3
MATLAB kan anv¨andas som en vanlig kalkylator, i likhet med Mathematica: >> 4711 + 1174 ans = 5885 >> 3^4, 5*(4+12) ans = 81 ans = 80 Resultatet anges som ans, om man inte s¨atter ut n˚ agot variabelnamn. Normalt g¨or man detta, och resultatet av en r¨akneoperation kan d˚ a se ut p˚ a detta s¨att: >> x =16 x = 16 >> y = 3*x y = 48 Introduktion till vetenskapliga ber¨akningar I, Tom Sundius 2009
JJ J I II ×
4
a Observera anv¨andningen av parenteser: Om man skriver a/b +c l¨aser MATLAB detta som + c, medan b a a/(b+c) uppfattas som . Detta beror p˚ a r¨akneoperationernas olika prioritet. b+c Fo¨r division finns det tv˚ a tecken, vanlig division (/) och v¨ansterdivision (\). 2/5 betyder detsamma som 5\2. I matrisr¨akning har v¨ansterdivisionen en speciell betydelse, som vi skall se senare. Normalt skriver MATLAB ut resultaten med 5 siffrors noggrannhet. Vill man ha mera, kan anv¨anda kommandot format long (eller format long e, som ger exponentformat). Man kan ocks˚ a behandla komplexa tal i MATLAB. Variablerna i och j anv¨ands f¨or att beteckna den imagin¨ara enheten (om de inte definierats p˚ a annat s¨att), och man kan d¨arf¨or utf¨ora r¨akneoperationer med dem: >> z=1+2i z = 1.0000 + 2.0000i >> conj(z) ans = 1.0000 - 2.0000i Introduktion till vetenskapliga ber¨akningar I, Tom Sundius 2009
JJ J I II ×
5
>> z^2 - z ans = -4.0000 + 2.0000i >> imag(z^2-z) ans = 2.0000 Observera, att conj(z) h¨ar betecknar den komplexa konjugaten av z , dvs z . Argumentet av ett komplext tal (dvs fasvinkeln i det komplexa planet) anges med angle, t.ex. >> arg=angle(z) arg = 1.1071 MATLAB k¨anner till alla de vanliga element¨ara funktionerna: >> cos(45) ans = 0.5253 Introduktion till vetenskapliga ber¨akningar I, Tom Sundius 2009
JJ J I II ×
6
>> pi ans = 3.1416 >> cos (pi/4) ans = 0.7071 sin(pi) ¨ar inte exakt 0, eftersom pi ¨ar en approximation f¨or π : >> sin(pi) ans = 1.2246e-016 Observera, att de trigonometriska funktionerna anges i radianer, men det ¨ar l¨att att konvertera dem till grader, eftersom MATLAB k¨anner till π (som vi sett). Vi skall ¨annu se p˚ a n˚ agra andra exempel. Den naturliga logaritmfunktionen betecknas med log, och den briggska med log10:
Introduktion till vetenskapliga ber¨akningar I, Tom Sundius 2009
JJ J I II ×
7
>> x=10; >> log(x) ans = 2.3026 >> log10(x) ans = 1 Nepers tal e f˚ ar man genom att r¨akna ut exponentialfunktionens v¨arde f¨or x = 1: >> e=exp(1) e = 2.7183 I princip k¨anner MATLAB till bara en matematisk storhet, en rektangul¨ar matris. Ofta tolkas dock matriser med en rad eller en kolumn som vektorer och matriser med en rad och en kolumn som skal¨arer. En kolumnvektor och en radvektor anges p˚ a olika s¨att i MATLAB:
Introduktion till vetenskapliga ber¨akningar I, Tom Sundius 2009
JJ J I II ×
8
>> vcol = [1; 2; 3; 4] vcol = 1 2 3 4 >> vrow = [ 5 6 7 8] vrow = 5
6
7
8
Man kan operera p˚ a en vektor (elementvis) med en funktion, t.ex. kvadratroten:
Introduktion till vetenskapliga ber¨akningar I, Tom Sundius 2009
JJ J I II ×
9
>> sqrt(vcol) ans = 1.0000 1.4142 1.7321 2.0000 Om man vill se vilka variabler som hittills definierats, anv¨ander man kommandot who: >> who Your variables are: a ans
b vcol
vrow x
y
Om man vill stryka variabeldefinitionerna, anv¨ander man kommandot clear. Man kan ocks˚ a alstra en vektor med kolonoperatorn: >> v = 0:8 Introduktion till vetenskapliga ber¨akningar I, Tom Sundius 2009
JJ J I II ×
10
v = 0
1
2
3
4
5
6
7
8
>> v2=0:0.5:2 v2 = 0
0.5000
1.0000
1.5000
2.0000
Som vi ser, uppfattas argumentet mellan de tv˚ a gr¨anserna som en tillv¨axtterm. Ifall den fattas, antas tillskottet vara 1. Man kan ocks˚ a utf¨ora operationer p˚ a en vektor: >> 2.^v ans = 1
2
4
8
16
32
64
128
256
Observera punkten framf¨or operatorn ^, som anger att operationen utf¨ors elementvis. Introduktion till vetenskapliga ber¨akningar I, Tom Sundius 2009
JJ J I II ×
11
Fo¨r att ber¨akna vektorprodukten av tv˚ a vektorer, kan man anv¨anda funktionen cross. Produkten av de tv˚ a vektorerna i + j och j + k (i = [100] etc.) kan s˚ alunda ber¨aknas som f¨oljer: >> a = [1 1 0]; >> b = [0 1 1]; >> c = cross(a,b) c = 1 -1 1 Matriser kan konstrueras p˚ a olika s¨att. Dels kan man konstruera dem direkt fr˚ an en elementlista, dels kan man alstra dem med hj¨alp av MATLAB-kommandon, och dels kan de l¨asas in fr˚ an filer. En matris kan konstrueras genom att r¨akna upp element f¨or element: >> B(1,1)=1; B(1,2)=2; >> B(2,1)=3; B(2,2)=4 B = 1 2 3 4 men vanligen anv¨ands en annan metod. Introduktion till vetenskapliga ber¨akningar I, Tom Sundius 2009
JJ J I II ×
12
En fyrradig matris kan man t.ex. l¨asa in genom att ange elementen radvis, ˚ atskilda av semikolon, och omge hela uttrycket med klamrar: >> A = [ 1 2 3 4 ; 2 3 4 1 ; 3 4 1 2 ; 4 1 2 3 ] A = 1 2 3 4 2 3 4 1 3 4 1 2 4 1 2 3 (om man inte vill att MATLAB skall skriva ut matrisen, kan man tillfoga ett semikolon (;) efter kommandot). En radvektor kan anges p˚ a liknande s¨att, elementen kan vara matematiska uttryck (observera, att pi betecknar symbolen π ): >> x = [-0.5 , 0.5 , 1, pi ] x = -0.5000
0.5000
1.0000
3.1416
Introduktion till vetenskapliga ber¨akningar I, Tom Sundius 2009
JJ J I II ×
13
Emedan x ¨ar en radvektor, kan man inte utf¨ora multiplikationen Ax, vilket observeras av MATLAB: >> b = A*x ??? Error using ==> * Inner matrix dimensions must agree. D¨aremot fungerar multiplikationen xA: >> b=x*A b = 16.0664
7.6416
7.7832
9.9248
Genom att transponera x till en kolumnvektor med apostrofbeteckningen >> x = x’ x = -0.5000 0.5000 1.0000 3.1416 kan man utf¨ora ber¨akningen: Introduktion till vetenskapliga ber¨akningar I, Tom Sundius 2009
JJ J I II ×
14
>> b = A*x b = 16.0664 7.6416 7.7832 9.9248 L˚ at oss definiera tv˚ a matriser A och B : >> A = [1 2 3; 4 5 6; 7 8 9] A = 1 2 3 4 5 6 7 8 9 >> B = [1 2 3; 4 5 6] B = 1 2 3 4 5 6 Introduktion till vetenskapliga ber¨akningar I, Tom Sundius 2009
JJ J I II ×
15
Vi kan l¨att bilda matrisprodukten BA: >> B*A ans = 30 66
36 81
42 96
men f¨ors¨oker vi r¨akna ut AB , lyckas det inte: >> A*B ??? Error using ==> * Inner matrix dimensions must agree. Detta beror p˚ a, att raderna i matrisen A har tre element, medan kolumnerna i B har endast tv˚ a element. Matrismultiplikationen kan allts˚ a inte utf¨oras.
Introduktion till vetenskapliga ber¨akningar I, Tom Sundius 2009
JJ J I II ×
16
Determinanten av matrisen A ¨ar noll, medan determinanten av B inte existerar: >> det(A) ans = 0 >> det(B) ??? Error using ==> det Matrix must be square. Detta leder till att matrisen A blir singul¨ar och saknar invers, medan B definitionsm¨assigt inte har n˚ agon invers eftersom den inte ¨ar en kvadratisk matris. Vi ser detta genom att f¨ors¨oka ber¨akna inverserna med MATLAB:
Introduktion till vetenskapliga ber¨akningar I, Tom Sundius 2009
JJ J I II ×
17
>> inv(A) Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.937385e-018 ans = 1.0e+016 * 0.3152 -0.6304 -0.6304 1.2609 0.3152 -0.6304
0.3152 -0.6304 0.3152
>> inv(B) ??? Error using ==> inv Matrix must be square. Om man obetydligt f¨or¨andrar n˚ agot av elementen i A, f˚ ar vi en matris som uppf¨or sig b¨attre: >> A2 = [1 2 3; 4 5 6; 7 8 9.01]; >> det(A2) ans = -0.0300 Introduktion till vetenskapliga ber¨akningar I, Tom Sundius 2009
JJ J I II ×
18
Determinanten ¨ar fortfarande ganska liten, men inversen existerar, ¨aven om konditionstalet rcond ¨ar litet: >> rcond(A2) ans = 9.2695e-005 >> inv(A2) ans = 98.3333 -199.3333 100.0000 -198.6667 399.6667 -200.0000 100.0000 -200.0000 100.0000 Ekvationssystemet Ax = b, d¨ar b ¨ar kolumnvektorn [456]0 kan dock lo¨sas med hj¨alp av MATLAB p˚ a f¨oljande s¨att: >> b = [4; 5; 6]; >> A\b
Introduktion till vetenskapliga ber¨akningar I, Tom Sundius 2009
JJ J I II ×
19
ans = -0.3333 -2.3333 3.0000 (observera ”v¨ansterdivisionen” \). Om vi byter ut matrisen A mot den modifierade matrisen A2, f˚ as b = [-3.3333, 3.6667, 0.0000], en r¨att stor skillnad. Om ett line¨art ekvationssystem Ax = b inneh˚ aller flere ekvationer ¨an obekanta, s¨ages det vara ¨overdetera man anpassar kurvor till experimentella data. Ekvationssystemet brukar d˚ a minerat. Detta intr¨affar ofta d˚ l¨osas med Gauss’ minsta kvadratmetod. MATLAB anv¨ander denna metod automatiskt d˚ a man anv¨ander operatorn \. Nedan visas ett exempel h¨arp˚ a: >> A=[1 2; 3 4; 5 6]; >> b=[7; 8; 9]; >> A\b ans = -6.0000 6.5000 Introduktion till vetenskapliga ber¨akningar I, Tom Sundius 2009
JJ J I II ×
20