Johdatus Mathematicaan

Mat-1.1110 Matematiikan peruskurssi C1

Viikko 41/ Hakula/Pursiainen

Mathematica

Mathematica on laaja ja voimallinen matemaattinen yleisohjelmisto eli numeeris-symbolinen ohjelmisto. Mathematica tukee lukuisia eri ohjelmointitapoja ja sisältää kattavan aliohjelmakirjaston. Kuten aina, laajuus ja voimallisuus merkitsee korkeaa käyttöönottokynnystä. Mathematican suunnittelijoiden valinta poiketa normaalista matemaattisesta notaatiosta tekee erityisesti vasta-alkajan elämän usein hankalaksi. Kuitenkin, jo muutamassa vuodessa ohjelmiston käyttö alkaa tuntua luonnolliselta ja käyttäjä alkaa ohjata ohjelmaa eikä päinvastoin, kuten valitettavan usein tuntuu alussa käyvän.

Mathematicaa voi käyttää joko työarkkirajapinnan tai voimakäyttäjän komentorivin kautta. Tässä rajoitumme vain työarkkeihin. Napauttelemalla tekstin oikealla puolella näyviin "hakasulkeisiin" kappaleet (tai solut) saa sulkeutumaan ja avautumaan.
Näin arkin selausta ylös ja alas voi halutessa vähentää.

Laskin

Kuten arvata saattaa, kaikki kuviteltavissa olevat laskintoiminnot sisältyvät Mathematicaan. Allaolevan laskutoimituksen voi evaluoida eli suorittaa siirtymällä soluun ja painamalla Shift-Return (tai Enter)

1 + 1

Jotain tämänkaltaista olisi pitänyt tapahtua! Olen poistanut mallivastauksesta mahdollisuuden suorittaa sitä uudestaan.

1 + 1

Entäpä sitten nämä?

2 + 3
2^3
2*3
2/3
2 3

Havaitaan, että Mathematica pyrkii säilyttämään tarkimman muodon mahdollisimman pitkälle.Lisäksi todetaan, että kertolasku voi olla myös implikoitu eli kirjoitetun notaation mukainen. Ominaisuutta on syytä käyttää harkiten, vai mitä?

2 3 23 23 2 3

Entä, jos todella halutaan liukulukuja? Tällöin on komento N paikallaan

4/5

N[4/5]

Muuttujat

Arvoja voi tallettaa muuttujiin.

a = 2

a  +  2

Toisaalta,

b := 2

b + 2

Näillä kahdella sijoitusoperaatiolla on merkittävä ero: ensimmäinen suoritetaan välittömästi, mutta jälkimmäinen aina uudestaan, kun muuttujaa käytetään. Seuraava esimerkki lienee valaiseva:

a = 2
b = 2
c = a
d := b
c
d
a = 3
b = 3
c
d

Lopussa c ja d ovat siis eri suuret! Sijoitus a=3 tapahtui c=a:n jälkeen, eikä sillä ole
mitään vaikutusta c:n arvoon. d:n arvoa määritettäessä sen sijaan tarkistetaan aina myös
b:n arvo.

Muutujan arvot poistetaan komennoilla Clear ja Remove.

k = 3

Seuraava sijoitus vastaa komentoa Clear[k].

k =.

k

Epätoivon hetkellä, kun kaiken haluaa alusta, on ydinlataus paikallaan:

Remove["Global`*"]

Kaikki normaalissa istunnossa tehdyt määritelmät on nyt poistettu!

"Aliohjelmat"

Mathematica on täysiverinen ohjelmointiympäristö. Sen tehokas käyttö edellyttää ns. funktionaalisen ohjelmoinnin ymmärrystä, mutta vähemmän tehokkaanakin se on silti tehokas...

Perussilmukka muodostetaan Do-lausekkeella.

Do[ 
      komentoja, 
{indeksi,  alku,  loppu,  
        askel}]

Lasketaan vaikkapa 2+2+2+2+2.

a = 2;  Do[a = a + 2, {i, 1, 4, 1}];  a

Tai

