MATLAB-perusteet

Ville-Pekka Eronen, Jarkko Peltomäki
Matematiikan ja tilastotieteen laitos, Turun yliopisto
Viimeksi päivitetty: 31.8.2017

Tähän dokumenttiin on koottu johdanto MATLAB-ohjelmiston (jäljempänä Matlab) käyttöön. Materiaalin on tarkoitus antaa riittävät perusvalmiudet ohjelmiston käyttöön Tero Aittokallion kurssille Tieteellinen laskenta. Materiaali ei ole johdatus ohjelmointiin tai ongelmanratkaisuun Matlabilla, vaan ainoastaan ohjelmiston käyttöön. Dokumentin materiaali on pääasiassa tarkoitettu itseopiskeluun – itse tekeminen ja kokeilu ovat parhaita tapoja oppia ohjelmistojen käyttöä ja ohjelmointia.

Contents

Matlab ja sen historia lyhyesti

Matlab eli matrix laboratory on kaupallinen ohjelmisto, joka on nimensä mukaisesti matriisipohjainen numeeriseen laskentaan tarkoitettu ohjelmisto. Sen kehitti 1970-luvun lopulla Cleve Moler, joka halusi opiskelijoilleen pääsyn Fortranilla ohjelmoituihin lineaarialgebrakirjastoihin ilman Fortranin osaamista. Moler kaupallisti ohjelmansa ja perusti yhteistyökumppaniensa kanssa MathWorks-yrityksen vuonna 1984. MathWorks kehittää Matlabia edelleen aktiivisesti; ohjelman uusin versio R2017b, järjestyksessä 50., julkaistiin vuoden 2017 jälkipuoliskolla. Matlabin ohjelmointikielen erikoisuus on se, että se on matriisipohjainen eli reaalimuuttujatkin ovat aidosti 1x1-matriiseja. Monet numeeriset ongelmat ovat luontevasti tulkittavissa matriisien kielelle, joten Matlab, pohjautuen optimoituun LINPACK-lineaarialgebrakirjastoon, on erittäin tehokas ja suosittu numeerisen laskennan ohjelmisto. Sen matriisioperaatiot ovat huomattavasti nopeampia kuin esimerkiksi Mathematican. Toisin kuin Mathematicassa, Matlabissa ei ole symbolisia toimintoja, vaikkakin ne on mahdollista saada käyttöön lisäpakettien avulla. Matlab on kaupallinen ohjelmisto ja sen lisenssi on varsin kallis. Opiskelijoiden on tosin mahdollista saada lisenssi paljon halvemmalla. Matlabin alkeiden opetteluun voi käyttää myös vapaita ohjelmistoja kuten Scilabia tai GNU:n Octavea. Nämä ohjelmat eivät ole täydellisiä Matlab-klooneja, mutta perusohjelmointi niillä on hyvän samankaltaista. Matlab on saatavissa kaikille järkeville alustoille: Unix, Mac OS X ja Windows.

Linkkejä muihin materiaaleihin

Käyttöliittymä ja tulkki

Matlab avattaessa avautuu seuraavankaltainen näkymä:

Tärkein on kuvassa keskellä oleva komentotulkki (command window), jossa voi interaktiivisesti ajaa Matlab-komentoja. Vasemmalla on esillä tämänhetkinen työskentelyhakemisto (current folder). Oikealla on tämänhetkinen tila (workspace) sekä komentohistoria (command history).

Tulkin kehote >> erottuu tulosteriveistä, joilla kehotetta ei ole. Komento suoritetaan painamalla ENTER-näppäintä. Esimerkiksi:

x = 5 + 1
x =

     6

phi = (1+sqrt(5))/2
phi =

    1.6180

Tästä huomataan, että Matlab tekee jokaisen komennon/rivin jälkeen tulosteen. Tämä interaktiivisuus ei aina ole tarpeen, joten tulostuksen voi vaimentaa päättämällä rivin puolipisteellä ;. Puolipisteen tai pilkun avulla voi myös kirjoittaa useamman komennon samalle riville, mikäli näin haluaa tehdä. Prosenttimerkki % toimii kommentin aloitusmerkkinä. Matlabin syntaksi on suoraviivaista ja tuttua kaikille C-tyyppisiä kieliä ohjelmoinneille. Syntaksiin ei tässä dokumentissa juuri puututa tämän enempää – se tullee selväksi esimerkkien ja lukijan omien kokeilujen kautta. Sanottakoon tässä kuitenkin, että Matlab on dynaamisesti tyypitetty kieli. Huomaa, että luomalla uuden muuttujan voi ylikirjoittaa Matlabin komentoja: esimerkiksi sin-nimisen muuttujan määrittely estää sin-funktion käytön; tätä tulee välttää.

