torstai 4. toukokuuta 2017

Funktion kaikkien paikalliset ääriarvojen haku Matlab-työkaluin

Lopullinen julkaisupäivä: 15.05.2017
Matlabin "Live editorilla" tuotettu, lyhennetty versio kirjoituksesta
(Lopussa mainittu "kaunis dokumentti").

Tiivistelmä:

Simo Kivelän blogissa huhtikuussa 2017 käsiteltiin tämän kevään lyhyen matematiikan ylioppilastehtävää, jossa "Kalle ja Leena tekevät biofysiikan perusteiden harjoitustyötä".  
Tässä tarina.

Yhdyn Simon mielipiteisiin. Ryhdyin kehittelemään diskreetin matematiikan ajattelutapoja, joiden kautta päästäisiin käsittelemään järkevällä tavalla tällaisia minimointi/maksimointi-tehtäviä. Kun nykyisessä signaalinkäsittelyssä jatkuvat funktiot digitoidaan, eli diskretoidaan riittävän hienojakoisesti, on ehkä perusteltua lähestyä koko problematiikkaa diskreeteistä funktioista lähtien. Tällä tavalla saadaan yhtenäisesti käsitellyksi diskreettien ja "sileiden" funktioiden tapaukset ja samantien numeeriset algoritmit ääriarvojen laskemiseen.

Matemaattisia pohjatietoja kirjoituksen ymmärtämiseksi tarvitaan varsin vähän. Lukion lyhyen linjan 1. luokan tiedot riittävät käsittääkseni vallan mainiosti.

Tarkoitus on kirjoittaa ensin tästä, ja jatkaa seuraavaksi käyrän sovittamisesta diskreettiin dataan, mikä on toinen osa yllä mainittua ylioppilastehtävää.

Erinomainen ohjelmisto näihin tarkoituksiin on Matlab, jonka ilmaisella "kloonilla" Octavella voidaan tehdä kaikki tässä esitettävät jutut yhtä hyvin. Editointiin ja julkaisuun tarvittavat työkalut ovat Matlab:ssa kehittyneempiä, mutta Octaven vaatimattomampi käyttöympäristö on loppujen lopuksi pieni haitta verrattuna siihen oivalliseen ajattelutapaan, johon korkean tason numeerinen vektorikieli antaa välineet. Octavella kehitetyt ohjelmat, m-tiedostot, toimivat sellaisenaan suoraan Matlab:ssa.

Kirjoitus on suunnattu myös lukijalle, jolla ei ole Matlab/Octave-kokemusta. Siksi annan yksityiskohtaisia opastustietoiskuja ohjelman käyttöön yksinkertaisin esimerkein, jotka Matlab-osaaja voi sivuuttaa. Ehkä antoisinta on lukea tarinaa avaamalla Matlab/Octave, ja suorittamalla samalla omia kokeiluja. Muistutan, että Octaven voi asentaa ohjelman webbisivulta eri käyttöjärjestelmissä itselleen, Ubuntu-Linux:ssa yksinkertaisimmin suoraan komennolla:
sudo apt-get install octave .

Matematiikan opetus ja ohjelmointi: Valitsemalla riittävän korkean tason kieli, jossa on matemaattisia käsitteitä ja sovelluksia tukevat rakenteet ja funktiot, voitaisiin matematiikan opetukseen tuoda ohjelmointiosuutta, joka tukisi itse matematiikan oppimista ja soveltamista ja samalla opettaisi ohjelmoinnin perusteita. Octave/Matlab voisi olla varteenotettava ehdokas.

Diskreetin ja sileän funktion paikalliset ääriarvot

 

Lähdetään liikkeelle diskreetissä pistejoukossa määritellystä funktiosta. Selvitetään juurta jaksain periaatteet, joilla päästään käsiksi ääriarvojen laskentaan, tarvitsematta mitään pohjatietoja asiasta. Kehitellään Matlab/Octave-työkalut ohjelmien perusominaisuuksia ja ajattelutapoja hyödyntäen, joilla periaatteet saadaan muutetuksi laskenta-algoritmiksi. Kun sileää (jatkuvasti derivoituvaa) funktiota f "digitoidaan", eli lasketaan arvoja rittävän tiheästi diskretoiden, saadaan samalla ohjelmalla lasketuksi halutulla tarkkuudella likiarvoja f:n ääriarvoille.

Kehitetään esimerkki, jossa on paljon sekä piikikkäitä että sileitä ääriarvoja ja sovelletaan kehittämäämme ohjelmaa.  Sen laskemia tuloksia alkuarvauksina käyttäen tarkennetaan tuloksia Matlab:n valmiin funktion fminsearch  tuottamiin, ja koetaan muutama yllätyskin.

Lopuksi valaistaan mahdollisuutta rinnakkaislaskentaan, josta tarkemmin tulevissa blogikirjoituksissa.

Matlab/Octave-oppia


Kirjoitus on suunnattu sellaiselle lukijalle, joka ei ole Matlab:ia aiemmin kohdannut. Siksi pyrin valaisemaan jokaista ohjelman käyttöä Matlabin perusteita opettavin selityksin ja esimerkein.
Matlab-istunnon muoto:
>> format compact   % Komento, jolla turhat välit tiivistetään.
>> 1+1              % Komento, ei sijoitusta muuttujaan.
ans =               % Ohjelman vastaus
     2              % Ohjelman vastaus
>> kaksi=1+1        % Komento, sijoitus muuttujaan
kaksi =             % Vastaus
     2              % Vastaus

Siis komentoa edeltää ohjelman kehote ( >> ) .
Jos seuraat juonta Matlab:n/Octaven äärellä, voit leikkaus/liimaus-tyyliin kopioida rivejä istuntoosi. Tietenkään ohjelman kehotetta ei kopioida. Suositeltavampa on kopioda kokonaisuus tekstieditoriin ja poistaa kehotteet ja tulosteet. Tiedosto talletetaan vaikkapa nimelle opiaariarvot.m .
Tiedoston komennot voi suorittaa kokonaan tai osissa (kts. oppaista "skriptit", toki myös copy/paste-tyyliin). Työskentelysuosituksia on tuossa lyhyessä oppaassa ja muissa. Kts. myös tämän kirjoituksen viimeistä kappaletta.


Perustapaukset ääriarvoille


Lasketaan  ja piirretään Matlab:lla

>> x=[0 pi/2 4 5 2*pi]
x =
         0    1.5708    4.0000    5.0000    6.2832
>> y=sin(x)
y =
         0    1.0000   -0.7568   -0.9589   -0.0000
>> plot(x,y,'.-.',x,y,'or');grid on

Miten haetaan min ja max diskreetille funktiolle ?


