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. Suositeltavampaa 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
|
Klikkaa kuva tarkemmaksi/suuremmaksi, jos siltä tuntuu. Samoin seuraavat kuvat.
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:
-
Syksyllä 2016 pitämäni kurssin luentokalvoja (englanniksi)
- Erään Matlab-oppaan ("YAGTOM") suomenkielinen mukaelma funktioaiheesta
- 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ä?.