Matlab-komennoista saa lisätietoa kirjoittamalla help komento tai kirjoittamalla komennon hakukenttään (ks. yo kuvassa oikealla ylhäällä). Kirjoittamalla pelkästään help saa näkyviin aihealueittain järjestetyn listan komennoista. Komento lookfor sana listaa kaikki komennot, joiden ohjeen ensimmäiseltä riviltä löytyy annettu sana (1. rivi antaa lyhyen kuvauksen koko komennosta. Lisäksi lookfor sana -all etsii haluttua sanaa koko ohjeesta. Myös Mathworksin sivuilta löytyy paljon tarpeellista tietoa.

Harjoitus. Tutustu Matlabin käyttöliittymään ja kokeile peruslaskutoimituksia. Etsi help-komennon avulla tietoa käskystä format ja tulosta kultainen leikkaus suuremmalla tarkkuudella. Onko tulostustarkkuudella merkitystä laskutarkkuuden kannalta? Selvitä itsellesi mitä tekee format compact.

Kokeilemalla esimerkiksi edeltävää esimerkkiä huomataan, että määritellyt muuttujat x ja phi jäävät muistiin ja ne näkyvät tämänhetkisessä tilassa (workspace). Komento who listaa kaikki määritellyt muuttujat, clear poistaa kaikki määritellyt muuttujat ja clear x poistaa muuttujan x. Kokeile myös komentoa whos. Matlabissa on myös laskimen tapaan muuttuja ans, joka muistaa viimeisimmän asetetun arvon tai tuloksen. Komento clc tyhjentää kehotteen.

Komennolla disp voi tulostaa. Esimerkiksi

disp(phi)
    1.6180

Myös merkkijonojen tulostus onnistuu (Matlabissa merkkijono aloitetaan ja päätetään heittomerkillä).

disp('Hello world!')
Hello world!

Merkkijonot ovat Matlabissa merkkitaulukoita, joten monimutkaisempien tulosteiden rakentaminen käy seuraavaan tapaan (ks. myöhemmin vektorien ja matriisien luonti).

disp(['Muuttujan phi arvo on ' num2str(phi) '.'])
Muuttujan phi arvo on 1.618.

Matlabissa on valmiiksi ohjelmoituna tavanomaiset alkeisfunktion, kuten abs, sin, cos, tan, exp, log, sinh, jne. Näistä löytyy lisätietoa kirjoittamalla help elfun. Erityisfunktiot, kuten Besselin funktiot, löytyvät komennolla help specfun.

Harjoitus. Vieraile Matlab-dokumentaation sivulla Complex Numbers ja selvitä miten kompleksiluvut toimivat. Kokeile alkeisfunktioden kutsumista kompleksilukuparametrilla.

Matriisit

Kuten edellä jo mainittiin, Matlab on ensisijaisesti matriisikieli. Yksinkertaisimmillaan matriisi luodaan antamalla riveittäin sen alkiot:

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

     1     2     3
     4     5     6
     7     8     9

Määrittelyn voi myös jakaa useammalle riville, jolloin puolipistettä ei tarvita:

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

     1     2     3
     4     5     6

Kolmella pisteellä voi pitkän rivin katkaista kahdelle eri riville:

A = [1 2 3 4 5 6 ...
7 8 9
10 11 12 13 14 15 16 17 18]
A =

     1     2     3     4     5     6     7     8     9
    10    11    12    13    14    15    16    17    18

Matriisin alkioiden määrittelyt voi erottaa myös pilkulla välilyönnin sijaan. Luonnollisesti alkiot voi määritellä laskutoimituksen avulla.

A = [1, 2*sin(pi/4), exp(2)]
A =

    1.0000    1.4142    7.3891

Vektorin pituuden saa selville komennolla length:

length(A)
ans =

     3

Matriisin dimensiot taas selvittää komento size, kokeile! Matlabissa on paljon valmiita komentoja matriisien luontiin. Esimerkiksi eye antaa identiteettimatriisin:

eye(3)
ans =

     1     0     0
     0     1     0
     0     0     1

Harjoitus. Tutustu help-komennon avulla komentoihin rand, ones, zeros.

Matriiseille on määritelty tutut laskuoperaation +, -, *, /, \, ^, ' eli yhteenlasku, vähennyslasku, kertolasku, oikeata ja vasemmalta jakaminen, potenssiinkorotus ja transponointi. Ohjelmoijan tulee luonnollisesti vaatia, että matriisit ovat kooltaan yhteensopivat kunkin laskutoimituksen kohdalla. Laskuoperaatioista pääsee helpoimmin jyvälle kokeilemalla; jakolaskuun palataan myöhemmin lineaarisen yhtälöryhmän ratkaisemisen yhteydessä. Mainitsemisen arvoinen lisätieto tässä on skalaarin vähentäminen (tai lisääminen):

A = eye(3)
A - 1
A =

     1     0     0
     0     1     0
     0     0     1


ans =

     0    -1    -1
    -1     0    -1
    -1    -1     0

Matlab osaa siis automaattisesti tulkita tässä vähennettävän skalaarin matriisilla ones(3) kerrottuna. Lisäksi on syytä huomauttaa, että ' on hermiittinen transponointi eli varsinaisen transponoinnin lisäksi suoritetaan kompleksikonjugointi. Tavallisen transponoinnin saa käskyllä transpose tai notaatiolla .'.

Kertolaskusta, potenssiinkorotuksesta ja jakolaskusta on myös pisteittäiset versiot .*, .^, ./. Siis vaikkapa .* laskee kahden matriisin tulon pisteittäin (alkoittain) eikä laske matriisituloa. Esimerkiksi

A = [1 2; 3 4], B = [1 1; 2 2]
A .* B
A =

     1     2
     3     4


B =

     1     1
     2     2


ans =

     1     2
     6     8

Monet Matlab-funktiot, kuten yllä mainitut alkeisfunktiot, toimivat pisteittäin. Esimerkiksi funktio sin laskee matriisin jokaisen alkion sinin:

sin(eye(2))
ans =

    0.8415         0
         0    0.8415

Kaksoispistenotaatio ja indeksointi

Usein on tarvetta tuottaa vektoreita, joiden alkioiden etäisyydet ovat samat. Tähän tulee avuksi seuraava kaksoispistenotaatio:

  start:step:end

joka tuottaa vektorin, jossa on muotoa luvut start, start + step, start + 2* step, ..., jotka ovat pienempiä tai yhtäsuuria kuin end. Esimerkiksi:

1:0.1:2
1:0.3:2
ans =

  Columns 1 through 7

    1.0000    1.1000    1.2000    1.3000    1.4000    1.5000    1.6000

  Columns 8 through 11

    1.7000    1.8000    1.9000    2.0000


ans =

    1.0000    1.3000    1.6000    1.9000

Kohdan step voi jättää välistä, jolloin sen arvoksi oletetaan 1. Tässä koodi, joka tuottaa saman tuloksen kuin yllä käyttäen for-silmukkaa. Tätä tapaa ei tietenkään tule käyttää, mutta tässä onkin tarkoitus esitellä for-silmukan syntaksi ja miten vektoriin lisätään alkioita.

X = [];
s = 1; e = 2; t = 0.1;
for n = 0:(e-s)/t
    X = [X s + n*t];
end
X
X =

  Columns 1 through 7

    1.0000    1.1000    1.2000    1.3000    1.4000    1.5000    1.6000

  Columns 8 through 11

    1.7000    1.8000    1.9000    2.0000

Harjoitus. Selvitä miten komento linspace toimii ja miten se eroaa kaksoispistenotaatiosta. Komento logspace saattaa myös olla myöhemmin hyödyksi.

Seuraavassa koodissa esimerkinomaisesti matriisin alkioiden, rivien ja sarakkeiden valinnasta.

A = [1 2 3 4; 5 6 7 8; 9 10 11 12]
A(2,3)                                   % Kohdassa (2,3) oleva alkio.
A(8)                                     % Kohdassa (2,3) oleva alkio.
A(:,1)                                   % 1. sarake
A(1,:)                                   % 1. rivi
A(:)                                     % Koko matriisi pystyvektorina.
A(2:3,2:3)                               % Alimatriisi, rivit 2-3, sarakkeet 2-3
A(:,[1,4])                               % Kaikki rivit, sarakkeet 1 ja 4.
B = logical([0 1 0 0; 1 0 1 0; 0 0 1 0]) % Luodaan "logiikkamatriisi", joka kertoo mitkä A:n alkiot valitaan.
A(B)                                     % Valitsee A:n ne alkiot, joita vastaavassa kohdassa B:ssä oli ykkönen, järjestys ylhäältä alas, vasemmalta oikealle.
A =

     1     2     3     4
     5     6     7     8
     9    10    11    12


ans =

     7


ans =

     7


ans =

     1
     5
     9


ans =

     1     2     3     4


ans =

     1
     5
     9
     2
     6
    10
     3
     7
    11
     4
     8
    12


ans =

     6     7
    10    11


ans =

     1     4
     5     8
     9    12


B =

     0     1     0     0
     1     0     1     0
     0     0     1     0


ans =

     5
     2
     7
    11

Matriisia voi myös muokata analogisesti:

A(2,2:3) = [0 0]    % Muutetaan 2. rivin sarakkeita 2 ja 3.
A(:,4) = [12; 8; 4] % Muutetaan 4. sarake.
A =

     1     2     3     4
     5     0     0     8
     9    10    11    12


A =

     1     2     3    12
     5     0     0     8
     9    10    11     4

Harjoitus. Määrittele jokin matriisi ja lisää sen 1. vaakarivi vakiolla 3 kerrottuna sen 2. vaakariviin. Kaksoispistenotaatiolla saa siis helposti toteutettua matriisien elementaariset vaaka- ja pystyrivioperaatiot. Näin on mahdollista ohjelmoida matriisin muokkaaminen vaikkapa redusoituun porrasmuotoon.

Hieman merkkijonoista

Matlabissa merkkijonoista voi hakea ja niitä voi muokata säännöllisillä lausekkeilla. Seuraava esimerkki etsii pidemmästä merkkijonosta merkkijonon aa alkamiskohdat (ei etsi päällekkäisiä esiintymiä).

regexp('abbabaabaaa', 'aa')
ans =

     6     9

Luonnollisesti aitojen säännöllisten lausekkeiden kirjoittaminen onnistuu myös. Seuraavassa etsitään pitkästä merkkijonosta kaikki osasanat, jotka alkavat kirjaimella k ja päättyvät kirjaimeen o.

s = 'koko kokoa koko kokko kokoon';
regexp(s, 'k[\w]*o')                 % palauta täsmäävien sanojen alkamiskohdat
regexp(s, 'k[\w]*o', 'match')        % palauta täsmäävät sanat
ans =

     1     6    12    17    23


ans = 

    'koko'    'koko'    'koko'    'kokko'    'kokoo'

Osasanat korvataan seuraavasti:

regexprep(s, 'k', 'o')    % vaihdetaan k-kirjaimet o-kirjaimiksi
ans =

oooo ooooa oooo ooooo ooooon

Lue lisää säännöllisistä lausekkeista sivulta Match regular expression (case sensitive).

m-tiedostot ja omien funktioiden kirjoittaminen

Skriptit

Suoritettavat komennot on interaktiivisen kokeilun jälkeen syytä tallettaa muistiin omaksi tiedostokseen. Matlabissa koodi on tapana tallettaa .m-päätteiseen tekstitiedostoon. Koodin voi kirjoittaa millä tahansa editorilla, mutta kätevintä se lienee tehdä Matlabin omalla editorilla. Matlabin editorin saa auki valitsemalla käyttöliittymästä New -> Script. Saman saa aikaan kirjoittamalla kehotteeseen edit tiedostonimi. Kopioidaan seuraava koodi editoriin ja talletetaan se tiedostoksi testi.m nykyiseen työskentelyhakemistoon (tämä on tärkeää, jotta Matlab löytää tiedoston).

x = randn(1, 100);
m = mean(x)
s = std(x)

Tämän jälkeen kirjoitetaan tulkissa testi (voi myös klikata editorissa Run).

testi
m =

    0.1883


s =

    1.0040

Lopputulos on siis sama kuin mikäli tiedoston testi.m komennot olisi kirjoitettu tulkkiin. Huomaa, että tulostus tapahtui aivan samoin kuin tulkissa, minkä lisäksi määritellyt muuttujat x, m ja s jäivät voimaan tiedoston komentojen suorittamisen päätyttyä. Tällä tavoin voi siis vaikkapa alustaa muuttujille sopivat arvot interaktiivisen istunnon aluksi. Huomaa, että m-tiedostossa voi samaan tapaan kutsua toisen m-tiedoston koodia. Muista, että kommenttimerkin % avulla voi (ja näin kannattaakin tehdä) kommentoida m-tiedoston sisältämää koodia. On hyödyllistä huomata, että editorissa maalatun koodin voi suorittaa tulkissa painamalla F9-näppäintä.

Tähän tapaan luotua m-tiedostoa kutsutaan skriptiksi, eikä se voi ottaa kutsuttaessa parametreja. Seuraavassa aliluvussa katsotaan miten m-tiedostoihin voi tallettaa funktioita, joita kutsuttaessa voi luonnollisesti antaa parametreja.

Matlabilla on myös mahdollista tehdä dokumentteja, jotka yhdistelevät koodia, tekstiä ja muuta muotoilua ja jonka voi kääntää html- tai pdf-tiedostoksi. Tämäkin ohje on vain yksi m-tiedosto, jonka publish-komento kääntää siistiin muotoon. Tästä lisätietoja sivulla Publishing Markup.

Funktiotiedostot

Funktiomäärittelyt kirjoitetaan samoin m-tiedostoihin, ja niiden syntaksi on seuraavanlainen:

function y = avg(x)
y = sum(x)/length(x);

Jotta Matlab tunnistaa m-tiedoston funktiomäärittelyksi eikä skriptiksi, tulee avainsanan function löytyä m-tiedoston ensimmäiseltä riviltä. Funktion sisältämälle tiedostolle tulee antaa sama nimi kuin määritettävälle funktiolle (tässä avg.m), jotta Matlab löytää halutun funktion. Huomaa, että funktion komennot on syytä päättää puolipisteeseen, muutoin syntyy paljon ylimääräistä tulostetta. Muista ohjelmointikielistä tutusti funktion sisältämät muuttujamäärittelyt ovat lokaaleja eivätkä jää voimaan funktion suorituksen päätyttyä.

Yllä y on funktion palauttama arvo, joka asetetaan tämän funktion viimeisellä rivillä. Funktion nimi on avg ja sen ottaman parametrin nimi on x. Funktio voi luonnollisesti ottaa useampia parametreja: avg(x,y) jne.

Funktio voi myös palattaa useamman arvon. Tällöin syntaksi on seuraava:

function [a,b] = foobar(x)
a = sum(x);
b = sum(x)/length(x);

Tätä funktiota kutsutaan seuraavasti:

[a,b] = foobar(x)

Mikäli ei sattumalta välitä kaikista paluuarvoista, voi niistä napata vain ensimmäisen:

a = foobar(x)

Edeltävällä hakasulkunotaatiolla voi valita yhden tai useamman palautusarvon useamman kuin kahden palautusarvon tapauksessa. Mikäli funktio ei palauta mitään, voi palautusarvon jättää kirjoittamatta:

function func(x)
...

Matlab tukee myös optionaalisia parametreja eli funktion kaikkia kutsuparametreja ei tarvitse antaa. Tällöin funktion ohjelmoinnista riippuu tuottaako tämä virheen vai muuttaako se funktion käyttäytymistä jotenkin (ks. jäljempänä yksityiskohtaisemmin omien funktioiden ohjelmoinnista).

Yhdessä m-tiedostossa voi määritellä useamman funktion, mutta nämä ylimääräiset funktiot ovat lokaaleja – tiedoston varsinainen "pääfunktio" voi kutsua näitä muita funktioita. Tässä tapauksessa kaikki funktiomäärittelyt tulee päättää avainsanalla end, muutoinhan olisi epäselvää milloin mikin funktio päättyy.

Anonyymit funktiot

Useat Matlab-funktion ottavat funktioita parametriksi, kuten funktion nollakohdan etsivä fzero-komento. Oletetaan, että tiedostossa fun.m on seuraava sisältö:

function y = fun(x)
y = 1./(exp(x) + sin(x)) - x.^2;

Tiedoston määrittelemä funktio laskee siis funktion $\frac{1}{e^x + \sin{x}} - x^2$ arvon. Etsitään funktion nollakohtaa välillä $[-0.5,1]$:

fzero(@fun, [-0.5, 1])
ans =

    0.6350

Tässä @-merkin rooli on siinä, että se luo viittauksen (funktiokahvan) funktioon fun. Viittauksia voi tallettaa muuttujiin, eli seuraavakin koodi toimii odotetusti:

f = @fun;
fzero(f, [-0.5, 1])
ans =

    0.6350

On kuitenkin vaivalloista kirjoittaa m-tiedosto jokaista pikku reaalifunktiota kohden, joten Matlab mahdollistaa ns. anonyymin funktion luonnin. Seuraava koodi tekee edeltävän asian ilman, että erillistä m-tiedostoa tarvitsee luoda:

f = @(x) 1./(exp(x) + sin(x)) - x.^2;
fzero(f, [-0.5, 1])
ans =

    0.6350

Anonyymin funktion voi määritellä myös parametrista riippuvaksi:

a = pi;
f = @(x) 1./(exp(x) + sin(x)) - a*x.^2;
fzero(f, [-0.5, 1])
ans =

    0.4090

Lue lisää anonyymeistä funktioista täältä.

Kutsuttavien m-tiedostojen sijainti

Oletuksena Matlab etsii kutsuttavia skriptejä tai funktioita komennon path listaamista hakemistoista. Näistä ensimmäisenä on nykyinen työskentelyhakemisto. Hakemistoja voi lisätä listaan addpath-komennolla. Hakemistoilla on merkitystä siinä vaiheessa, kun niissä on määritelty useampi samanniminen m-tiedosto, jolloin Matlab valitsee ensimmäisen hakemistolistassa järjestyksessä vastaantulevan tiedoston. Komennolla which saa aina selville mitä tiedostoa Matlab milloinkin pyrkii käyttämään. Huomaa, että m-tiedostolle ei kannata antaa samaa nimeä kuin jollekin Matlabin sisäiselle (builtin) funktiolle, sillä sisäisiä funktioita kutsutaan aina ennen kuin m-tiedostoja.

Matlabin työkalupakeista (Toolbox)

Matlabissa on käytettävänä erilaisia työkalupakkeja (toolbox). Nämä ovat typpillisesti kokoelma funktioita eli m-tiedostoja, mutta ne voivat sisältää myös sovelluksia ja dataa.

Istunnon aikana sinä hetkenä käytettävissä olevat (Mathworksin) työkalupakit löytää komennolla ver. Komennon suorittamalla nähdään, että käytössä on esimerkiksi "optimization toolbox". Näin ollen kaikki tämän pakin funktiot ovat käytettävissä. Mathworksin sivuilta https://mathworks.com/products.html löytyy lisää paketteja, mutta ne ovat maksullisia. Mathworksin itse tekemien pakkien lisäksi löytyy kolmansien osapuolien (matlabin käyttäjien) kokoamia pakkeja. Näitä löytyy esimerkiksi Matlab Centralista. Näiden käyttöönotto selviää esimerkiksi pakettien readme-tiedostoista. Jos pakki on vain kansio m-tiedostoja, addpath komennon käyttö riittää. Esimerkiksi kansoissa G:\ELBM on nippu m-tiedostoja. Nämä sadaan käyttöön seuraavasti.

help ELBM          % Matlab ei tunne tätä funktiota.
addpath('G:\ELBM')
help ELBM          % Nyt löytyy.
path

Yksityiskohtaisemmin omien funktioiden ohjelmoinnista

Matlab-ohjelmoinnissa on luonnollisesti käytössä jokaisesta ohjelmointikielestä tutut kontrollirakenteet if, for, while jne. Näitä ei tässä materiaalissa käydä läpi tarkemmin; lisätietoa asiasta löytyy sivulta Control Flow. Kontrollirakenteiden yhteydessä joutuu usein käyttämään loogisia operaatioita, ks. help ops. Monia Matlabin komentojen toteutuksia voi tutkia type-komennolla. Kokeile vaikkapa komentoa type('quad').

%
% Seuraavaksi esittelemme funktion, joka laskee piin likiarvon sarjan
% $4(\frac11-\frac13+\frac15-\frac17+\cdots)$ avulla (suppeneminen on
% erittäin hidasta). Tämän funktion käyttäytyminen riippuu sille annettujen
% parametrien lukumäärästä; se myös antaa virheilmoituksen mikäli annettu
% parametri ei kelpaa.
%
% <include>pi_approx.m</include>
%
% Kutsujan pyytämien paluuarvojen lukumäärän palauttaa |nargout|-funktio.
% Tämä on joskus hyödyllinen tieto, kun jonkin palautusarvon laskeminen on
% työlästä; mikäli sitä ei pyydetä, niin sitä ei tarvitse laskeakaan.
%
% Syötteen pyytäminen käyttäjältä tapahtuu |input|-funktion kautta.
% Funktiolle annetaan parametrina teksti, joka tulostetaan näytölle.
% Funktio palauttaa arvon, jonka käyttäjä antaa.
% Alla on edellinen funktio muokattuna, jossa käyttäjältä kysytään erikseen
% iteraatiomäärää.
%
% <include>pi_approx_inter.m</include>
%
% Dynaamisen tyypityksen henkeen funktiot olisi hyvä ohjelmoida niin, että
% ne kelpuuttavat erityyppisiä parametreja. Esimerkiksi reaalifunktion
% arvon annetussa pisteessä laskeva funktio voi aivan hyvin kelpuuttaa
% parametrikseen vektorin tai matriisin – arvo lasketaan vain pisteittäin
% matriisin alkioille. Tämä on Matlabissa jokseenkin vaivatonta
% pisteittäisten laskutoimitusten avulla, kuten seuraava esimerkki näyttää.
% Tämän esimerkin funktio laskee sinin arvon parametrikseen saamilleen
% luvuille Taylorin sarjan avulla.
%
% <include>sin_taylor.m</include>
%

Oletusiteraatiomäärällä saadaan neljän desimaalin tarkkuudella samat tulokset kuin Matlabin sin-funktiolla.

sin_taylor([0 0.1; 0.5 1]), sin([0 0.1; 0.5 1])
ans =

         0    0.0998
    0.4794    0.8415


ans =

         0    0.0998
    0.4794    0.8415

Harjoitus. Ohjelmoi funktio, joka ottaa syötteeksi neliömatriisin A (tarkasta ettei käyttäjä anna väärän kokoista matriisia) ja laskee matriisin $e^A$ käyttämällä sarjakehitelmästä $e^A = A + 1/2!A^2 + 1/3!A^3 + \ldots$ niin montaa termiä kuin käyttäjä pyytää.

Debuggauksesta

Jos koodi päättyy virheeseen, Matlabin tulosteesta näkee miksi ja millä rivillä tämä tapahtui. Tälle riville pääsee Matlabin tulostamasta linkistä. Nämä voi myös selvittää lasterror-funktiolla.

testi = cell(1,2);  testi{1} = 'a'; testi{2} = 2;
cell2mat(testi);
err = lasterror
err.message
err.stack
err.identifier

Koodia voi debugata ja ajaa rivi kerrallaan käyttämällä Matlabin debuggeria. Debuggerin käynnistys onnistuu laittamalla koodiin pysähdyksen. Tämä onnistuu avaamalla debugattava m-tiedosto ja klikkaamalla rivinumeron vieressä olevaa viivaa. Tällöin viivan korvaa punainen pallo. Toinen tapa on käyttää dbstop komentoa. Testataan tätä Matti Vuorisen matematiikan menetelmäkurssin luentomonisteen sivulta 64 löytyvällä m-koodilla.

dbstop in e310.m at 15
e310

Koodi asetti pysähdyksen riville 15. Kun tiedoston ajoi, Matlab meni debugger-moodiin, jonka tunnistaa kehotteen edessä löytyvästä K-kirjaimesta. Matlabin komentoja voi käyttää normaalisti debugger-moodissa. Koodia voi nyt suorittaa rivi kerrallaan. Seuraava rivi suoritetaan dbstep-komennolla. Komento dblcont suorittaa koodia seuraavaan pysähdykseen asti. Debuggaus-moodista pääsee pois dbquit-komennolla. Tämä ei suorita koodia loppuun vaan keskeyttää koodin ajamisen.

who
n
dbstep
who
dbstep
who
dbcont
n
dbquit
n

Debuggerin avulla voidaan siis tarkistaa, että muuttujilla on koodin edetessä ne arvot kuin pitääkin.

Joskus oikein toimiva koodi on turhan hidas. Koodin eri osien kestoa voi mitata tic- toc-komentoparin kanssa.

tic
for i=1:10000
  v(i)=i;
end;
toc
tic
v=1:10000;
toc
Elapsed time is 0.001650 seconds.
Elapsed time is 0.000280 seconds.

Lisää aikaan liittyviä funktioita löytyy komennolla help timefun.

Datan tuominen Matlabiin

Datan tuomiseen on monta funktiota kuten load, fscanf, textscan, importdata ja xlsread. Lisää funktioita sekä näiden dokumentointia löytyy sivulta Supported File Formats for Import and Export.

Tutkitaan aluksi sisäistä funktiota load. Funktio pystyy lataamaan mat-tiedoston tai ASCII-muotoisen tiedoston, joka sisältää ainoastaan numeerisia arvoja. Mat on matlabin oma tiedostotyyppi, johon voi tallentaa muuttujia työtilasta save-komennolla.

m = 1;
M = [1,2; 3,4];
save('test.mat', 'M');  % Tallennetaan matriisi M.
who;                    % Tarkistetaan mitä työtilasta löytyy.
clear;
who;                    % Nyt työtilassa ei ole mitään.
load('test.mat');
who;                    % M ladattu.

Jos save-komennolle ei anna erikseen tallennettavia muuttujia, se tallentaa koko työtilan. Jos komennolle ei anna tiedoston nimeä, muuttujat tallennetaan tiedostoon matlab.mat.

Tiedostossa plotdata.txt on kaksi saraketta numeerisia arvoja. Datan saa ladattua matriisiksi komennolla

M = load('plotdata.txt');

Joskus data sisältää muutakin kuin numeerisia arvoja, jolloin voi turvautua muihin funktioihin. Tarkastellaan seuraavan xslx-muotoisen datan tuomista xlsread-funktiolla.

xlsx-tiedosto

Data sisältää kolme lehteä (sheet), joissa ensimmäisessä on yksinkertainen taulukko. Tarkistetaan komennolla help xlsread, kuinka tämä voidaan lukea.

Ladataan taulukon tiedot seuraavalla komennolla

M = xlsread('xlsxtest.xlsx')

Havaitaan, ettei datan merkkijonosaraketta eikä merkkijonoriviä luettu ja data tuli tallennetuksi muuttujaan M. Merkkijonot saadaan tallennetuksi muuttujaan R seuraavasti:

[M,R,U] = xlsread('xlsxtest.xlsx')

Muuttuja U sisältää "raakadatan". Sen tyyppi on solutaulukko(cell array) ja jokaisen solun sisällä on yksi luettu alkio. Yleisesti solu voi sisältää minkä tahansa datatyypin. Solusta saa ulos siihen tallennetun tiedon viittaamalla siihen aaltosulkeilla tai käyttämällä funktiota cell2mat.

class(U)
class(U(1,2))
class(U{1,2})
N=cell2mat(U)
N=cell2mat(U(2:4,2:5))
sum(sum(abs(N-M)))

Siirrytään seuraavaan lehteen. Tässä datassa on yksi NA-arvo ja kaksi puuttuvaa arvoa. Tämän lisäksi on uusi merkkijonosarake. Katsotaan kuinka xlsread selviää tästä.

[M,R] = xlsread('xlsxtest.xlsx',2)

Havaitaan, että edelliset puuttuvat arvot oli korvattu (double) arvolla NaN ja merkkijonorivi oli luettu oikein muuttujaan R. Kolmannessa lehdessä on desimaalipilkkuja. Nämä ovat epätriviaaleja siinä mielessä, että Matlabissa desimaalit ovat pisteitä.

[M,R] = xlsread('xlsxtest.xlsx',3)

Havaitaan, että xlsread osasi muttaa desimaalipilkut pisteiksi.

Tehtävä: kokeile lukea seuraavat datat käyttäen importdata-funktiota. (spoiler:ongelmia luvassa)

pistedata erikoisdata pilkkudata

Matlabin import wizard (File->Import Data...) käyttää funktiota importdata. Kuten tehtävästä voi todeta se ei aina ole mutkatonta ja joskus on tarpeen käyttää muita funktioita.

Yksi helppo tapa tuoda tehtävän datat on lukea ne excelissä, tehdä xlsx-tiedosto ja käyttää xlsread-funktiota. Toinen tapa on käyttää textscan tai fscanf-funktiota. Tällöin data luetaan rivi kerrallaan ja tehdään riveille tarvittavia muutoksia. Esimerkiksi tiedoston dotdata voi lukea seuraavasti:

fo=fopen('G:\matlab\dotdata.txt'); % Avataan tiedosto.
otsake=fscanf(fo,'%s %s %s %s',4); % otsikkorivi erilliseen muuttujaan.
m=textscan(fo,'%s %f %f %f %f');   % loput riveistä on samaa muotoa.
fclose(fo);                        %Tiedoston sulkeminen.
M=[m{2},m{3},m{4},m{5}];           % Numeerisen matriisin muodostaminen

Tiedoston commadata.txt voi myös ensin muokata pistedataksi fscanf-funktiolla.

fo=fopen('G:\matlab\commadata.txt');
pilkkudata=fscanf(fo,'%c')           % '%c' säilyttää välilyönnit
fclose(fo);                          % Tiedoston sulkeminen.
pistedata=regexprep(o,',','.')       % Korvataan merkkijonossa pilkut pisteiksi.
fc=fopen('G:\matlab\commadatamod.txt','w'); % Avataan ja tarvittaessa luodaan tiedosto.
fprintf(fc,'%c',pistedata);
fclose(fc);

Lisää apua datan lukemisen ongelmiin löytyy sivulta https://mathworks.com/help/matlab/import_export/recommended-methods-for-importing-data.html.

Datan vienti Matlabista

Edellä käsitellyn save-komennon lisäksi datan vientiin on useita read-funktioita vastaavat write-funktiot. Näistä löytyy tietoa edellisistä kahdesta linkistä. Tallennetaan esimerkiksi matriisi

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

excel-tiedostoksi. Oikean syntaksin löytää komennolla help xlswrite.

xlswrite('xlswritetest.xlsx', M);

Tehtävä. Tallenna yo. matriisi M tekstitiedostoon, siten että sarakkeilla on otsikot "yksi, "kaksi" ja "kolme" käyttäen funktiota fprinf.

2D-grafiikka

Plot

Komento plot piirtää annetut pisteet koordinaatistoon ja yhdistää pisteet viivalla (oletusarvoisesti). Esimerkiksi

x = [0 1 2];
y = [1 -1 0.5];
plot(x, y);

Matlab asettaa x- ja y-akselien skaalat oletusarvoisesti piirrettävän datan mukaan. Piirrettävää aluetta voi säätää seuraavasti:

axis([-1,1,-2,2]) % x-akseli -1 - 1, y-akseli -2 - 2

Akselien skaalaa säädetään esimerkiksi komennoilla axis equal (sama jakoväli kaikilla akseleilla) ja axis square (akselien näytettävät osuudet yhtä pitkät graafisesti). Perustilaan pääsee kirjoittamalla axis normal. Huomaa, että nämä komennot tulee antaa plot-komennon jälkeen (Matlab ei voi tuntea dataa ennen piirtämistä). Komento grid näyttää kuvassa hilaviivat. Kokeile!

Komennot axis ja grid ja muutkin muokkaavat implisiittisesti tämänhetkistä kuvaa (figure). Matlabissa voi olla samanaikaisesti auki useampi kuvaikkuna; uuden kuvaikkunan saa auki komennolla figure. Komento luo uuden kuvaikkunan, johon viitataan juoksevalla järjestysnumerolla. Komento figure(n) vaihtaa nykyistä kuvaikkunaa, mikä tarkoittaa, että komennot axis, grid jne. kohdentavat toimintansa uuteen ikkunaan. Seuraava esimerkki piirtää sini- ja kosinifunktion eri ikkunoihin eri skaaloilla.

clear;
x = linspace(0, 2*pi);
y = sin(x); z = cos(x);
figure;                  % Luo uusi kuvaikkuna, nyt käsitellään implisiittisesti kuvaa nro 1.
plot(x, y);
axis equal;
figure;                  % Nyt käsitellään implisiittisesti kuvaa nro 2.
plot(x, z, 'x');         % Piirretään vain datapisteet, ei yhdistäviä janoja.
axis([0, 2*pi, -2, 2]);
figure(1);               % Nyt käsitellään jälleen kuvaa nro 1.
grid;
figure(2);

Kuvaikkunan sisällön voi tallettaa tiedostoon yleisimmissä formaateissa komennolla print. Esimerkiksi print('kuva', '-dpdf') tallettaa tämänhetkisen kuvan PDF-muodossa työskentelyhakemiston tiedostoon kuva.pdf. Parametri -deps tallettaa kuvan EPS-muodossa. Tallennus onnistuu myös kuvaikkunan valikosta: File -> Save As.

Samaan kuvaan on luonnollisesti mahdollista piirtää useampi kuvaaja samaan kuvaan. Komennon plot kutsuminen uudelleen pyyhkii edellisen kuvan. Tämän estämiseksi tulee käyttää käskyä hold on. Tämän jälkeen kaikki plot-komennot piirtävät vanhan kuvan päälle kunnes annetaan käsky hold off. Plot-komennolle voi antaa myös useamman datan samanaikaisesti seuraavalla tavalla:

  plot(x1, y1, 'mjono1', x2, y2, 'mjono2', x3, y3, 'mjono3, ...)

Annetuilla merkkijonoilla (LineSpec) on mahdollista muokata käyrän väriä ja muotoa. Kokeile esimerkiksi arvoja '--or' ja ':^g'. Lisätietoa sivulta Plot.

Harjoitus. Anna plot-käskylle parametriksi vektori kompleksilukuja. Selvitä edeltävän linkin kautta miten piirtää datapisteet ilman yhdistäviä viivoja.

Fplot

Edellä esitelty plot-komento piirsi kahden vektorin määräämät tasopisteet ja yhdisti ne janoilla. Tämä komentohan piirtää vain annetun datan, jonka ei edes tarvitse esittää funktiota. Mikäli haluaa piirtää jonkin Matlabin määrittelemän funktion tai itse ohjelmoidun funktion, ei ole tarpeen toimia kuin edellä, vaan voi käyttää komentoa fplot:

fplot(@(x) x.^2 .* sin(x), [0, 4*pi]);

Matlab hoitaa nyt annetun piirtovälin diskretisoinnin automaattisesti.

Parametrisen funktion piirto tasoon onnistuu myös komennolla fplot:

fplot(@(t) cos(3*t), @(t) sin(2*t));

Napakoordinaatistossa piirtäminen onnistuu luonnollisesti myös:

t = linspace(0,2*pi);
polar(t, sin(2*t).^2);

Monta kuvaajaa samaan kuvaan

Komento subplot jakaa kuvan (figure) useaan alikuvaan (subplot). Komento subplot(m, n, p) jakaa kuvan m*n alikuvaan, jotka on järjestetty m*n-ruudukoksi. Parametri p valitsee nykyiseksi alikuvaksi järjestyksessä p:nnen alikuvan (lasketaan riveittäin). Seuraavassa esimerkissä luodaan kaksi vaakatasossa vierekkäistä alikuvaa, piirretään niihin kuvaajat ja asetetaan kuville vielä otsikot (tutki tarkemmin käskyä title).

close all;                    % sulkee kaikki kuvaikkunat
subplot(1, 2, 1);
fplot(@(x) sin(x), [0, 2*pi]);
title('Sinikäppyrää');
subplot(1, 2, 2);
fplot(@erf, [-1,1]);
title('Jotain aivan muuta');

Muuta hyödyllistä tietoa

Tässä vielä muutama esimerkki tekstin lisäyksestä kuviin, viivapaksuuden muuttamisesta jne. Kokeile!

figure()
x = linspace(-2*pi, 2*pi); y1 = sin(x); y2 = cos(x);
plot(x, y1, x, y2)
title('Otsikko');
legend('y = sin(x)', 'y = cos(x)');
xlabel('teksti x-akselin alle');
ylabel('teksti y-akselin viereen');
axis([-5, 5, -1.5, 1.5]);
set(gca, 'FontSize', 15);           % gca viittaa tämänhetkisen kuvan akseleihin
set(gca, 'FontWeight', 'bold');

3D-grafiikka

Analogisesti 2D-tapauksen kanssa plot3-komento piirtää annetut pisteet kolmiulotteiseen avaruuteen.

x = [-1:0.5:1]; y = [2 2 2 3 3]; z = [1:5];
plot3(x, y, z);
xlabel('x'); ylabel('y'); zlabel('z');
grid;
axis([-2,2,2,3,0,5]);

Avaruuskäyrän piirto onnistuu myös samoin:

fplot3(@(t) sin(t), @(t) cos(t), @(t) t, [0, 8*pi]);

Tarkastellaan seuraavaksi pinnan (eli funktion tasolta suoralle) piirtämistä. Tämä tapahtuu surf-komennolla. Katsotaanpa kuitenkin ensin komentoa meshgrid. Tämä komento ottaa parametreikseen kaksi vektoria x ja y ja muodostaa näistä kaksi matriisia X ja Y, jotka vastaavat vektoreiden x ja y määräämiä hilapisteitä:

x = [0:4]; y = x;
[X, Y] = meshgrid(x, y)
X =

     0     1     2     3     4
     0     1     2     3     4
     0     1     2     3     4
     0     1     2     3     4
     0     1     2     3     4


Y =

     0     0     0     0     0
     1     1     1     1     1
     2     2     2     2     2
     3     3     3     3     3
     4     4     4     4     4

Nyt kohdassa (i,j) olevat alkiot X(i,j) ja Y(i,j) muodostavat alkiota x(j) ja y(i) vastaavan hilapisteen. Näissä hilapisteissä onkin helppo laskea pisteittäin jonkin funktion arvo. Kun komennolle surf annetaan hilapisteet X ja Y sekä lasketut arvot, piirtyy ruudulle hilapisteitä yhdistävät tasonkappaleet. Piirretään esimerkkinä funktion $\sin{\sqrt{x^2+y^2}}/\sqrt{x^2+y^2}$ kuva.

x = linspace(-10,10); y = x; % x- ja y-akselien piirtoalueet samat
[X, Y] = meshgrid(x, y);
d = sqrt(X.^2 + Y.^2);
Z = sin(d)./d;
surf(X, Y, Z);

Harjoitus. Piirrä edeltävä funktio (tai oma suosikkifunktiosi). Tutustu kuvaikkunan Tools-valikon sisältöön, erityisesti kohtiin Rotate 3D ja Data Cursor.

Harjoitus. Suorita edellisillä arvoilla komento contour(X, Y, Z). Mitä saatu kuva merkitsee? Tutki myös komentoa surfc.

Muita graafisia ominaisuuksia

Tässä esitellään hyvin lyhyesti tavanomaisia datan esitysmuotoja kuten histogrammeja ja piirakkadiagrammeja. Aiheesta pidemmin tällä sivulla.

Generoidaan tuhat normaalisti jakautunutta lukua ja piirretään niistä histogrammi:

x = randn(1000, 1);
histogram(x);

Piirretään vektorin x määräämistä luvuista piirakkadiagrammi (suhteelliset osuudet kiertyvät vastapäivään klo 12 alkaen):

x = [0.3 5 2.3 0.88 1];
pie(x)

Lineaarialgebraa

Tarkastellaan ensin lineaarisen yhtälöryhmän ratkaisua. Tässä tulevat käyttöön Matlabin operaatiot / ja \ eli oikealta ja vasemmalta jakaminen. Yhtälöryhmän $Ax = y$ ratkaiseminenhan tapahtuu kertomalla vasemmalta matriisin A käänteismatriisilla. Matlabissa tämä hoidetaan "jakamalla" vasemmalta matriisilla A eli

A = [1 2 3; 0 -4 1; 0 3 -1];
y = [1 2 3]';                % pystyvektori
x = A\y
A*x                          % tarkistus
x =

    65
    -5
   -18


ans =

     1
     2
     3

Luonnollisesti oikealta jakaminen ratkaisee yhtälöryhmän $xA = y$. Tämä jakamisoperaatio ei ole Matlabissa sama asia kuin käänteismatriisilla kertominen, vaan operaatio liittyy nimenomaan yhtälöryhmän numeeriseen ratkoisemiseen, ks. mldivide ja mrdivide. Jakamisoperaatio antaa ratkaisun jopa siinä tapauksessa, että yhtälöryhmässä on enemmän yhtälöitä kuin muuttujia:

A = [1 2 3; 0 -4 1; 0 3 -1; 0 2 1];
y = [1 2 3 4]';
x = A\y
x =

   -6.7419
    0.6774
    2.1290

Tässä tapauksessa Matlab etsii ratkaisun pienimmän neliösumman mielessä.

Kerroinmatriisin tulisi jaettaessa luonnollisesti olla säännöllinen. Katsotaanpa mitä tapahtuu epäsäännöllisessä tapauksessa:

A = [1 -2 4; 0 0 6; 0 0 1];
det(A)                      % A:n determinantti on 0.
y = [1 2 3]';
A\y
ans =

     0

Warning: Matrix is singular to working precision. 

ans =

  -Inf
  -Inf
     3

Kuten tulosteesta näkyy, Matlab antaa älyttömän tuloksen, joka viittaa siihen, että yhtälöryhmällä ei ole lainkaan ratkaisua. Jos yhtälöryhmällä olisi ratkaisu, niin se kuuluisi matriisin A pystyriviavaruuteen. Lisäämällä vektori y matriisin A pystyriviksi huomataan, että aste kasvaa:

rank(A)
rank([A y])
ans =

     2


ans =

     3

Täten vektori y ei kuulu matriisin pystyriviavaruuteen eli yhtälöryhmä on ristiriitainen. Jos yhtälöryhmällä $Ax = y$ on ratkaisu, kun A ei ole säännöllinen, niin ratkaisun löytää Mooren-Penrosen käänteismatriisin avulla:

A = [1 3 7; -1 4 4; 1 10 18];
y = [5 2 12]';
x = pinv(A)*y
A*x
x =

    0.3850
   -0.1103
    0.7066


ans =

    5.0000
    2.0000
   12.0000

Mikäli tarkkaa ratkaisua ei ole, antaa yo. menetelmä ratkaisun pienimmän neliösumman mielessä. Lue tarkemmin lineaaristen yhtälöryhmien ratkaisusta täältä.

Matriisin käänteismatriisin saa etsittyä komennolla inv:

A = [1 1; 0 2];
B = [1 0; 5 2];
B*inv(A)
B/A             % myös tämä käy
ans*A           % tarkistus
ans =

    1.0000   -0.5000
    5.0000   -1.5000


ans =

    1.0000   -0.5000
    5.0000   -1.5000


ans =

     1     0
     5     2

Edellä esiteltiinkin jo komennot det ja rank, jotka laskevat matriisin determinantin ja asteen. Ominaisarvot- ja vektorit saa etsittyä komennolla eig:

A = [0 -6 -1; 6 2 -16; -5 20 -10];
eig(A)                             % etsii vain ominaisarvot (nopeampi)
[V, D] = eig(A)                    % etsii ominaisvektorit (tuloksen pystyrivit) ja ominaisarvot diagonaalimatriisissa
V*D*inv(V)                         % tämä antaa A:n (diagonaaliesitys)
ans =

  -3.0710 + 0.0000i
  -2.4645 +17.6008i
  -2.4645 -17.6008i


V =

  -0.8326 + 0.0000i   0.2003 - 0.1394i   0.2003 + 0.1394i
  -0.3553 + 0.0000i  -0.2110 - 0.6447i  -0.2110 + 0.6447i
  -0.4248 + 0.0000i  -0.6930 + 0.0000i  -0.6930 + 0.0000i


D =

  -3.0710 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i  -2.4645 +17.6008i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i  -2.4645 -17.6008i


ans =

   0.0000 + 0.0000i  -6.0000 + 0.0000i  -1.0000 + 0.0000i
   6.0000 - 0.0000i   2.0000 + 0.0000i -16.0000 + 0.0000i
  -5.0000 - 0.0000i  20.0000 + 0.0000i -10.0000 + 0.0000i

Jos matriisi ei ole diagonalisoituva, sen Schurin esityksen (unitaarinen konjugointi yläkolmiomatriisiksi) saadaan komennolla schur:

A = [13 8 8; -1 7 -2; -1 -2 7];
[U, T] = schur(A)
U =

    0.9428   -0.3333         0
   -0.2357   -0.6667   -0.7071
   -0.2357   -0.6667    0.7071


T =

    9.0000  -12.7279   -0.0000
         0    9.0000   -0.0000
         0         0    9.0000

Muita tyypillisiä matriisin esitysmuotoja esitellään sivulla Matrix Decomposition.

Numeerinen integrointi ja derivointi

Komento diff palauttaa syötevektorin alkioiden erotukset:

x = [4 3 0 -1 5];
diff(x)
ans =

    -1    -3    -1     6

Komennon avulla voi siis naiviisti laskea numeerisen derivaatan:

step = 0.1;                      % askelpituus
x = 0:step:5;
f = x.^2 + 3*x;
df = diff(f)/step;               % erotusosamäärä
plot(x, f, x(1:length(df)), df); % vektori df on yhtä alkiota lyhyempi

Helpompaa on käyttää gradient-komentoa, joka osaa laskea myös useampiulotteisen funktion gradientin.

gradient(f, step);

Numeerinen integrointi tapahtuu integrate-komennolla. Tämä komento on saatavilla vasta Matlabin versiosta R2012A lähtien. Vanhan version tapauksessa pitää käyttää alla olevaa trapz-komentoa.

integral(@(x) x.^2, 0, 4)
ans =

   21.3333

Harjoitus. Tutustu komentoon integral2 ja laske funktion $\sin{\sqrt{x^2+y^2}}/\sqrt{x^2+y^2}$ integraali alueessa $[-1, 1] \times [0, 2 \cdot \pi]$.

Integral-komento osaa laskea numeerisen integraalin vain Matlab-funktiolle - se ei kelpuuta datavektoria. Tuntemattoman funktion numeerisen integraalin voi datapisteiden avulla etsiä komennolla trapz (puolisuunnikassääntö). Lasketaan seuraavaksi tällä tavalla edeltävän funktion $x^2$ numeerinen integraali välillä $[0,4]$.

step = 0.1;
x = 0:step:4;
f = x.^2;     % generoidaan data
trapz(x, f)
ans =

   21.3400

Differentiaaliyhtälöiden numeerinen ratkaiseminen

Seuraavaksi katsotaan miten Matlabilla voi ratkaista numeerisesti differentiaaliyhtälöitä ja differentiaaliyhtälöryhmiä. Ratkaisuun käytetty komento on tässä ode45, jota suositellaan ensimmäiseksi kokeiltavaksi. Matlabissa on paljon muitakin ratkaisijoita, jota saattavat tulla tarpeen riippuen ratkaistavan differentiaaliyhtälön tyypistä, ks. Choose and ODE Solver. Matlab ratkaisee periaatteessa vain muotoa $y' = f(t, y)$ olevia differentiaaliyhtälöitä, korkeamman kertaluvun yhtälöt tulee muunnoksella saattaa tähän muotoon, ks. differentiaaliyhtälöiden kurssimateriaali tai edeltävä linkki. Komento ode45 ottaa parametreikseen yo. yhtälön oikean puolen määrittävän funktion, integrointivälin ja alkuarvon. Mikäli halutaan etsiä numeerinen ratkaisu differentiaaliyhtälölle $y' = 3ty$ aikavälillä $[0,1]$ ja alkuarvolla $y(0) = 1$, tulee siis kirjoittaa

[t, y] = ode45(@(t, y) 3*t*y, [0 1], 1);

Funktio palauttaa kaksi vektoria: ratkaisun y ajan t funktiona. Ratkaisukäyrän voi siis piirtää välittömästi:

plot(t, y, '-o')

Kysessä olevan differentiaaliyhtälön yleinen ratkaisu on $y = C e^{3t^2/2}$. Tässä tapauksessa Matlabin tuloste näyttää siis oikeanlaiselta. Komento ode45 käyttää ratkaisukäyrän etsintään Rungen-Kuttan menetelmää.

Harjoitus. Piirrä differentiaaliyhtälön $y' = \frac{A}{B}ty$ ratkaisukäyriä lukujen $A$ ja $B$ eri arvoilla. Käytä ohjelmoinnissa parametrista riippuvaa anonyymiä funktiota.

Tarkastellaan sitten vielä differentiaaliyhtälöryhmän ratkaisua ja käsitellään esimerkinomaisesti toisen kertaluvun differentiaaliyhtälön $y'' = -sin(y)$ ratkaisua. Tavanomaisella muunnoksella $y = y_1$, $y_2 = y_1'$ saadaan ratkaistavaksi ensimmäisen kertaluvun ryhmä $y_1' = y_2$, $y_2' = -sin(y_1)$. Luodaan tätä ryhmää varten seuraava funktio (jälleen voisi kirjoittaa anonyyminkin funktion):

function r = dy(t, y)
  % Huomaa, että parametri y on kahden alkion vektori, koska ryhmässä on kaksi
  % yhtälöä. Komento ode45 vaatii, että palautusarvo on pystyvektori, joten
  % funktiomäärittelyssä ei voi lukea esim. 'function [y1, y2] = dy(t, y)',
  % sillä tällöin funktio palauttaa vaakavektorin.
  r = zeros(2, 1);
  r(1) = y(2);
  r(2) = -sin(y(1));
end

Ratkaisu saadaan analogisesti:

[t, y] = ode45(@dy, [0 5], [0 1]);
plot(t, y(:,1), t, y(:,2));

Funktion minimoinnista

Etsitään sopivia funktioita minimoimiseen optimization toolboxista. Tämä sisältö löytyy dokumentoinnista tai avaamalla helpin ja kirjoittamalla "Optimization toolbox".

Yritetään löytää funktion $sin(x)+cos(x)$ minimiarvo. Työkalupakin funktioiden läpikäynnin jälkeen voidaan todeta, että ainakin fminsearch ja fminbnd -funktiot soveltuvat tehtävän ratkaisemiseen. Dokumentoinnista selviää funktiota fminsearch kutsutaan seuraavasti:

[x, fval, exitflag, output] = fminsearch(fun, x0, options);

Tässä kohdefunktio on ensimmäinen parametri. Toinen parametri on aloituspiste. Parametrin options avulla voi määritellä esimerkiksi paljonko tietoa tulostetaan.

opt = optimset('Display', 'iter'); % Tulostetaan lisätietoa iteraatioista.
sincossum = @(x) sin(x)+cos(x);
[x, fval, eflag, output] = fminsearch(sincossum, 0, opt);
sincossum(x)
output.funcCount
 
 Iteration   Func-count     min f(x)         Procedure
     0            1                1         
     1            2                1         initial simplex
     2            4           0.9995         expand
     3            6         0.998499         expand
     4            8         0.996494         expand
     5           10         0.992472         expand
     6           12          0.98438         expand
     7           14         0.968009         expand
     8           16         0.934527         expand
     9           18         0.864728         expand
    10           20         0.714808         expand
    11           22         0.382525         expand
    12           24        -0.333554         expand
    13           26         -1.34737         expand
    14           28         -1.38509         contract outside
    15           30         -1.41225         contract inside
    16           32         -1.41225         contract inside
    17           34         -1.41412         contract inside
    18           36         -1.41412         contract inside
    19           38          -1.4142         contract inside
    20           40         -1.41421         contract inside
    21           42         -1.41421         contract inside
    22           44         -1.41421         contract inside
    23           46         -1.41421         contract inside
    24           48         -1.41421         contract inside
    25           50         -1.41421         contract inside
    26           52         -1.41421         contract inside
    27           54         -1.41421         contract inside
 
Optimization terminated:
 the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-04 
 and F(X) satisfies the convergence criteria using OPTIONS.TolFun of 1.000000e-04 


ans =

   -1.4142


ans =

    54

Edellisessä koodissa muuttujien opt ja output luokka on struct. Oleellisesti se on nimetty listaus alkioita. Edellä output.funcCount tulostaa ne muuttujan output alkiot, jotka ovat nimen funcCount alla. Tällä kertaa alkioita oli yksi, mutta yleisesti niitä voi olla useita, jolloin yksittäiseen alkioon viittaaminen käy kuten alla.

structure = struct('nimet', {'x','y';'z','w'}, 'arvot', {1,2;3,4});
class(structure)
structure.nimet
structure(2).nimet
ans =

struct


ans =

x


ans =

z


ans =

y


ans =

w


ans =

z

Funktiota fminbound kutsutaan lähes samalla tavalla kuin funktiota fminsearch. Tässä aloituspisteen sijaan tarvitaan väli, josta minimiä etsitään.

[x, fval, eflag, output]=fminbnd(sincossum, 0,2*pi, opt);
sincossum(x)
output.funcCount
 
 Func-count     x          f(x)         Procedure
    1        2.39996   -0.0618786        initial
    2        3.88322     -1.41286        golden
    3        4.79993    -0.908745        golden
    4        3.88982     -1.41324        parabolic
    5        3.92941     -1.41421        parabolic
    6          3.927     -1.41421        parabolic
    7        3.92696     -1.41421        parabolic
    8        3.92703     -1.41421        parabolic
 
Optimization terminated:
 the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-04 


ans =

   -1.4142


ans =

     8

Funktion nollakohtien hausta

Yhden muuttujan funktion erään nollakohdan hakeminen onnistuu komennolla fzero. Syötteenä annetaan viittaus funktioon ja aloituspiste, jonka lähistöltä nollakohtaa etsitään. Tätä komentoa käsiteltiin aiemmin anonyymien funktioiden yhteydessä.

fsini = @(x) sin(x);
[x, fval, ef, output] = fzero(fsini, 2)
fsini(x)                                % tarkistus
%
x =

    3.1416


fval =

   1.2246e-16


ef =

     1


output = 

    intervaliterations: 10
            iterations: 6
             funcCount: 27
             algorithm: 'bisection, interpolation'
               message: 'Zero found in the interval [0.72, 3.28]'


ans =

   1.2246e-16

Funktiolle voi antaa aloituspisteen sijaan aloitusvälin, jolta nollakohta löytyy, ks. komennon ohjeet.

Funktio fzero etsii vain yhtä nollakohtaa. Polynomifunktion kaikki nollakohdat voi etsiä komennolla roots.

a = [1, 2, 3, 4, 5, 6]; % Polynomin kertoimet, x^5 + 2*x^4 + ... + 6
x = roots(a)
polyval(a, x)      % Laske polynomin arvo löydetyissä pisteissä.
x =

   0.5517 + 1.2533i
   0.5517 - 1.2533i
  -1.4918 + 0.0000i
  -0.8058 + 1.2229i
  -0.8058 - 1.2229i


ans =

   1.0e-13 *

   0.3464 + 0.1821i
   0.3464 - 0.1821i
  -0.1332 + 0.0000i
   0.1421 + 0.1998i
   0.1421 - 0.1998i

Tehtävä. Etsi edellisen polynomin jokin nollakohta komennolla fzero.

Usean funktion yhteisen nollakohdan hakeminen eli epälineaarisen yhtälöryhmän ratkaisuun sopii funktio fsolve.

nlineq = @(x) [sum(x.^2)-4, [-1,2]*transpose(x.^2)];
opt = optimset('Display', 'iter');
x = fsolve(nlineq, [1, 1], opt)
nlineq(x)                                            % tarkistus
%
                                         Norm of      First-order   Trust-region
 Iteration  Func-count     f(x)          step         optimality    radius
     0          3               5                             6               1
     1          6        0.929784       0.849837           4.99               1
     2          9      0.00256326       0.189768          0.235               1
     3         12     2.83323e-08        0.01091       0.000777               1
     4         15     3.53384e-18    3.64426e-05       8.68e-09               1

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.




x =

    1.6330    1.1547


ans =

   1.0e-08 *

    0.1329   -0.1329

Datan sovituksesta

Tutkitaan seuraavaksi miten annettuun pisteistöön voidaan sovittaa käyrä. Sopivia menetelmiä voi etsiä hakukoneella tai kirjoittamalla tulkkiin lookfor fit (tämä listaa komennoit, joiden kuvaus sisältää sanan "fit"). Tulostetuista komennoista ainakin polyfit ja lsqcurvefit näyttävät sopivan käyrän sovitukseen. Generoidaan data ja testataan näitä kahta.

x = 5*rand(200, 1);                       % 200 pistettä tasaisesta jakaumasta väliltä [0,5].
y = 2*x.^2 - 5*x + 1 + 0.5*randn(200, 1); % 2. asteen polynomi 2x^2-5x+1 ja kohina.
pcoef = polyfit(x,y,2)                    % sovitetaan toisen asteen polynomia, saadaan kertoimet
xsorted = sort(x);                        % järjestetään data piirtoa varten, Matlab piirtää
                                          % yhdistävät viivat datapisteiden antamassa järjestyksessä
ysorted = polyval(pcoef, xsorted);        % evaluoidaan kertoimien määräämä polynomi pisteissä x
plot(xsorted, ysorted, '-', x, y, '.');
title('Datan sovitus polyfit-komennolla', 'Color', 'k','Fontsize', 15);
plottext=sprintf('Sovitettu käyrä: \ny=%3.2fx^{2}%3.2fx+%3.2f', pcoef);
text(0.5,5,plottext);
pcoef =

    1.9834   -4.8765    0.8545

Tehdään sama sovitus käyttäen lsqcurvefit-funktiota.

fitf = @(a,x) a(1)*x.^2 + a(2)*x + a(3);
pcoef = lsqcurvefit(fitf, [0 0 0], x, y)
Local minimum found.

Optimization completed because the size of the gradient is less than
the default value of the optimality tolerance.




pcoef =

    1.9834   -4.8765    0.8545

Havaitaan, että funktio löysi täsmälleen saman polynomin kuin polyfit. Näin pitääkin olla, sillä molemmat käyttävät pienimmän neliösumman menetelmää. Sovitetaan seuraavaksi uuteen dataan käyrää a+b*exp(c*x). Funktiota polyfit ei voi suoraan käyttää, mutta funktion lsqcurvefit käyttö on yhtä suoraviivaista kuin edellä. Tutkitaan samalla mitä muuta lsqcurvefit-funktiosta saa irti.

x = 5*rand(200, 1);
y = 2*exp(0.4*x) + 1 + 0.5*randn(200, 1);
fitf = @(a, x) a(1) + a(2)*exp(a(3)*x);
opt = optimset('display', 'iter');
[fcoef,resnorm,res,exitflag,output] = lsqcurvefit(fitf, [0,0,0], x, y, [], [], opt);
xsorted = sort(x);
ysorted = fcoef(1) + fcoef(2)*exp(fcoef(3) * xsorted);
plot(xsorted, ysorted, '-', x, y, '.');
title('Datan sovitus lsqcurvefit-komennolla' ,'Fontsize', 15);
plottext=sprintf('Sovitettu käyrä: \ny=%3.2f+%3.2fe^{%3.2fx}', fcoef);
text(0.5, 8, plottext);
figure
resnorm
sum(res.^2)
plot(x, res, '.');
title('Residuaalit', 'Fontsize', 15);
exitflag
output.funcCount % funktioevaluointien määrä
output.algorithm % käytetty algoritmi
                                         Norm of      First-order 
 Iteration  Func-count     f(x)          step          optimality
     0          4         13077.9                      1.45e+03
     1          8         2552.35         5.1297        3.6e+03      
     2         12         715.011        4.10192       4.55e+03      
     3         16         94.0662         1.3572       2.25e+03      
     4         20         70.6755       0.779054       1.98e+03      
     5         24         41.1678       0.116167           51.7      
     6         28         41.1503     0.00100448          0.034      
     7         32         41.1503    3.64641e-05       6.77e-06      

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the default value of the function tolerance.




resnorm =

   41.1503


ans =

   41.1503


exitflag =

     3


ans =

    32


ans =

trust-region-reflective

Funktiolla lsqcurvefit voi sovittaa myös usean muuttujan funktioita. Konvergointi voi olla hidasta ja iteraatioiden oletusarvoja voi joutua kasvattamaan.

w = 5*rand(200,1);
z = 5*rand(200,1);
y = 2*exp(0.4*w) + exp(0.2*z) + 1 + 0.1*randn(200,1);
fitf = @(a,x) a(1) + a(2)*exp(a(3)*x(:,1)) + a(4)*exp(a(5)*x(:,2));
fcoef = lsqcurvefit(fitf, [0,0,0,0,0], [w,z], y);
opt = optimset('MaxFunEvals', 30000, 'MaxIter', 5000);
[fcoef,resnorm,res,ef,output] = lsqcurvefit(fitf, [0,0,0,0,0], [w,z], y, [], [], opt);
output.funcCount
output.iterations
Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the default value of the function tolerance.




Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the default value of the function tolerance.




ans =

   102


ans =

    16

Havaintoaineiston mallintamista varten kannatta tutustua Statistics toolbox:n tarjoamiin metodeihin, kuten glmfit ja mvregress. Ratkaistaan ensimmäinen esimerkki funktiolla glmfit.

x = 5*rand(200,1);
y = 2*x.^2 - 5*x + 1 + 0.5*randn(200,1);
X = [x, x.^2];                           % Mallimatriisi. Vakiotermiä ei tarvitse kirjoittaa.
pcoef = glmfit(X, y, 'normal')
xsorted = sort(x);
ysorted = polyval(flipud(pcoef), xsorted);
plot(xsorted, ysorted, '-', x, y, '.');
title('Datan sovitus glmfit-komennolla', 'Color', 'k','Fontsize', 15);
plottext=sprintf('Sovitetty käyrä: \ny=%3.2fx^{2}%3.2fx+%3.2f', flipud(pcoef));
text(0.5,5,plottext);
pcoef =

    1.1702
   -5.2219
    2.0432

Interpoloinnista

Interpolointi on hyvin läheistä sukua käyrän sovitukselle. Interpoloinnissa on tarkoitus estimoida funktion arvoja kahden havaintopisteen väliltä. Interpoloiva käyrä kulkee jokaisen havaintopisteen kautta. Tutkitaan seuraavaa dataa.

x = 1:8; y = exp(0.2*x) + 0.2*transpose(rand(8, 1));

Naivi tapa etsiä interpoloiva funktio on sovittaa riittävän korkea-asteinen polynomi. Korkea-asteinen polynomi saataa kuitenkin heilahdella suuresti pisteiden välillä, minkä vuoksi parempi ratkaisu on käyttää splinejä.

coef = polyfit(x, y, 7);
xx = 0.6:0.2:8.4;
yy = polyval(coef, xx);
plot(x, y, 'k.', xx, yy, 'b-');
hold;
yy = interp1(x, y, xx, 'spline'); % Löydetty komennolla lookfor interpolation.
plot(xx, yy, 'r-', x, y, '.');
legend('data', 'polynomi', 'splini', 'Location', 'NorthWest');
title('Interpolation', 'Color', 'k','Fontsize', 15);
hold off;
Current plot held

Tehtävä. Tee edellä oleva spline-interpolointi käyttäen fit-funktiota.