Funktion kuvaajan muodostavat punaiset 'rinkulapisteet'.
>> [x;y]  % Vektorit x ja y allekkain:
ans =
         0    1.5708    4.0000    5.0000    6.2832
         0    1.0000   -0.7568   -0.9589   -0.0000
 Havainnot:
  •  Välillä $(x_1,x_2)$ funktio kasvaa, eli $y_2 - y_1 > 0$ 
  •  Välillä $(x_2,x_3)$ funktio pienenee, eli $y_3 - y_2 < 0$.
SIIS: $x_2$ on maksimikohta.
  •  Välillä $(x_2,x_3)$ funktio pienenee, samoin välillä $(x_3,x_4)$. SIIS:  $x_3$ ei ole min- eikä maxkohta.
  •  Välillä $(x_3,x_4)$ funktio pienenee ja välillä $(x_4,x_5)$ kasvaa. SIIS:  $x_4$ on minimkohta.
Lasketaan Matlab:lla (Octavella):

Avaimen onneen tarjoaa funktio diff, joka laskee argumenttivektorin peräkkäisten komponenttien erotusvektorin.
Ennenkuin jatketaan esimerkkiämme, katsotaan, miten diff toimii.
>> v=[-1 1 1 3 0]
v =
    -1     1     1     3     0
>> diff(v)
ans =
     2     0     2    -3
Jatketaan juonta, jossa siis
>> x=[0 pi/2 4 5 2*pi];  % Puolipiste estää tulostuksen
>> y=sin(x)      
   y =
         0    1.0000   -0.7568   -0.9589   -0.0000 
>> diff(y)              
ans =
    1.0000   -1.7568   -0.2021    0.9589
Erotukset mittaavat x-vektorin sisäpisteissä $x_2,x_3,x_4$ tapahtuvia muutoksia.
>> [x;y]
ans =
         0    1.5708    4.0000    5.0000    6.2832
         0    1.0000   -0.7568   -0.9589   -0.0000
diff(y) 
ans =
    1.0000   -1.7568   -0.2021    0.9589
Luemme diff(y)-vektorista sen, mikä näkyy y-vektorista (ja kuvasta)
  • Pisteessä $x_2$ erotus muuttuu $+ \rightarrow -$, siis max-piste.
  • Pisteessä $x_3$ $- \rightarrow -$ , siis ei ääriarvoa.
  • Pisteessä $x_4$ $- \rightarrow +$, siis min-piste.
Huomaa, että vektori diff(y)./diff(x) on erotusosamäärävektori, se edustaa siten diskreetin funktion derivaattaa. Kun selvitämme "derivaatan" merkinvaihteluita, ei tarvitse jakaa, koska diff(x) koostuu positiivisista luvuista, useissa tapauksissa vakioaskelista $\Delta x$.
(Miksi yllä jakolaskussa (./ ) on piste ? No tässähän jaetaan kaksi vektoria keskenään. Piste viittaa vastinalkioittain ("pisteittäin") tapahtuvaan jakoon, palataan kohta lähemmin.)
Olemme kiinostuneita erotusten merkeistä, joten kannattaa soveltaa vielä sign-funktiota:
>> sdy=sign(diff(y))
sdy =
     1    -1    -1     1
>> d2y=diff(sdy)  % Erotusvektorin (merkkivektorin) erotukset:
d2y =
    -2     0     2
Tästäpä näkyy suoraan, että $x_2$ on maksimipiste, $x_3$ ei ole ääriarvo, $x_4$ on minimipiste.

Huom: Jos erotusvektorissa diff(y) esiintyy 0, niin kaksi vierekkäistä y-arvoa yhtyy, jolloin "tangentti" on vaakasuora. Tällöin vastaavat y-arvot ovat "epäaitoja" ääriarvoja. Jätämme nämä tapaukset tässä ulkopuolelle, koska ne eivät käytännössä tule esiin sileiden funktioiden, eivätkä piikikkäisten ääriarvojen yhteydessä. Sinänsä niiden mukaan ottaminen ei ole vaikeaa, mutta lienee parempi suojella esitystä liiallisilta rönsyiltä. Jääköön lukijalle (Matlab-opettavaiseksi) harjoitustehtäväksi.
Jotta voisimme kehittää koodin, joka tekee valinnan puolestamme, tarvitsemme vielä yhden tehokkaan Matlab-työkalun, loogisen valinnan.

Indeksointi ja looginen valinta Matlab:ssa


 Kts. Lyhyen Matlab-oppaan opetusta 

Periaate yksinkertaisimmillaan valaistukoon tällä, oppaasta poimitulla esimerkillä: 

>> v=-3:3
v =
    -3    -2    -1     0     1     2     3
>> v>0  % v:n jokaiselle alkiolle tehdään vertailu 0:aan.
ans =
  1×7 logical array
   0   0   0   0   1   1   1   % 1 niille, joille ehto toteutuu.
>> v(v>0)                      % Valitaan ne.
ans =
     1     2     3             % Helppoa, mahtavaa!

Jatketaan Matlab-työtä esimerkkimme parissa


>> maxind=d2y < 0
maxind =
  1×3 logical array
  1   0   0
>> minind=d2y > 0
minind =
  1×3 logical array
   0   0   1  
Sisäpisteiden $x_2,x_3,x_4$ ja vastaavien y-pisteiden muodostamista vektoreista ....
>> xin=x(2:end-1)
xin =
    1.5708    4.0000    5.0000
>> yin=y(2:end-1)
yin =
    1.0000   -0.7568   -0.9589
... valitaan loogisten vektorien maxind ja minind avulla ne, joiden kohdalla on 1 (true).
>> xmax=xin(maxind)
xmax =
    1.5708
>> ymax=yin(maxind)
ymax =
     1
>> xmin=xin(minind)
xmin =
     5
>> ymin=yin(minind)
ymin =
   -0.9589
Kaksi viime vaihetta voidaan yhdistää helposti luettavaksi "idiomiksi" tyyliin:
>> xmin=xin(d2y>0)  % Valitaan ne, joissa toinen differenssi > 0
xmin =
     5
HUOM! Edellä suoritetut operaatiot ovat kaikki vektorioperaatioita. Niin ollen ne toimivat ilman muutoksia miten pitkille x- ja y-vektoreille tahansa, ja antavat kaikki (x,y)-vektoriparin määrittelemän diskreetin funktion max- ja min-pisteet.

Tehtävä: Määritä kaikki min-ja max-pisteet datalle:
>> x=[0 pi/2 4 5 10 15 19 23 8*pi]; % Puolipiste estää tulostuksen.
>> y=sin(x);                        % Hyödyllistä, kun dataa on paljon.
>> plot(x,y,'.-',x,y,'or');grid on;




