K. Kumericki 2008-01-21

3. Napredno programiranje u Mathematici

3.1 Liste

Liste i elementi

Liste su vazan dio Mathematica jezika. Srodne su poljima (array) iz standardnih programskih jezika (C, Fortran), ali imaju bitno više svojstava. Elementi liste mogu biti praktički bilo koji Mathematica objekti.

l1 = List[1, 3, 5, 7, 9, 11]
l2 = {Sin, Cos, Log}
l3 = {{}, {"jedan"}, %}

Vitičaste zagrade su samo skraćeni oblik funkcije List. Puni oblik se rijetko koristi.

Elementima liste se pristupa putem funkcije Part, odnosno njenog skraćenog oblika [[ ]]  koji se moze ukucati i putem sekvence Esc-[-[-Esc što onda mozda ljepše izgleda.

Part[l1, 2]
l1[[2]] = "tri"
l1[[2]]

l1[[{3, 
        6}]]
Take[l1, {3, 6}]

l1

Listama se moze manipulirati na skoro svaki zamislivi način. First, Last, Take, Drop, Rest, Extract, Length, Prepend, PrependTo, Append, AppendTo, Insert, Delete, ... (vidi Help).

Zadatak: Izbrišite u gornjoj listi l1 treći element  i dodajte na kraj još dva elementa, brojeve 13 i 15.

Elementi liste mogu biti i druge liste, pa onda imamo višedimenzionalne liste. Dvodimenzionalne liste (listu lista jednake duljine) smo već sreli kao matrice. Evo primjera trodimenzionalne liste:

l4 = {{{a, b}, {h}, {d, e}}, {{f}, {g, h}, {i, j, k}}}

l4[[1]]
l4[[1]][[3]]
l4[[1, 
      3]]
l4[[1, 3, 2]]

Prilikom programiranja u Mathematici se razmjerno često javlja potreba za smanjivanjem dimenzije liste tj. za micanje zagrada. Funkcija Flatten bez dodatnog argumenta pretvara svaku listu u jednodimenzionalnu tj. miče sve unutarnje zagrade, dok dodatnim argumentom mozemo zatraziti da to ide samo do neke razine.

Flatten[l4]
Flatten[l4, 1]

Konstruiranje lista

Range[10]

Range[2, 5, 0.5]

Range radi samo s brojevima. Univerzalnija funkcija za konstruiranje lista je Table.

Table[((x^n)/(n!)), {n, 1, 5}]

tbl = Table[a[i, j], {i, 1, 3}, {j, 1, 2}]

tbl  // TableForm

Nekoliko Mathematica funkcija (Integrate, Sum, ...) poput Table kao argumente imaju i objekte koji se zovu iteratori i dolaze u slijedećim oblicima:

({{{n, nmin, nmax, step}, izraz se ponavlja,a n ide od nmin do nmax, uz korak step}, {{n, nmin, nmax}, izraz se ponavlja,a n ide od nmin do nmax, uz korak 1}, {{n, nmax}, izraz se ponavlja,a n ide od 1  do nmax, uz korak 1}, {{n}, izraz se ponavlja n puta}})

3.2 Funkcije

Mogućnost definiranja novih funkcija je osnovna stepenica k naprednijem programiranju. Nove funkcije u Mathematici  definiramo na slijedeći način

f[x_] := (x^2)

f[2]

f[(((2  x + 3))^2)]

Vazno je uociti da je prilikom definicije funkcija potrebno upotrebljavati tzv.  odgodjeno pridruzivanje ":=" sto za posljeedicu ima da se definicija funkcije izvrijednjava tek kad je pozovemo s nekim konkretnim argumentom. Nadalje, primijetite simbol "_"  (podvlaka) na lijevoj strani, tzv.  blank. Izrazi oblika "x_"  se zovu  uzorci (patterns).

({{x, jedinstveni simbol x}, {x_, \(bilo\ koji\ izraz, \ 
            privremeno\ \((za\ vrijeme\ izvrijednjavanja\ funkcije)\)\ \ zvan\ x\)}})

Naravno, funkcija moze biti i od više varijabli:

f[x_, a_] := (x^a)

f[4, 2]

Simplify[%]

Moguće je i ograničiti domenu neke funkcije pomoću funkcije Condition, koja se obično zapisuje u skraćenom infiks obliku "/;" (čitaj npr.  "uz uvjet"):

p[x_ /; x > 0] := (x^3)

"p(x), uz uvjet da je x veći od nula, je (x^3)."

p[3]

p[-3]

Tako mozemo definirati funkciju po dijelovima

step[x_ /; x ≥ 0] := 1
step[x_] := 0    /;  x < 0

Za negativni dio domene smo, demonstracije radi, upotrijebili alternativnu poziciju operatora "/;". U starijim verzijama Mathematice  samo je taj način bio dozvoljen.

Plot[step[x], {x, -1, 1}, 
      PlotStyle -> Thickness[0.01]]

Zadatak: Definirajte funkciju  plotpoint[z_] koja crta kompleksnu ravninu i na njoj točku koja odgovara kompleksnom broju z. (Treba vam funkcija ListPlot, te funkcije Re i Im za uzimanje realnog i imaginarnog dijela.)

Domaći zadatak 3-1 (relativistička kinetička energija)

Formule za nerelativističku i relativističku kinetičku energiju tijela mase m koje se giba brzinom v su

(Subscript[K, nr]) = (1/2) m (v^2)                         (Subscript[K, rel])=m (c^2)[(1/(Sqrt[1 - ((v^2)/(c^2))])) - 1]

(a) Definirajte odgovarajuće funkcije knr i krel te usporedite na istom dijagramu ponašanje tih funkcija za brzine od 0 do 3×(10^8) m/s za neku vrijednost mase. Označite koordinatne osi!
(b) Odredite brzinu pri kojoj je greška nerelativističke formule  tisućinku promila.

3.3 Transformacijska pravila

Transformacijska pravila smo već sreli prilikom provjere rješenja jednadzbi, ali njihova je primjena univerzalna.

izr = a + b Log[(x^2)] + c x Sin[x];

izr  /.  {p_^2 -> "nesto na kvadrat", x -> y,  
        c -> x}

Primijetite da se pravila primjenjuju po redu kako su navedena i da u njima smijemo koristiti uzorke.

Ako transformacijsko pravilo involvira izvrijednjavanje neke funkcije  najcesce je potrebno umjesto strelice "->" koristiti  tzv. odgodjeno transformacijsko pravilo  ":>" (Mathematica ga tipografski u biljeznici prikazuje kao ":>") koje se izvrijednjuje tek prilikom upotrebe pravila:

2 + 3  I  /. z_ :> plotpoint[z]

{3 + I, 7 - 3  I} /. z_ :> plotpoint[z]

Dogodilo se to da se uzorak  z_ "spario" s cijelom listom {3 + I, 5 + I, 7 - 3  I}, a funkcija plotpoint ima smisla samo ako joj je argument broj.

Potrebno je ograniciti sparivanje uzoraka. Postavit cemo uvjet da je z broj:

{3 + I, 
        7 - 3  I} /. (z_ /; NumericQ[z]):>plotpoint[
          z]

NumericQ je specijalna funkcija koja testira da li je neka velicina numericka ili nije.

NumericQ[d]

NumericQ[2/3]

Postoje i mnoge druge funkcije koje kao rezultat vracaju True ili False pa se mogu koristiti za specificiranje uvjeta djelovanja transformacijskih pravila:

?*Q

Uvjete se moze kombinirati funkcijama And, Or, Not, ...

{2, 3, 6, 
        10} /. (x_ /; And[EvenQ[x], x < 8]) -> "pogodak"

Zadatak: Kreirajte listu brojeva 1,2,3,....,100 i u njoj sve slozene brojeve (dakle koji nisu prim-brojevi, vidi PrimeQ), zamijenite nulom.

3.4 Proceduralno programiranje

Osnovna paradigma proceduralnog programiranja je da se kreće od nekog definiranog stanja programa (vrijednosti varijabli i nizova varijabli u memoriji računala). To stanje  se specificira početnim naredbama pridruzivanja (x=0, i=1)  (assignment) i onda se algoritam izvodi primjenom raznih grananja, petlji i potprograma koji svi mijenjaju to stanje i vode na konačno stanje u kojem vrijednosti nekih varijabli odgovaraju rješenju problema.

If-then grananje

Grananje  se u večini proceduralnih jezika izvodi nekim oblikom if-then-else petlje. Tu mogućnost nudi i Mathematica putem funkcije If

If[test, 
        <komande koje se izvršavaju ako je test True ("then")>,
         <komande koje se izvršavanju ako je test False ("else")>]

Naglasimo da je If normalna Mathematica funkcija. Rezultat njenog izvrijednjavanja je rezultat zadnje u nizu komandi sadržanih u argumentu br. 2 ili 3 (koji od tih nizova će se izvršiti ovisi o testu).

Na primjer, funkcija step koju smo sreli u prošlom poglavlju može se realizirati i ovako

step[x_] := If[x < 0, 0, 1]

Plot[step[x], {x, -1, 1}, 
      PlotStyle -> Thickness[0.01]]

Ako želimo konstrukciju s višestrukim grananjem if-elif-elif-else, onda moramo koristiti "ugniježđeni"  (nested) If:

f[x_] := If[x < -1 , (x^2), If[x > 1, (x^2), x]]

Plot[f[x], {x, -3, 3}, 
      PlotStyle -> Thickness[0.01]]

Spomenimo usput da ovako definiranu funkciju slobodno deriviramo.

Plot[f'[x], {x, -3, 3}, 
      PlotStyle -> Thickness[0.01]]

Primijetite da smo gore mogli izbjeći ugnjezdjeni If korištenjem logičkih funkcija, danih u slijedećoj tablici (Četvrti stupac navodi niz tipaka koje daju znak u trećem stupcu. Ako zagrade ne kazu drugačije, prvenstvo izvršavanja je odozgo prema dolje.):

({{Not, !, ¬, Esc - not - Esc}, {And, &&, ∧, Esc - and - Esc}, {Or, ||, ∨, Esc - or - Esc}})

f[x_] := If[x < -1  || x > 1, (x^2), x]

Or[1 > 3, -1 > 0, ¬ 5 > 10]

Petlje

Sve standardne vrste petlji, do-petlja, while-petlja i for-petlja, postoje i u Mathematici. Sintakse su

Do[komande, <iterator>]

gdje je <iterator> objašnjen u prošlom poglavlju.

While[test, komande]

gdje se komande (niz Mathematica komandi odvojenih točka-zarezom ";") izvršavaju sve dok je god test  True.

For[start, test, incr, komande]

gdje se prvo izvrši start, a zatim se ponavlja izvršavanje komandi i incr  sve dok je god test True (vrlo slično sintaksi for-petlje u C programskom jeziku).

Primjer upotrebe while-petlje smo vidjeli prilikom definiranja funkcije faktorijelProc u uvodu ovog poglavlja. Evo kratkih primjera za do- i for-petlje:

Do[(2^n), {n, 0, 1, 0.5}]

Nema rezultata! Sve ove petlje, premda po sintaksi normalne Mathematica funkcije, imaju svojstvo po kojem se razlikuju od većine ugrađenih funkcija, a to je da ne vraćaju nikakav rezultat. Preciznije, vraćaju izraz Null. (Isti izraz vraćaju i one komande nakon kojih stavimo točka-zarez.)

FullForm[%]

Ukoliko želimo neko ispisivanje rezultata koraka petlje, koristimo Print

Do[Print[(2^n)], {n, 0, 1, 0.5}]

Za iskakanje iz petlje možemo koristiti Break ili Return. Npr. slijedeći algoritam pronalazi prvi cijeli broj čiji rastav sadrži više od pet različitih prostih faktora.

For[
      i = 1, 
True, 
i = 
        i + 1, 
If[Length[FactorInteger[i]] > 5, 
        Break[], ]]; i

(Primijetite da If nema treći argument tj. da je on Null.) Ovo je relativno neelegantan program za pronalaženje najmanjeg broja koji ima više od pet različitih faktora. Bolje bi bilo

i = 1; While[Length[FactorInteger[i]] < 6,  i = i + 1]; i

a najbolje

Times   @@   Table[Prime[n], {n, 6}]

ali ovo zadnje spada već u funkcionalno programiranje.

Moduli

Jedan od problema s proceduralnim programiranjem je obilje popratnih efekata koje procedure izazivaju, poput pridjeljivanja vrijednosti pomoćnim varijablama. Na primjer, gornji račun je pridijelio vrijednost varijabli i

i

što je zapravo nepotrebno i neželjeno, a često dovodi i do grešaka. Funkcionalno programiranje takve stvari obično automatski izbjegava (vidi zadnji primjer), dok se kod proceduralnog programiranja to može izbjeći korištenjem modula tj. funkcije Module

Module[{x,  y = (Subscript[y, 0]),  …}, komande]

koja specificira da su varijable x, y, … lokalne tj. da ne postoje izvan modula, te koja vraća kao rezultat rezultat posljednje iz niza komandi. Primijetite da se prilikom specificranja lokalnih varijabli one mogu i inicijalizirati na neku vrijednost. Primjenu modula smo vidjeli u uvodu, pri definiranju funkcije faktorijelProc.

a = 1

f[n_] := Module[{a}, a = (((x + y))^n); Expand[a]]

f[3]

a

Varijabla a je zadržala svoju globalnu vrijednost bez obzira na to što joj je unutar modula pridjeljivana neka druga vrijednost.

Clear[a]

Zadatak: Isprogramirajte proceduralno funkciju fibProc[n_]koja računa n-ti Fibonaccijev broj.

Domaći zadatak 3-2: Isprogramirajte  funkciju kredit[iznos_, kamata_, rok_]koja daje mjesečni anuitet i cijenu kredita za zadani iznos, kamatu i rok kredita u godinama. Cijena kredita je razlika između dobivenog i konačno otplačenog novca.

Nadalje, isprogramirajte funkciju kreditTablica[iznos_, kamata_, {start_,stop_,step_}]koja ispisuje tablicu za niz rokova kredita zadanih putem naznačenog iteratora. Stupci tablice neka budu: rok, anuitet, cijena.

Zadatak: Konstruirajte funkciju brown[n_] koja daje listu {{(Subscript[x, 0]),  (Subscript[y, 0])},{(Subscript[x, 1]), (Subscript[y, 1])},…,{(Subscript[x, n]),(Subscript[y, n])}} pozicija čestice koja se giba u ravnini slučajnim gibanjem koje je definirano rekurzijom (Subscript[x, i])=(Subscript[x, i - 1])+rnd1  i  (Subscript[y, i])=(Subscript[y, i - 1])+rnd2  gdje su rnd1 i rnd2 slučajni brojevi između -1 i 1. Nacrtajte odgovarajuću putanju čestice.

Naputak:
1. Programirajte proceduralno koristeći Module
2. Definirajte varijablu randomVec koja će pri svakoj upotrebi (pozivu) poprimati vrijednost novog slučajnog vektora (x, y). Za to će biti potrebno koristiti odgođeno pridruživanje vrijednosti toj varijabli.
3. Inicijalizirajte listu pozicija tako da sadrži (0,0) kao prvi vektor pozicije i putem petlje dodajte u svakom koraku toj listi novi slučajni vektor.
4. Za crtanje možete koristiti ListPlot

3.5 Pakiranje koda u pakete

Mathematica paketi su nacin da se racunalni kod koji obavlja neku specijaliziranu funkciju organizira u jednu datoteku koja se ucitava po potrebi. Npr. u biljeznici 2-MatematikaS.nb se koristi paket za vektorsku analizu koji se ucitava na slijedeci nacin:

<< Calculus`VectorAnalysis`

Nakon ucitavanja mozemo koristiti njegove funkcije, npr. za izracunavanje gradijenta, divergencije ili rotacije polja. Npr. polje sile koje odgovara centralnom potencijalu V(r) = -(1/r):

SetCoordinates[Cartesian[x, y, z]];
-Grad[-(1/(Sqrt[(x^2) + (y^2) + (z^2)]))]

Dobili smo vektor kojeg prepoznajemo kao (Overscript[F, →])(r)=- ((Overscript[r, →])/(r^3))

Zadatak: Pronadjite unutar instalacije Mathematice datoteku VectorAnalysis.m, u kojoj je implementiran gornji paket, i promotrite njenu strukturu.

Uobicajena struktura datoteke "ime paketa.m" koja implementira neki paket je:

(* Opci podaci (autor, verzija, ...) u obliku Mathematica komentara *)

BeginPackage[" ime paketa`"];

(* Popis "eksportiranih" funkcija dostupnih korisniku i njihov opis *)

funkcijaF::usage = "funkcijaF[var1,var2,...] radi to i to ..."
funkcijaG::usage = "funkcijaG[var1,var2,...] racuna to i to ..."

Begin["`Private`"];

(* definicija raznih lokalnih pomocnih varijabli i funkcija te definicija eksportiranih funkcija *)

End[];

EndPackage[];

Kljucno je uociti da ce korisniku biti izravno dostupne samo funkcije navedene u "::usage" linijama prije Begin["`Private`"]. Te "::usage" linije sluze i za dokumentiranje koda. Kad korisnik ukuca ?funkcijaF dobit ce kao odgovor desnu stranu "::usage" linije. Imena eksportiranih funkcija treba odabrati tako da nije vjerojatno da iste funkcije postoje u nekim drugim paketima, ali pomocne funkcije koje ne eksportiramo mogu imati proizvoljna imena.

Ukoliko paket za svoje funkcioniranje treba neke druge pakete (npr. VectorAnalysis), ispravan nacin njihovog pozivanja je ukljucivanje u komandu BeginPackage

BeginPackage[" ime paketa`", "Calculus`VectorAnalysis`"]

Zadatak: Isprogramirajte funkciju sila[potencijal_] koja na gornji nacin racuna silu iz zadanog potencijala u kartezijevom sustavu i to prvo kao Module, a zatim kao paket u datoteci sila.m.

Zadatak: Modificirajte gornji paket tako da radi striktno u x-y ravnini (eliminirajte z-koordinatu). Testirajte s potencijalima (Subscript[V, 1])(r)=- (1/r) = -(1/(Sqrt[(x^2) + (y^2)])) i (Subscript[V, 2])(r)=sin(x y). Dodajte paketu funkciju poljeSile[potencijal_, {(Subscript[x, min]),(Subscript[x, max])}, {(Subscript[y, min]),  (Subscript[y, max])}] koja crta vektorsko polje sile koristenjem paketa Graphics`PlotField`

Zavrsni zadatak: Isto kao u proslom zadatku, dakle prvo kao Module, a zatim kao paket u datoteci gibanje2D.m, isprogramirajte funkciju gibanje2D[potencijal_,{initr_,initv_},iter_] koja crta putanju gibanja cestice jedinicne mase u proizvoljnom 2D potencijalu za proizvoljne pocetne uvjete (Subscript[(Overscript[r, →]), 0]) = {(Subscript[x, 0]), (Subscript[y, 0])}=initr i (Subscript[(Overscript[v, →]), 0]) = {(Subscript[vx, 0]), 
            (Subscript[vy, 0])}=initv, za vrijeme zadano u iteratoru iter={t,t0,tmax}. Predati (emailom) trebate paket gibanje2D.m i demonstracijski notebook <Vase prezime>.nb u kojem cete opisati sto Vas paket radi, imajuci u vidu obicnog citatelja (npr. nekog studenta fizike), a ne mene!  (15 bodova)

Naputak: Postupajte kao u zadatku s kraja biljeznice 2-MatematikaS.nb. Iscrtavanje mozete napraviti u nekom fiksiranom dijelu x-y ravnine. Rezultat za centralni potencijal treba biti npr.

gibanje2D[-(1/(Sqrt[(x^2) + (y^2)])), {{-3, -3}, {0.6, 
          0.7}}, {t, 0, 7}]

Neobavezni bonus: (a) Dodajte svojstvo da program crta samo onaj dio ravnine u kojem se cestica giba (Dakle, morate prije iscrtavanja odrediti maximalne i minimalne vrijednosti x i y koordinata putanje za zadani vremenski interval {t, (Subscript[t, 0]),(Subscript[t, max])}). (5 bodova)
(b) Napravite funkciju koja proizvodi animaciju gibanja koristeci  paket Graphics`Animation`. (3 boda)