a = 2;  Do[a = a + 2, {4}];  a

Tai

a = 2;  Do[a += 2, {4}];  a

Tapoja on siis monia. Missään vaiheessa tätä ei ole väitetty yksinkertaiseksi!
Välituloksia voi kaiuttaa Print-komennolla.

a = 2;  Do[a = a + 2;  Print[a], {i, 1, 4, 1}];  a

Omia funktioita voi tietenkin määritellä. Mathematica poikkeaa tässä kohtaa monista muista "perinteisemmistä" ohjelmointikielistä ja -ympäristöistä siten, että satunnainen käyttäjä ei useimmiten tule käyttäneeksi aitoja funktioita vaan konstruktioita, jotka erehdyttävästi näyttävät funktioilta.

Muodostetaan ensin polynomi.

x^3 + 2  x + 1

Olisi mukava voida evaluoida polynomin arvoja eri x:n arvoilla. Toimitaan seuraavasti:

poly[x_] := x^3 + 2  x + 1
poly[1]

Perinteisiin ohjelmointikieliin perehtynyt hätäinen päättelijä arvelisi, että on tullut määritellyksi funktio poly, jolla on yksi argumentti, x, jota merkitään x_. Hämmästyttävää kyllä, hätäinen päättelijämme voi käyttää Mathematicaa pitkään, ilman että hän koskaan aavistaa harhautuneensa.

Entäpä paikalliset muuttujat? Toimitaan seuraavasti:

poly2[x_] := 
        Module[{vakio = 1}, 
x^3 + 2  x + 
            vakio
];
poly2[1]

Jottei kukaan tämän tekstin lukija jäisi hätäisten päättelijöiden joukkoon, lienee paikallaan kertoa, että edelliset konstruktiot määrittävät toiminnan evaluoitavalle hahmolle, poly[x_] tai poly2[x_]. Hahmot eivät kuulu tämän johdannon piiriin.
Oikea funktio määritellään seuraavasti:

poly3 = Function[x, x^3 + 2  x + 1];
poly3[1]

Kuvanpiirto

Visualisointiominaisuudet Mathematicassa ovat kohdallaan. Ohessa esimerkkejä yleisistä piirtorutiineista. Uutena piirteenä mainittakoon optiot, eli valinnaiset komennot.
Ne saa selville komennolla Options[komennonnimi]. Valinnainen komento annetaan muodossa optio->arvo.

Plot[Sin[x], {x, -Pi, Pi}];

Options[Plot]

Plot[Sin[x], {x, -Pi, Pi}, Frame -> True, 
        Axes -> None];

Plot[{Sin[x], Sin[2  x], Sin[4  x]}, {x, -Pi, Pi}];

Derivaatta ja integraali

Korkeampi symbolinen laskenta on tekniikan ylioppilaalle aina ilon aihe. Derivointi ja integrointi sujuvat Mathematican avulla leikiten. Peruskomennot ovat D ja Integrate.

f = Sin[2  x] Exp[4  x]

D[f, x]
D[f, {x, 2}]

Ja sitten vielä muutama integraali:

Integrate[x, x]
Integrate[x, {x, 0, 1}]
Integrate[f, {x, 0, 1}]

Mathematicalle ei tuota vaikeuksia käsitellä parametria integraalissa:

Integrate[x Exp[k x], x]

Onko vastaus oikein, kun k=0?

Jos halutaan integroida numeerisesti, voidaan suorittaa numeerinen integrointi.
Integrate korvataan yksinkertaisesti komennolla NIntegrate.

NIntegrate[f, {x, 0, 1}]

Maplesta poiketen:

Timing[NIntegrate[f, {x, 0, 1}]]
Timing[N[Integrate[f, {x, 0, 1}]]]

Eivät todellakaan ole sama asia!

Esimerkki: Monte-Carlo integrointi

Luennolla käsiteltiin Monte-Carlo integroinnin perusidea. Nyt kun meillä on tehokas ja taipuisa laskinlaite käytössämme, voimme tehdä uskoa vahvistavia kokeita.