Ohje: Tee yllä olevan mallin mukaan (kopioimalla koodirivit). Suorita sitten piirto tähän tapaan:
>> clf     % clear figure
>> plot(x,y,'.-.',xmin,ymin,'*r',xmax,ymax,'ob') % 'r'=red, 'b'=blue
>> legend('Datapisteet yhdistettyinä','minimit(red))','maksimit(blue)')

Kirjoitetaan yleispätevä Matlab-funktio (ohjelma)


Pitkän päälle koodirivien kopiointi käy turhauttavaksi. Siksipä kannattaa tässä vaiheessa kirjoittaa homman hoitava funktio, jolle annetaan argumentteina x ja y-vektorit. Funktion tulee palauttaa vektorit xmin, ymin, xmax, ymax

Jos tunnet Matlab:n perusteita oman funktiomäärittelyn verran, voit hypätä seuraavan kohdan yli.

Matlab-funktion määrittely ja käyttö


Tässä muutama johdatusviite aiheeseen:
  1. Syksyllä 2016 pitämäni kurssin luentokalvoja (englanniksi)
  2. Erään Matlab-oppaan ("YAGTOM") suomenkielinen mukaelma funktioaiheesta
  3. Ennen mainitun lyhyen oppaan tekstiä aiheesta
Otetaan tähän näkyville viitteen (1) esimerkkifunktio vielä hiukan yksinkertaistaen. Jospa haluamme laskea vektorin x alkioiden keskiarvon. Matlab-lauseke voidaan kirjoittaa seuraavan esimerkin tavoin:
>> x=1:10;
>> avg=sum(x)/length(x)  % Ei tarvita (./), koska jaetaan skalaareja,  
                            % ... eli lukuja, ei vektoreita.
Haluamme funktion "keskiarvo", jota voitaisiin kutsua millä tahansa vektorilla x tähän tapaan:
>> tulos=keskiarvo(x)
Tekstitiedostoon keskiarvo.m kirjoitetaan koodi, jonka 1. rivillä on avainsana function seuraavaan tapaan:
function y=keskiarvo(x)
% Lasketaan x-arvojen keskiarvo.
% Input: vektori x
% Output : x:n alkioiden keskiarvo
% Kutsuesimerkki: tulos=keskiarvo(1:10)
%
y=sum(x)/length(x);  % Ainoa koodirivi
Nyt voidaan komentaa Matlab-istunnossa vaikka:
>> help keskiarvo  
  Lasketaan x-arvojen keskiarvo.
  Input: vektori x
  Output : x:n alkioiden keskiarvo
  Kutsuesimerkki tulos=keskiarvo(1:10)
 
>> tulos=keskiarvo(1:10)
tulos =
    5.5000
>> keskiarvo(rand(1,100))
ans =
    0.5184
>> keskiarvo(rand(1,1000))
ans =
    0.5153
>> keskiarvo(rand(1,1000))
ans =
    0.5046
Toisinaan on tarvetta funktiolle, jolla on useita argumentteja lähtö- ja/tai tulospuolella. Voisimme haluta funktion, joka palauttaa keskiarvon lisäksi x-vektorin arvoalueen, eli vektorin [min(x),max(x)]. Edellistä muokattaisiin näin:

 function [ka,arvoalue]=ka_minmax(x)
 % Lasketaan x-arvojen keskiarvo ja arvoalue.
 % Input: vektori x
 % Out1 : x:n alkioiden keskiarvo
 % Out2 : [min(x) max(x)]
 % Kutsuesimerkki [ka,range]=keskiarvo(1:10)
 %
 ka=sum(x)/length(x);  
 arvoalue=[min(x),max(x)];
 end


HUOMAA: Funktiotiedosto on aivan tavallinen tekstitiedosto. Sen kirjoittamiseksi ei ole mitenkään pakko käyttää viitteessä (1) neuvottua klikkailua, ei edes Matlabin editoria, vaan mitä tahansa tekstieditoria. Tämä on tärkeää (ja helpottavaa) huomata erityisesti Octave-työskentelyä ajatellen.

Käyttöesimerkki:
>> x=linspace(-pi,4*pi); % Väli [-pi,4*pi] jaetaan 99 osaväliin.
>> y=x.*sin(x); % Pisteittäinen, vastinalkioittainen vektorikertolasku
>> [k,mima]=ka_minmax(y);
>> k
k =
   -0.5928
>> mima
mima =
  -11.0250    7.9160
>> plot(x,y)
>> title('x sin x')
>> hold on
>> plot(x([1 end]),[k k])
>> grid on
>> legend('x sin x','keskiarvo')
>> yrajat=mima+[-0.5 0.5]
yrajat =
  -11.5250    8.4160
>> ylim(yrajat)        % Väljennetään y-aluetta hiukan.

 

Ääriarvojuoni jatkuu



Kokoamme esitetyt palaset, funktiomäärittelysyntaksin ja kehittelemämme koodirivt yhteen. Tuloksena voimme kirjoittaa seuraavanlaisen funktion:

function [xmin,ymin,xmax,ymax] = lok_minmax(x,y)
% Etsitään diskreetillä x-akselilla annettujen y-arvojen paikalliset
% minimit ja maksimit
% Esim 1) 
% x=0:0.1:7; y=sin(5*x); [xmin,ymin,xmax,ymax]=lok_minmax(x,y)
% plot(x,y,xmin,ymin,'x');ylim([-1.1,1.1]);grid on;shg
% Esim 2)
% x=linspace(0,7);y=sin(5*x); [xmin,ymin,xmax,ymax]=lok_min(x,y)
%
% Voidaan siis käyttää myös sileän funktion ääriarvojen numeeriseen
% laskentaan.

dysign=sign(diff(y)); % Differenssivektorin (y(2)-y(1), y(3)-y(2), ...)
                      % ... "merkkivektori"
d2y=diff(dysign);     % 2. differenssi (+-2 tai 0)
xin=x(2:end-1);       % x-sisäpisteet
yin=y(2:end-1);       % vastaavat y-pisteet
xmin=xin(d2y>0);      % Poimitaan sisäpistevektorista xin ne, joissa 
                      % 2. differenssi > 0.;
ymin=yin(d2y>0);      % vastaavat y-pisteet
xmax=xin(d2y<0);      % Poimitaan sisäpistevektorista xin ne, joissa 
                      % 2. differenssi < 0.;
ymax=yin(d2y<0);      % vastaavat y-pisteet
(Jos haluat ottaa käyttöön, niin COPY/PASTE ja sijoita työhakemistoosi tiedostoksi lok_minmax.m ) tai kts. linkki viimeisessä kappaleessa.)

Nyt voidaan keskittyä rakentamiemme työkalujen tarjoamiin nautintoihin.

Lukuisa määrä sileitä ja piikikkäitä ääriarvoja



Kehitellään nopeasti värähtelevä funktio, paikoitellen sileä, toisin paikoin taas piikikäs, vaikkapa näin: $$f(x)=0.1\, x + (\cos x) |\sin(5\,x)|$$
Määritellään $f$ Matlab-funktiokahvaksi ("function handle"), kts. yllä olevat referenssit.

>> f=@(x)0.1*x+cos(x).*abs(sin(5*x))
f =
  function_handle with value:
    @(x)0.1*x+cos(x).*abs(sin(5*x))
Miksi (.*). Kts: Lyhyesti: cos(x) ja abs(sin(5*x)) ovat yhtä pitkiä vektoreita, jotka halutaan kertoa vastinalkioittain (eli "pisteittäin") keskenään. Tulossa 5*x, pistettä ei tarvita, koska 5 on skalaari.
Jatketaan Matlab-istuntoa:
>> vali=[0 10];
x=linspace(0,10,301);
plot(x,f(x));title('Piikkejä ja sileitä ääriarvoja');
grid on;ylim([-0.7 1.7]);shg  % shg:"show graphics"




Jaetaan väli [0,10] 300:aan osaväliin, joten osavälin pituus:
>> h=10/300
h =
    0.0333
(Kuva piirrettiin muutaman kokeilun jälkeen tällä "resoluutiolla", sen laadusta päätellen pisteitä on riittävästi, jotta päästään kaikkiin ääriarvoihin käsiksi.)
>> x=linspace(0,10,301);  % Kerrataan x:n määrittely.
>> y=f(x);                % Puolipiste estää tulostuksen.
x ja y ovat nyt 301:n alkion vektoreita, jotka yhdessä edustavat diskreettiä funktiota, jonka ääriarvoja haemme.
>> [xmin,ymin,xmax,ymax]=lok_minmax(x,y);
>> length(xmin)
ans =
    17
>> length(xmax)
ans =
    17
Sekä minimejä että maksimeja on 17 kpl.
Jatketaan piirtämistä edelliseen kuvaan: Piirretään minimipisteisiin punainen tähti ja maksimipisteisiin sininen rinkula.
>> hold on  % Jatka piirtoa edelliseen grafiikkaruutuun.
>> plot(xmin,ymin,'*r',xmax,ymax,'ob')
>> ylim([-0.7 1.7])  % Väljennetään y-aluetta.


Kuvan perusteella näyttää siltä, että olemme löytäneet approksimaatiot kaikille ääriarvoille. Approksimaatioiden absoluuttinen tarkkuuskin on tiedossamme: |virhe| on korkeintaan h= 0.0333
Katsotaan nyt numeerisia arvojakin, ehkä vaikka 7 ensimmäistä:
>> [xmin(1:7);ymin(1:7)]
ans =
    0.6355    1.2709    1.7391    2.2408    2.8428    3.4448    4.0134
    0.0922    0.1481    0.0623   -0.3835   -0.6686   -0.6085   -0.2023
>> format long    % Täysi näyttötarkkuus (=laskentatarkkuus)
>> [xmin(1:3);ymin(1:3)]  % Nyt vain 3 arvoa näytiksi.
ans =
   0.635451505016722   1.270903010033445   1.739130434782609
   0.092242289698521   0.148144547093171   0.062294731898389
>> format       % Takaisin oletusnäyttötarkkuuteen.


Onko kaikki löydetty, miten parannetaan tarkkuutta?


Helpoin tapa olisi hienontaa jakoa. Ensimmäiseen kysymyksen riittäisi kaksinkertaistaa kerran pari, tyyliin
>> N=300;
>> %% Toistetaan seuraavia rivejä:
>> N=2*N; x=linspace(0,10,N);y=f(x);   
>> [xmin,ymin,xmax,ymax]=lok_minmax(x,y);
>> minlkm=length(xmin);   % Muuttuuko edellisestä?
>> maxlkm=length(xmax);   %  ------- " ----------
>> minlkm                 % No katsotaan. Muuttujan arvon saa näkyviin
                          % kirjoittamalla sen nimen komentokehotteeseen.
minlkm =
    17                     
>> maxlkm                 
maxlkm =
    17
>> N                      % Tarkistetaan vielä N:n arvo
N =
   600

Pistemäärän kaksinkertaistus ei muuttanut lukumäärää. On siis hyvä syy uskoa, että kaikki ääriarvot on mainitulla tarkkuudella löydetty.
Tässä virheraja puolittui. Tarkkuutta saataisiin yksi dekadi lisää 10-kertaistamalla jakopisteiden lukumäärä (N=10*N).

Tässä voisi olla eräs luonteva lopetuskohta. Mutta luontevaa on myös jatkaa hiukan.

Tuloksien tarkentaminen kehittyneillä hakumenetelmillä


Numeerisia optimointimenetelmiä, samoin kuin nollakohdan hakumenetelmiä on kehitetty erityisesti yhden paikallisen ääriarvon tai nollakohdan laskentaan. Tällaisia ovat Fibonacci- ja kultaisen leikkauksen haku, sekä Newtonin menetelmä (ääriarvon tapauksessa derivaatan) nollakohdan etsinnässä. Yleensä nämä menetelmät keskittyvät yhden minimin hakuun, joka sijaitsee "riittävän lähellä" annettua alkuarvoa $x_0$
Matlab:ssa on valmis funktio fminsearch, joka toimii myös usean muuttujan funktioille. Se ei käytä derivaattaa, joten voidaan soveltaa myös piikikkäisiin tapauksiin. Yksinkertaisimmillaan funktion kutsu on
 >> xmin = fminsearch(Fun,x0)

Jos halutaan hakea kaikki ääriarvot jollain alueella, pitää ensin löytää sopivat alkuarvaukset. No nythän meillä on ohjelma juuri siihen tarkoitukseen. Voimme jatkaa seuraavanlaiselĺa koodilla, joka on luontevinta suorittaa skriptitiedostosta, siksi jätimme kehotteet pois.
 N=length(xmin);
 tarkennetutminimit=zeros(1,N);
 for k=1:N
   tarkennetutminimit(k)=fminsearch(f,xmin(k));
 end
 xminT=tarkennetutminimit;
 yminT=f(xminT);
Piirretään:
plot(xminT,yminT,'ok',xmax,ymax,'ob');ylim([-0.7 1.7]);
title('Rengastetut pisteet löytyvät, 3 jää saamatta')


Tutkitaanpa asiaa