Olkoon esimerkkinä integraali (∫_0^1) x Sin[x] dx.

Integrate[x Sin[x],  {x, 0, 1}]

ref = NIntegrate[x Sin[x],  {x, 0, 1}]

Satunnainluvut generoidaan komennolla Random. Oletusarvoisesti satunnaisluvut lasketaan välillä [0,1]. 

Muodostetaan ensin taulukko satunnaislukuja:

x = Table[Random[], {10}]

Funktion voi nyt evaluoida suoraan.Huomaa, että taulukon alkioiden laskutoimitukset menevät pisteitäin.

x Sin[x]

Monte-Carlo kvadratuuri vaatii laksettujen arvojen summan. Hätäinen päättelijä sortuisi silmukkaan (tai pahempaan!), mutta luonnollisempaa on operoida koko listaan kerralla.

Apply[Plus, x Sin[x]]
Plus @@ (x Sin[x])

Ja siis ensimmäinen arvio:

Plus @@ (x Sin[x]) / 10

Ei huono!

Parannetaan

Olemme siis menestyksellä saaneet laskutoimitukset suoritettua. Seuraava askel on automatisoida askeleet. Asetetaan parametriksi solmujen (pisteiden) lkm.

MCIntegrate[n_Integer] := 
      Module[{x}, 
x = 
          Table[Random[], {n}]; 
Plus @@ (x Sin[x])/
          n
]

Pikatesti:

MCIntegrate[10]
MCIntegrate[20]
MCIntegrate[30]
MCIntegrate[40]

Parasta on laskea keskiarvot yli testijoukon:

Plus @@ Table[MCIntegrate[10], {10}]/10
Plus @@ Table[MCIntegrate[20], {10}]/10
Plus @@ Table[MCIntegrate[30], {10}]/10
Plus @@ Table[MCIntegrate[40], {10}]/10

Ymmärrystä

Tarkastellaan lopuksi kvadratuurin käyttäytymistä n:n funktiona.

KAMCIntegrate[
        n_Integer] := 
Plus @@ Table[MCIntegrate[n], {10}]/
        10

Talletetaan vektoriin n kokonaislukuja väliltä [10,200] askelpituuden ollessa 10.
Lasketaan integraalien keskiarvot jokaisessa vektorin n pisteessä. Lyhenne f /@ v voidaan lausua muodossa Map[f,v] eli tuloksena saadaan vektori {f[v[[1]], f[v[[2]]],...}. (Mathematican terminologiassa vektori lausutaan lista.)

n = Range[10, 200, 10];
v = KAMCIntegrate /@ n;

Piirretään integraalin arvo jokaisella n.

ListPlot[MapThread[List, {n, v}], 
        PlotJoined -> True];

Talletimme tarkan ratkaisun muuttujaan ref, joten seuraavaksi voimme piirtää virheen itseisarvon kuvaajan.

ListPlot[MapThread[List, {n, Abs[ref - v]}], 
        PlotJoined -> True];

Jotta voisimme varmistua asymptotiikan oikeellisuudesta, tarvitsemme logaritmisen kuvaajan. Tämä löytyy aliohjelmakirjastosta, joka on ladattava erikseen.

<< Graphics`Graphics`

Piirretään kuvaaja, joka ei paljoa sano (vielä!) ja talletetaan kuva muuttujaan g1.

g1 = 
        LogLogListPlot[MapThread[List, {n, Abs[ref - v]}], 
          PlotJoined -> True];

Voimme nimittäin piirtää referenssitrendin.

g2 = 
      LogLogListPlot[MapThread[List, {n, 1/Sqrt[n]}], 
        PlotJoined -> True]

Ja lopulta voimme varmistua siitä, että oma rutiinimme toimii niin kuin pitääkin.

Show[{g1, g2}, PlotRange -> All];

Kommentti

Tämä esimerkki vastaa varsin hyvin normaalia, ammattimaista työskentelyä. Oppimiskynnys on alussa korkea, mutta integroituna yli koko opiskeluajan (ja elämän?) saatava hyöty on ilmeinen.