Sijoitetaan allekkain kiinnostavat komponentit:
>> [xmin([7 8 9 13 14 15]);tarkennetutminimit([7 8 9 13 14 15])]
ans =
    4.0134    4.5151    5.0167    7.5251    8.0268    8.5284
    4.0296    3.4393    5.0265    9.1190    8.5254    8.5254
Ylärivillä on siis omalla ohjelmallamme saadut pisteet, joita käytetään alkuarvoina fminsearch-funktiolle. Kuten kuvasta näkyy, 3 minimiä jää saavuttamatta tässä "tarkennnuksessa". Kun katsotaan komponentteja 7 ja 8, huomataan, että alkuarvauksen kasvattaminen johtaa tuloksen pienenemiseen. Myös komponentti 13 (arvo 7.5251) suppenee "väärään" minimiin. Huomataan siis, että tarkentamalla saadaan oikeita arvoja, mutta osa suppenee samoihin arvoihin ja siten tässä tapauksessa 3 arvoa jää saavuttamatta.

Katsotaan vielä kokonaisuudessaan:
>>[xmin;tarkennetutminimit]
ans =
  Columns 1 through 7
    0.6355    1.2709    1.7391    2.2408    2.8428    3.4448    4.0134
    0.6283    1.2566    1.7324    2.2422    2.8359    3.4393    4.0296
  Columns 8 through 14
    4.5151    5.0167    5.6522    6.2876    6.9231    7.5251    8.0268
    3.4393    5.0265    5.6549    6.2832    6.9115    9.1190    8.5254
  Columns 15 through 17
    8.5284    9.1304    9.7324
    8.5254    9.1191    9.7225
>> diff(ans,1)    % Erotukset rivi-indeksin suunnassa
ans =
  Columns 1 through 7
   -0.0071   -0.0143   -0.0067    0.0014   -0.0069   -0.0055    0.0163
  Columns 8 through 14
   -1.0757    0.0098    0.0027   -0.0045   -0.0116    1.5939    0.4986
  Columns 15 through 17
   -0.0031   -0.0114   -0.0099
Virheraja näyttää pitävän paikkansa, paitsi komponenttien 8, 9, 13, 14, 15 kohdalla. MUTTA: Niistä lähtien fminsearch suppeneekin kohti "väärää" minimiä. Meidän menetelmämme on luotettavampi, joskin epätarkempi. Parannusta saadaan suorittamalla tarkennettu käsittely niille alkuarvoille, jotka fminsearch ajaa virherajan ulkopuolelle.
Jos tarkkuusvaatimus ei ole kovin suuri, voidaan tyytyä oman ohjelmamme käyttöön, jossa saadaan tarkkuutta lisää joko jakamalla koko väli tiheämmin (jokaista dekadia kohti dekadi lisää, aika tehotonta, tosin Matlab:ssa viuhuu lujaa). Parempi tapa: Jaetaan kukin osaväli (xmin(k)-h,xmin(k)+h) esim. 100:aan osaan.

Harjoitustehtävä


Tee vastaava selvitys maksimipisteille. Huomaa, että fminsearch sopii funktion $f$ maksimien etsintään tyyliin >> xmax=-fminsearch(@(x)-f(x),x0) , koska $\max f(x) = - min(-f(x))$

Alkuarvauksen valinnasta

Käsittelemme alkuarvon valintaproblematiikkaa tarkemmin seuraavissa blogikirjoituksissa funktion nollakohdan määritystehtävän ja Newtonin menetelmän yhteydessä. Senverran voidaan paljastaa, että pieni muutos alkuarvauksessa voi johtaa suppenemiseen eri ratkaisua kohti. Jos tutkitaan niiden pisteiden joukkoa, joista alkaen suppeneminen tapahtuu johonkin tiettyyn ratkaisuun, johdutaan usein fraktaaliseen rakenteeseen. Tämä näkyy erityisen havainnollisesti tutkittaessa esim. polynomifunktioiden nollakohtia kompleksitasossa.
Pääsääntönä toki on, että yleensä nopeasti suppenevat menetelmät toimivat, kun valitaan alkuarvo "riittävän läheltä" ratkaisua, joten jaon riittävä hienontaminen johtaa yleensä hyvään tulokseen.

Rinnakkaislaskentaa, parallell toobox


Viimeisen vuosikymmenen aikana on algoritmien suunnittelussa yleistynyt enenevässä määrin rinnakkaisuuden hyödyntäminen. Tämä tarkoittaa sitä, että ohjelma pyritään jakamaan toisistaan riippumattomiin osiin, jotka voidaan suorittaa samanaikaisesti. Fyysisesti sen tekee mahdolliseksi tietokonearkkitehtuuri, jossa ison muistin ja yhden tehokkaan prosessorin sijasta on jaettu muisti ja monta prosessoria. Tällaisia suuria "klustereita" on mm. Aalto-yliopiston Tritonus klusteri ja CSC-superlaskentakeskuksen Taito,
joissa kummassakin on Matlab ja seuraavassa puheena oleva parallel toolbox käytettävissä.
Matlabiin on siis saatavissa rinnakkaislaskentaan erikoistunut "parallell toolbox", jolla voidaan suorittaa tämänkaltainen rinnakkaistus hyvin helposti. Harjoittelu kotikoneella on vaivatonta, monessa tapauksessa koodia ei tarvitse muuttaa lainkaan siirryttäessä isompaan klusteriin. Kotiläppärin 2-ydinprosessori antaa vain 2 rinnakkaisprosessia, mutta se mahdollistaa alkutestaamisen, vaikkei vielä sanottavaa tehonlisäystä tuokaan.
Edellä esitetty algoritmi soveltuu erinomaisesti rinnakkaislaskentaan, kunkin alkuarvon laskemisen jälkeinen tarkentaminen jakautuu luonnollisesti riippumattomiin osiin.
Palaan aiheeseen seuraavissa blogikirjoituksissa. Viittaan tässä kurssiin, jonka pidin syksyllä 2016 yhdessä Juha Kuortin kanssa: Matlab course 2, "advanced", parallel

Optimointia ja nollakohtien hakua usemman muuttujan funktioille:

Tehokkuuden tarve lisääntyy huomattavasti, kun muuttujia tulee lisää, jolloin rinnakkaislaskennan edut moninkertaistuvat. Näistäkin aiheista lisää jatkossa.

Matlab/Octave-tiedostot, Matlab:n julkaisuvälineiden tuotokset


Viimeksi mainittu pdf-dokumentti on hiukan lyhennettynä sama, joka tässä on kirjoitettu Blogger-välineillä. Matlabin Live editor on kunnianhimoinen ja monin tavoin hieno työväline, jolla myös matemattiset kaavat voidaan editoida suoraan Latex-muodossa, kuvat tulevat automaattisesti ym. Omat pikku hankaluutensakin siinä toki on, ja suositeltavinta on monestakin syystä rakentaa Matlab-dokumentti yllä olevan esimerkin mukaisesti strukturoiduksi m-skriptiksi, ja ajaa (tai olla ajamatta) julkaisuksi. Käytön tietyistä epämukavuuksista kertovat omat korkeat versionumerot tiedostojen nimissä. Toisaalta tuloksena saatava dokumentti on kyllä vieläkin kauniimpi, kuin Bloggerilla tuotettu, vai mitä?.

maanantai 23. marraskuuta 2015

Alexander Grothendieck, Laurent Schwartz, funktionaalianalyysi,distribuutioteoria

Alexander Grothendieck


Sekaannuksen välttämiseksi siirrän viime vuoden joulukuulta peräisin olevan kirjoituksen yhteen ja samaan aktiiviseen "blogivirtaan". Huomautan samalla, että Grothendieck'n funktionaalianalyysin tutkijan uraan suuresti vaikuttanut Laurent Schartz sai Fieldsin mitalin vuonna 1950. 
Ylihuomenna sattuu Helsingin yliopiston matematiikan laitoksen kollokviosarjassa olemaan aiheeseen liittyvä esitys:
The next colloquium of the Dept. of Mathematics and Statistics will take place on next Wednesday Nov. 25th, at 4.15 pm in CK112. The guest speaker is Prof. Jesper Lützen (University of Copenhagen) Abstract:
In the fall of 1947 Laurent Schwartz visited Copenhagen and Lund invited by Harald Bohr, who was very impressed and enthusiastic about Schwartz's new theory of distributions. This was Schwartz's first longer trip abroad and he continued to have fond memories of his visit. There is a direct link from this visit to Schwartz's reception of the Fields Medal in 1950. In the talk I shall explain the circumstances around Schwartz's visit to Denmark and Sweden in the context of the development of the theory of distributions.

Grothendieckiin:


Alexander Grothendieck


Tämän päivän (7.12.2014) Hesarissa oli nekrologipalstalla  Osmo Pekosen kirjoitus otsikolla:
Alexander Grothendieck 1928 – 2014, Nero uudisti matematiikan ja katosi.

Facebookin Matematiikka-ryhmässä on ollut ihmettelyä siitä, kuka tämä tuntematon nero oikein oli.

Koska olen itse matematiikan tutkijana toiminut funktionaalianalyysin alalla,
erityisesti Grothendieck'n töistä lähtöisin olevien kysymysten parissa, ajattelin
hiukan valottaa asiaa tältä kannalta.

Ryhdyin kokeilemaan blogikirjoittelua, kun Facebookin kirjoitukset pian katoavat datavirran syövereihin.

Viitteitä

The Grothendieck Circle mourns the passing away of Alexandre Grothendieck on November 13, 2014.
With the agreement of Grothendieck's family, the work of the Circle to bring Grothendieck's unique story and writings to the public has resumed.

Muita linkkejä

Funktionaalianalyysia ja algebrallista geometriaa ym.

Useimmissa yllä olevissa viitteissä painotetaan eniten algebrallisen geometrian alan töitä, joiden perusteella AG:lle myönnettiin Fieldsin mitali, josta hän kieltäytyi. (Tulee etsimättä mieleen tapaus Grigori Perelman v. 2006).
Wikipediasta löytyy lyhyt ja ytimekäs jaottelu, jossa funktionaalianalyysin osuus myös näkyy:
On the advice of Cartan and André Weil, he moved to the University of Nancy where he wrote his dissertation under Laurent Schwartz and Jean Dieudonné in functional analysis, from 1950 to 1953.[16] At this time he was a leading expert in the theory of topological vector spaces.[17] By 1957, he set this subject aside in order to work in algebraic geometry and homological algebra.

Funktionaalianalyysi, topologiset vektoriavaruudet, topologiset tensoritulot

Näihin aiheisiin minut johdatti Kalus Vala, joka oli ranskalaisen matemaatikkokoulukunnan ja "Bourbakismin" suuri ihailija ja tuntija. Hän myös perusti seminaarin, jota itsekin vedin pitkään, ja joka jatkaa edelleen toimintaansa nuorempien sukupolvien johdolla, pisimpään toimineena nykyisenä vastuuhenkilönä Hans-Olav Tylli:

Grothendieck-turhautumaa musiikilla ilmaistuna 

Funktionaalianalyysin alalla Grothendieck'n teokset vaativat kokonaisen uuden koulukunnan selittämään hänen tuloksiaan. Ehkä huomattavin koulukunta tällä alalla oli A. Pietsch:n tutkimusryhmä Jenassa DDR:ssä.
Alan vaihdos algebrallisen geometriaan ei ainakaan helpottanut asiaa, kuten seuraava laulu kertoo.

 Grothendieck-song

(Video jatkuu esitelmällä Jean-Pierre Serre: How to write mathematics badly.)

sunnuntai 4. lokakuuta 2015

Taustaa, Matlab-perusteita

Viimeksi kuluneiden runsaan 30:n vuoden kuluessa yhtenä intohimonani on ollut matemaattisten ohjelmistojen ja korkean tason kielien tarjoamien mahdollisuuksien kokeilu ja kehittäminen matematiikan opetuksessa ja tutkimuksessa. Nykyisin aihe on taas erityisen ajankohtainen etenkin koulumaailmassa, jossa ohjelmointia ollaan kovalla ryminällä tuomassa opetukseen mukaan.

Tässä "Mattie-esittelyä" varten tekemäni kooste.

Siinä palautetaan mieleen symbolilaskennan ja numeeristen ohjelmistojen historiaa, ja mukana on laaja linkkilista omiin kursimateriaaleihini vuosien varrelta.

Nostan tässä aikajanani alkupisteeksi esiin kauan ennen Internet-aikoja Helsingin yliopiston matematiikan laitoksella aloittamaamme seminaariin liittyvän näytteen, joka toimitettiin  tuonaikaisella viestintäohjelmalla nimeltään Portacom.

                                 https://math.aalto.fi/opetus/Mattie/Blogi/HYseminaari83-84.pdf

Aloitan nyt taas Matlab:sta, samasta ohjelmasta, josta silloin  syksyllä 1983 niinikään aloitettiin!! Toden totta, näissä IT-alan  myrskyissä on jotain näin pysyvää. Toki ohjelmassa on tapahtunut valtavasti kehitystä noista alkuajoista, mutta perusrakenteet ja ajattelutapa ovat samat kuin silloin.  Matlab on tänään mitä ajankohtaisin matemaattinen ohjelma ja korkean tason ohjelmointikieli, jolla on erittäin laaja käyttäjäkunta.

Toinen pääkohde on symbolilaskentaohjelma Maple. Maplen kanssa tasavertainen ohjelma on Mathematica, jonka käsittelyn jätän mahdollisten muiden blogikirjotteluuni osallistuvien henkilöiden varaan. "Vierailukirjoittajiksi" ovat alustavasti lupautuneet ainakin Aalto-yliopiston yhteistyökumppanini Juha Kuortti ja Pekka Alestalo .

Eräs intohimoni on näiden vuosien aikana ollut numeerisen korkean tason kielen (kuten Matlab, APL) ja symbolilaskentakielen sopiva yhdistäminen ongelmia ratkottaessa.  Tähän on nyt erillisten ohjelmien "seurustelun" lisäksi tarjolla Matlab:n "symbolic toolbox", jota ohjelman valmistaja Mathworks on viime vuosina tuntuvasti kehittänyt. Mainittu intohimo näkyy "Symbolinen ja numeerinen laskenta"-nimisinä erikoiskursseina, joiden materiaaleihin on linkkejä yllä. Tämä on yksi "toolbox", jota takuuvarmasti tullaan käsittelemään blogeissamme.


Mattie - matematiikkaa tietokoneella 

https://math.aalto.fi/opetus/Mattie/ 

 

Jarmo Malisen aloitteesta ryhdyimme Eerikki Lehtisen ja Juha Kuortin kanssa rakentamaan Aalto yliopiston matematiikan laitoksella tehtäväportaalia, johon oli määrä kerätä ja edelleen kehittää tehtäviä, ohjeita, opetusmateriaalia, ratkaisuja ym. matematiikan tietokoneharjoitusten tarpeisiin ennen kaikkea sen materiaalin pohjalta, jota oli tällaisia kokeiluja harrastaneiden opettajien toimesta kehittynyt. Simo Kivelältä saatiin alkajaisiksi huomattava määrä käyttövalmiita Mathematica-tehtäviä. Jatkossa painopiste on ollut Matlab- ja Maple-tehtävillä, koska näihin liittyvää materiaalia on itselleni  kertynyt eniten.

Teknisestä toteutuksesta tässä senverran, että tehtävien näkyvillä oleva lähdekoodi on $\LaTeX$:ia, opettaja voi ottaa suoraan tehtävän, ja muokata sitä tarpeen mukaan. Tehtäviä voi selata html-sivulla, jossa matemaattiset kaavat näytetään MathJax:llä: https://www.mathjax.org/ , aivan kuten tässäkin blogissa teemme. (Edellyttää parin "loitsun" kirjoittamista blogin html-koodin alkuun.)


Näytteeksi vaikka eksponenttifunktion sarjakehitelmä:

$$ e^x = \sum_{k=0}^\infty \frac{x^k}{k!}= 1+x + \frac{x^2}{2!} + \ldots $$

Matkan varrella on syntynyt lukuisa määrä oppimateriaalia, kuten Maple-oppikirja, monia Matlab-verkkomateriaaleja, ohjelmistoja hyödyntävää luentomateriaalia, hienoja oppilastöitä niin peruskursseihin kuin erikoiskursseihin liittyen.
Näitä nostamme esiin tämän alkavan "blogimatkan" varrella, "jalostamme" sekä ajanmukaistamme. Samalla täydennämme MattieT - osiota uusilla tehtävillä. Joukkoon sijoitamme myös ns. demotehtäviä, joiden ratkaisut näkyvät kaikille, korjaamme virheitä, ja pyrimme parantamaan käytettävyyttä palautteen mukaan.

Tätä kirjoittaessa syntyi ajatus luokitella tehtäviä myös sen mukaan, soveltuuko tämä matematiikan tason perusteella koulumaailmaan. Niinpä osassa tehtäviä kuvaavaa luetteloa on attribuutti: "KOULU":
https://math.aalto.fi/opetus/Mattie/MattieT/html/contentsmlBasic.html
Näitä tulee eri aihepiirien sisällysluetteloihin, kunhan ehtii.

Blogin pyrkimyksiä 

 

Blogi voisi toimia myös eräänlaisena etäopiskelukurssina täydentäen ja laajentaen edellä mainittuja lyhytkursseja, joissa jouduimme ajoittain juoksemaan asioita läpi liian vauhdikkaasti vähäisen tuntimäärän vuoksi.
Pyrimme myös tarjoamaan yhtä näkökulmaa koulujen koodausbuumiin esittelemällä ohjelmoinnin ja matematiikan rinnakkaiseloa tehtävillä, joissa koulumatematiikan pohja riittää. 

Aloitamme alkeista, annamme runsaasti viitteitä omatoimiseen opiskeluun, sekä tehtäväehdotelmia, joiden ratkaisut laitamme aikanaan näkyville.

Blogi ei pyri olemaan johdonmukainen opasmateriaali, saatamme hyppiä aiheesta, vaikeusasteesta ja ohjelmasta toiseen. Pyrkimys on joka tapauksessa myös ohjata itseopiskelun polulla. Toisaalta tyypillisiä kirjoituksia voivat olla jotkut kekseliäät ratkaisut, ajattelutavat, tekniset ratkaisut,  ohjelmistojen uudet piirteet, "toolboxit", ohjelmisovalmistajien omista blogeista esiin nostetut ideat, uusien ohjelmistoversioiden tarjoamat tekniikat, mahdollisuudet ovat miltei rajattomat.

Osa blogiteksteistä on tarkoitus kirjoittaa englanniksi, aluksi ainakin niin, että ensin suomeksi, ja sitten sama tai vastaava käännettynä englanniksi, katsotaan, mihin kaikkeen resurssit riittävät.


Matlab:n perusteita, osa 1

 

Miten pääsen ohjelman ääreen


Aalto-yliopistossa sekä henkilökunta että oppilaat saavat Matlab:n omalle koneelleen ohjelmistojakelusta. Opiskelijoille ja oppilaitoksille on saatavissa/neuvoteltavissa edullisia lisenssejä. Sama pätee muihinkin tässä blogissa käsiteltävin ohjelmiin.

Lisäksi  saatavilla on  ilmaiset "kloonit", jotka voi ladata linkeistä

https://math.aalto.fi/opetus/Mattie/MattieO/octave.html
https://math.aalto.fi/opetus/Mattie/MattieO/scilab.html


Matlabin hieno editori julkaisuominaisuuksineen, monet toolboxit, ym. eivät ole näillä klooneilla käytettävissä, mutta huomattava määrä Matlab-työskentelyn perusteista toimii samalla tavoin, joten opetteluun ne sopivat vallan hyvin.

Mistä aloitan


 MattieO:n Matlab-sivulla on kattava kokoelma viitteitä

https://math.aalto.fi/opetus/Mattie/MattieO/matlab.html 

Erittäin suositeltava suomenkielinen opas on Timo Mäkelän
https://sites.google.com/site/laskenta/matlab 
(Samalta tekijältä on myös hyvä Maple-opas.)

Itselläni on erinäisiä html-tekstejä, ne ovat sikäli käteviä tässä yhteydessä, että voin antaa täsmälinkkejä haluamiini kohtiin ja työstää tekstejä lennossa paremmiksi (tai huonommiksi) .
https://math.aalto.fi/~apiola/matlab/opas/mini/Index.html
https://math.aalto.fi/~apiola/matlab/opas/lyhyt/
https://math.aalto.fi/~apiola/matlab/opas/mini/tutoriaali.html


Jos pidät video-opetuksesta, voit kokeilla videotutoriaaleja sivulta:
http://se.mathworks.com/academia/student_center/tutorials/launchpad.html
Suositeltavia tässä vaiheessa: (1) Getting started, (2) Writing a Matlab program
Parasta tietysti jälleen, jos Matlab/Octave on samalla auki. Muita videoita ja koodiesimerkkejä on samalla sivulla.

Videoita löytyy Google-haulla suoraan Youtubesta. Seuraavassa on samoja asioita työskentelytyylistä ja  m-tiedostoista (skriptit ja funktiot) hiukan eri esimerkein.
https://www.youtube.com/watch?v=NPd13u4i6fM

Eräs tapa aloittaa, on käydä läpi kevään 2014 kurssille tekemäni ensimmäisen luennon (Beamer)kalvot:

https://math.aalto.fi/opetus/MatOhjelmistot/2014kevat/L/MatlabPerusteet1.pdf

Niissä on mukana myös relevantteja linkkejä.

Lyhyen pohjatiedon hankkimisen jälkeen, vaikkapa viimemainittujen luentokalvojen tai joidenkin muiden yllä olevien viitteiden  perusteella voit käydä esimerkkeihin käsiksi. Parasta on tietenkin jälleen, jos käsillä on Matlab tai Octave (tai Scilab), jolloin voit tehdä sopivia omia muunnelmia.
Matlab-komennot on talletettu tekstitiedostoon, ns. m-tiedostoon, tässä tapauksessa https://math.aalto.fi/opetus/Mattie/Blogi/Matlab/ExpTaylor.m
Tämä tiedosto on ajettu Matlab:n publish-komennolla (klikkauksella) tiedostoksi
https://math.aalto.fi/opetus/Mattie/Blogi/Matlab/html/ExpTaylor.html , 
joka samalla näyttää, miten Matlab-työstä saa kauniin html-dokumentin kuvineen. Poimin komennot kuvineen tähän suoraankin, ettei jää pelkän klikkailun varaan. Tämä on samalla blogialustan käyttökokeilua, miten Matlab-julkaisuvälineet ja blogialusta saataisiin parhaiten yhteistyöhön.
Sopivia komentokokonaisuuksia voit suoraan kopioida ("copy/paste") ohjelman komentoikkunaan tai omaan m-tiedostoosi.

Esimerkki 1

Exponenttifunktion Taylorin polynomeja


Matlabin publish-työkalun generoimassa (kauniimmassa) muodossa tässä: 
https://math.aalto.fi/opetus/Mattie/Blogi/Matlab/html/ExpTaylor.html

Kuten yllä jo muistutettiin, sarjakehitelmä on tällainen:

$$ e^x = \sum_{k=0}^\infty \frac{x^k}{k!} $$

Matlab-istunto:


close all                         % Suljetaan mahd. avoinna olevat
                                  % grafiikkaikkunat
x = [-2:.01:2];                   % Diskreetti x-akselin väli
y = exp(x);                       % exp-funktion arvot x-pisteissä

Lasketaan Taylorin polynomien arvoja x-pisteissä

p0 = ones(size(x));               % Ykkösistä koostuva vektori
p1 = p0 + x;                      % Lisätään x-vektori
p2 = p1 + (x.^2)/2;               % Seuraava termi, huomaa .^
p3 = p2 + (x.^3)/6;               %  Ja vielä yksi.

Piirretään kaikki samaan kuvaan

plot(x,y)
hold on                    % Seuraavat samaan kuvaan pyyhkimättä vanhaa.
plot(x,p0)
plot(x,p1)
plot(x,p2)
plot(x,p3)
grid on
legend('e^x','p0','p1','p2','p3','Location','NW')  % NW=North West
shg  % "Show Graphics" Ei välttämätön, mutta hyvä.

Entä, jos halutaan laskea ja piirtää korkeamman asteisia?

Otetaan käyttöön for-silmukka, tehdään alusta saakka:
p{1}=ones(size(x)); % Matlab ei salli 0-indeksiä. Huomaa aaltosulut {}.
                % p{k} voidaan ajatella indeksoituna (vektori)muuttujana.
N=10;           % Muuttele tarpeen mukaan
for k=1:N
    p{k+1}=p{k}+(x.^(k))/factorial(k);
end
%%
hold off   % Seuraavat piirtokomennot pyyhkivät edelliset piirrokset 
% Valitaan värit:  'r' = red, 'g'=green, 'k'=black 
plot(x,y,'r--',x,p{3},'g',x,p{5},'b',x,p{9},'k') % Parilliset polynomit
hold on    % Pidetään edelliset
plot(x,y,'r--','LineWidth',2)  % Punainen, paksumpi katkoviiva, jotta
                               % erottuisi seuraavassa p{9}:n kuvaajasta.
grid on
title('e^x ja parillisia Taylorin polynomeja')
legend('e^x','T_2','T_4','T_8','Location','NW')
shg
 
 
 
 
% Zoomaamalla gragfiikkaikkunan suurennuslasilla, nähdään, missä vaiheessa
% T_8 :n (eli p{9}:n) ja exp(x):n kuvaajat erkanevat.
%

Approksimointivirhe

figure            % Avataan uusi grafiikkaikkuna.
plot(x,y-p{9})    % Piirretään erotus exp(x)-T8(x)
title('Virhe: e^x - T_8(x)')
grid on
shg

Ja lasketaan maksimivirhe:
maksvirhe=max(abs(y-p{9}))
%
maksvirhe =

    0.0018

Jospa nyt maltan lopettaa tällä kertaa tähän.

Lisää esimerkkejä ja tehtäväehdotuksia seuraavalla kerralla.