2008-11-12

4. Programiranje u Mathematici

Uvod

Isti matematički problem je često moguće riješiti na različite načine, putem algoritama iza kojih stoje sasvim različiti načini razmišljanja. Razložimo to na primjeru funkcije faktorijel.

7!

Klasični način programiranja, upotrebljavan od dana prvih računala, je tzv. proceduralno (ili imperativno) programiranje kod kojeg  na program gledamo kao na niz naredbi koje obično postupno mijenjaju vrijednosti nekih varijabli pohranjenih u memoriji računala:

faktorijelProc[n_]:= Module[{i=1,fac=1},
While[i<n,
i=i+1;fac=fac*i];
fac];

faktorijelProc[7]

Zanemarite zasad precizno značenje komande Module koju ćemo objasniti kasnije, te preciznu sintaksu komande While. Ovdje je očito riječ o standardnom algoritmu kakvog bismo isprogramirali u bilo kojem proceduralnom jeziku poput Fortrana ili C-a: Imamo petlju koja se prolazi n puta i svaki puta množimo rezultat prošlog prolaza sa za 1 većim brojem.

Međutim, Mathematica nam omogućuje i drugačije pristupe. Tako je  kod funkcionalnog programiranja naglasak na izvrijednjavanju funkcija tj. primjeni funkcija na izraze. Faktorijel prirodno zamišljamo kao funkciju koja daje "umnožak svih brojeva od 1 do zadanog broja":

faktorijelFunkc[n_]:=Apply[Times, Range[n]]

faktorijelFunkc[7]

(Primijetite da nam ovdje nije trebala pomoćna varijabla (poput i iz faktorijelProc), ali nam je trebalo više memorije jer smo trebali pohraniti čitavu listu {1, 2, …, n} koju daje Range[n].

Programiranje transformacijskim pravilima u Mathematici  je svojevrsna nadgradnja funkcionalnog, gdje funkcije promatramo kao skupove pravila za transformaciju svojih argumenata. Uz ime funkcije je vezan niz pravila koja čine da se funkcija ponaša kao složeno transformacijsko pravilo. Tu smo filozofiju u prošlom poglavlju primijenili za definiciju funkcije step. Faktorijel možemo unutar te paradigme definirati  kao skup od dva pravila:  jedno koje transformira n-> n (n-1)! i koje Mathematica primjenjuje rekurzivno, te jedno trivijalno koje transformira 1! -> 1 i tako "zaustavlja" rekurziju

faktorijelTrans[1]:=1
faktorijelTrans[n_]:=n faktorijelTrans[n-1]

faktorijelTrans[7]

Primjetite da se ovdje koristi svojstvo da Mathematica uvijek prvo izvrijednjava specifična pravila pa tek onda prelazi na opća.

U ovom konkretnom slučaju, funkcionalni program se čini najelegantniji od sva tri. No, promotrimo funkciju koja izračunava n-ti član Fibonaccijevog niza (niz kod kojeg  je  svaki član definiran kao zbroj prethodna dva, a javlja se pri analizi idealizirane populacije zečeva). Ovdje se kao najelegantniji pokazuje pristup putem transformacijskih pravila:

fibTrans[0]:=0
fibTrans[1]:=1
fibTrans[n_]:=fibTrans[n-1]+fibTrans[n-2]

fibTrans[7]

U posljednje vrijeme je vrlo popularno i tzv. objektno  (object oriented) programiranje. Mathematica nije sasvim prilagođena istom  pa ga nećemo razmatrati. (Na internetu se može naći dodatni paket Class.m koji omogućuje neku vrstu objektnog programiranja u Mathematici.)

4.1 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 pridruživanja (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 ugnježđeni 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 kažu 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 specificiranja 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 4-1: Isprogramirajte proceduralno funkciju fibProc[n_]koja računa n-ti Fibonaccijev broj.

Zadatak 4-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 4-3: 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

Zadatak 4-4: Konstruirajte funkciju walk1D tako da simulira tzv. jednodimenzionalnog nasumičnog šetača. Svaki šetačev korak treba biti iste, jedinične duljine, ulijevo ili udesno. Dakle, umjesto vektora (rnd1, rnd2) iz gornjeg Brownovog gibanja trebamo slučajan odabir koraka +1 ili -1. Provjerite ispravnost tvrdnje da je očekivana udaljenost nasumičnog jednodimenzionalnog šetača od ishodišta nakon n koraka (Sqrt[n]).

4.2 Funkcionalno programiranje

Funkcionalno programiranje je stil programiranja koji stavlja naglasak na izvrijednjavanje izraza,  a ne na sukcesivno izvršavanje komandi. Kod čistih funkcionalnih programa obično nema pomoćnih varijabli i nepotrebnih popratnih efekata izvrijednjavanja funkcija. U tom smislu funkcionalno programiranje je pomalo slično radu s tabličnim kalkulatorom (npr. MS Excel): redosljed izračunavanja čelija je nebitan tj. očekujemo automatsku konzistenciju. Isto, vrijednosti čelija su dane izrazima, a ne slijedovima komandi. (Više informacija o funkcionalnom programiranju nudi  Functional Programming FAQ   ili   WWW stranice  Haskell funkcionalnog jezika.)

Takav stil programiranja često rabi nekoliko specijalnih funkcija (mogli bismo ih zvati "meta-funkcije") čija je uloga upravljanje primjenjivanjem drugih funkcija koje dolaze kao argumenti ovih meta-funkcija. Upoznajmo ih.

Map i Apply

Map[funkcija, izraz] primjenjuje funkciju na sve elemente izraza. (Preciznije, na sve pod-izraze prvog nivoa izraza — vidi odjeljak "Sve su izrazi" iz prošlog poglavlja za objašnjenje nivoa izraza, ali trebalo bi biti jasno iz slijedećih primjera.) Tako ako je izraz lista, onda se funkcija primijenjuje na sve elemente liste.

Map[fun,{a,b,{c,d}}]

Zgodna je i infiksna forma ove funkcije koja se dobiva simbolom /@

fun /@ (a (x^2)+b)

Druga važna "meta-funkcija" je Apply[funkcija, izraz] koja zamjenjuje glavu izraza s funkcijom. Ona isto ima često korištenu infiks formu @@.

Apply[fun,{a,b,{c,d}}]

fun @@ izr

Dakle, ako izraze promatramo kao u prošlom poglavlju tj. kao generičku formu  (Subscript[i, 0])[(Subscript[i, 1]), (Subscript[i, 2]),…] onda Map djeluje na argumente (Subscript[i, 1]), (Subscript[i, 2]),…, a Apply na glavu izraza (Subscript[i, 0]).

Naravno, kao argumente ovih meta-funkcija možemo koristiti i funkcije koje smo sami definirali

faktorijelFunkc[4]

faktorijelFunkc /@ Range[6]

Recimo da sad želimo ispisati tablicu s nizom prirodnih brojeva i njihovih faktorijela. Možemo prvo postupiti tako da prvo definiramo funkciju koja generira jedan red tablice i onda je pomoću Apply primijenimo na niz brojeva:

faklist[n_]:={n,n!}

faklist /@ Range[6] //TableForm

Međutim, baš kao i pomoćne varijable, tako ni pomoćne funkcije nisu u duhu funkcionalnog programiranja. Zato postoje tzv. čiste funkcije ("pure functions") koje su bezimene i definiraju se i upotrebljavaju na slijedeći način:

{#,#!}&[4]

Dakle, "&" nakon izraza je samo oznaka da je ono što prethodi funkcionalni oblik funkcije, a "#" stoje na mjestu argumenata.

Tako gornju tablicu dobijemo ovako:

{#,#!}& /@ Range[6] //TableForm

Ako čista funkcija treba biti funkcija više argumenata, koristimo #1, #2, itd. Dakle, drugim riječima, ako nam negdje treba pomoćna funkcija koju bismo definirali kao:
aux[x1, x2, …]:= (izraz koji uključuje x1, x2, …)
Možemo na mjesta na kojem bismo je upotrijebili staviti čistu funkciju:
(izraz koji uključuje #1, #2, …)&

Nest i Fold

Ove dvije funkcije  služe repetitivnom primjenjivanju neke funkcije na dani argument.

Nest[funkcija, argument, n]  primjenjuje funkciju na argument pa onda na dobiveni rezultat opet primjenjuje funkciju i tako n puta ("ugnježđivanje", nesting).

Nest[Log,x,3]

(((Sqrt[5])+1)/2)//N

Npr. tzv. zlatni omjer  (((Sqrt[5])+1)/2)  se može izraziti i kao beskonačni razlomak

GoldenRatio = 1+(1/(1+(1/(1+(1/(1+(1/("…"))))))))=1.618034...

Takav razlomak možemo izvrijedniti na slijedeći način:

Nest[(1+(1/(#)))&,x,5]

gdje smo upotrijebili čistu funkciju koja odgovara funkciji f(x)=1+(1/x).

% /. x->1 //N

To isto, ali elegantnije:

Nest[(1+(1/(#)))&,1,20] //N

Funkcija NestList[funkcija, argument, n]  nam daje listu s rezultatima primjene funkcije na argument  i to od nula do n puta (dakle kao da daje sve međukorake primjene funkcije  Nest):

NestList[fun,x,3]

NestList[(1+(1/(#)))&,1,8] //N

Funkcija FixedPoint je slična Nest — ona isto "ugnježđuje" funkciju, ali ne neki zadani broj puta već sve dok nema promjene rezultata (do na određenu točnost koja se može specificirati).

FixedPoint[(1+(1/(#)))&,1] //N

Ovo se vrti beskonačno! Probajmo ovako:

FixedPoint[(1+(1/(#)))&//N,1]

Zadatak 4-5: U čemu je bio problem tj. zašto se prva verzija vrtila beskonačno?

(Inaće, test koji FixedPoint provodi da bi provjerio da nema promjene rezultata (po defaultu je to ugrađeni predikat SameQ) se može i posebno specificirati putem opcije SameTest. Vidi donji program mandelbrot.)

Funkcija Fold[fun, (Subscript[a, 0]), {(Subscript[a, 1]),(Subscript[a, 2]), ⋯, (Subscript[a, n])}], gdje je fun funkcija od dva argumenta, vraća fun[[⋯fun[fun[(Subscript[a, 0]), (Subscript[a, 1])], (Subscript[a, 2])]⋯]].

Fold[fun,a,{b,c,d}]

Npr. faktorijel možemo implementirati i ovako

faktorijelFunkc2[n_]:=Fold[Times,1,Range[n]]

faktorijelFunkc2[7]

Ovdje imamo isto i funkciju FoldList koja nam daje listu međukoraka

FoldList[Times,1,Range[7]]

Ove se funkcije na prvi pogled čine egzotične, ali često se mogu korisno primijeniti u funkcionalnom programiranju. Npr. funkcija brown koja daje listu pozicija čestice koja se giba nasumičnim Brownovim gibanjem se može implementirati i ovako

randomVec:=RandomReal[{-1,1},2]
brown[n_]:=FoldList[#1+#2&,{0,0},Table[randomVec,{n}]]

ListPlot[brown[300],Joined->True,AxesOrigin->{0,0},AspectRatio->1]

Primijetite da je program elegantniji od proceduralnog. Nema pomoćnih varijabli, izuzev randomVec kojeg bismo se isto mogli riješiti tako da njegovu definiciju stavimo u definiciju funkcije brown i tako dobijemo one-liner. S druge strane, program je nešto neoptimalniji od proceduralnog jer koristi više memorije—trebaju mu dvije liste veličine n: jedna da drži n vektora pomaka randomVec, a druga da drži pozicije čestice. Proceduralni program smo izveli samo s jednom listom.

Rekurzivno definiranje funkcije (*)

Jedna od značajki funkcionalnog programiranja je i upotreba rekurzivne definicije funkcija. Već smo je sreli u prošlom poglavlju pri definiciji faktorijela i Fibonaccijevih brojeva. Evo još jednog simpatičnog primjera, enormno brzorastuće Ackermannove funkcije (važne u matematičkoj teoriji složenosti) definirane putem dvostruke rekurzije:

ackermann[0,n_]:=n+1
ackermann[m_,0]:=ackermann[m-1,1]
ackermann[m_,n_]:=ackermann[m-1,ackermann[m,n-1]]
$RecursionLimit=(10^4);
$IterationLimit=(10^8);

{ackermann[2,0],ackermann[2,1],ackermann[2,2]}

{ackermann[3,0],ackermann[3,1],ackermann[3,2]}

Pažnja, slijedeći red zamrzne moj notebook (da se to ne bi dogodilo, slijedeća čelija nije izvršna: u izborniku Cell->Cell Properties->Cell Evaluatable je isključeno).

{ackermann[4,0],ackermann[4,1]}

ackermann(4,1)=65533, a ackermann(4,2) je 19729-znamenkasti broj (2^6536)-3    (http://www.kosara.net/thoughts/ackermann42.html). Ti brojevi još nisu sami po sebi preogromni za Mathematicu, ali potreban broj rekurzija da se do njih dođe je prevelik.  S druge strane ackermann(5,2) se ne da zapisati u decimalnom zapisu, čak i u računalu napravljenom od svih čestica našeg svemira.

Primjeri razvoja funkcionalnih programa

Promotrimo razvoj programa koji ispisuje tablicu frekvencija pojavljivanja elemenata u listi. Načinimo prvo jednostavnu testnu listu

lst={a,b,a,c,a,d,b,c,c};

Funkcija koja broji broj pojavljivanja je Count[lista, element]

Count[lst,c]

Da bismo ispisivali tablicu trebamo listu oblika {{element1, frekvencija elementa1}, {element2, frekvencija elementa2},...}. Element te liste (red tablice) se može dobiti ovako:

el[x_]:={x,Count[lst,x]};
el[c]

Tablicu frekvencija za različite elemente dobijemo korištenjem čiste funkcije analogne el i njenom primjenom pomoću funkcije Map na popis elemenata

Map[{#,Count[lst,#]}&,{a,b,c,d}] //TableForm

To je sad to, jedino još želimo izbjeći da sami moramo ručno identificirati koji se sve elementi pojavljuju u listi. To nam može raditi funkcija Union.

Union[lst]

Dakle, imamo

Map[{#,Count[lst,#]}&,Union[lst]]

Kao zadnju stvar, poželjno je listu sortirati po frekvencijama. Funkcija Sort  je već definirana tako da zna sortirati liste različitih objekata, no ovdje bi po defaultu sortirala parove po prvom elementu i to po abecednom redu, a nama treba sortiranje po drugom elementu i to po veličini. Stoga moramo kao drugi argument funkciji Sort dati sortirajuću funkciju koju ćemo sami definirati. Sortirajuća funkcija je  predikat koji prima dva argumenta i vraća True ili False ovisno o tome da li prvi argument treba biti sortiran ispred ili iza drugog.

Razmislimo kako definirati funkciju koja kao argument prima dva para i vraća True ili False ovisno o tome da li je drugi element prvog para veći od drugog elementa drugog para ili ne. Može ovako:

frekventnijiQ[par1_,par2_]:= par1[[2]]>par2[[2]]

Sort[Map[{#,Count[lst,#]}&,Union[lst]] ,frekventnijiQ]

Naravno, elegantnije je to izvesti pomoću čiste funkcije:

Sort[Map[{#,Count[lst,#]}&,Union[lst]] ,#1[[2]]>#2[[2]]&]

I to je to.  Sad definiramo kompletnu funkciju:

frekvencije[l_]:=Sort[Map[{#,Count[l,#]}&,Union[l]] ,#1[[2]]>#2[[2]]&]

frekvencije[lst]

Primijenimo to na frekvenciju pojavljivanja slova u Shakespeareovom Mletačkom trgovcu (zapravo koristimo engleski original The Merchant of Venice, skinut sa WWW stranica projekta Gutenberg)

venice=ReadList["files/venice.txt",Character];

Length[venice]     (* ukupni broj slova *)

Take[frekvencije[venice],10]//TableForm

Vidimo da je "e" daleko najčešće slovo (poslije razmaka), a ta je činjenica vrlo važna pri dešifriranju engleskih tekstova.

Kao slijedeći primjer, isprogramirat ćemo funkciju mandelbrot koja crta fraktal poznat pod imenom Mandelbrotov skup. Ovaj skup je definiran kao skup svih točaka c kompleksne ravnine za koje iteracija

(Subscript[z, 0])=0
 (Subscript[z, n+1])=((Subscript[z, n])^2)+c

ne divergira.

Promotrimo najprije prvih n iteracija za neki proizvoljni kompleksni broj, koje možemo dobiti funkcijom NestList kojom iterativno n puta primjenimo funkciju f(z)=(z^2)+c na (Subscript[z, 0])=c kao početnu vrijednost (čime samo štedimo jedan korak u odnosu na (Subscript[z, 0])=0 kao početnu vrijednost)

iter[c_,n_]:=NestList[((#)^2)+c&,c,n]

iter[1+I,5]

Pojavljivanje divergencije možemo uočiti gledajući apsolutne vrijednosti (udaljenost od ishodišta kompleksne ravnine)

iter[c_,n_]:=Abs /@ NestList[((#)^2)+c&,c,n] //N

iter[1+I,5]

Uočite gore primjenu funkcije Map za distribuciju funkcije Abs preko liste. Alternativno se može iskoristiti činjenica da je Abs listabilna funkcija i samo definirati:

iter[c_,n_]:=Abs[NestList[((#)^2)+c&,c,n]] //N

U svakom slučaju, očito je da ovo divergira i da 1+i ne pripada Mandelbrotovom skupu. Pokušajmo s nekim drugim brojem:

iter[-1+0.2I,50]

Izgleda da -1+0.2i, ne divergira tj. da pripada skupu.

Kao nekakav prvi test divergencije možemo uzeti da ako apsolutna vrijednost ne pređe 2 nakon 50 koraka, vjerojatno nikad ni neće. (Matematičar Fatou je metodama kompleksne analize pokazao da kad jednom prođemo 2, divergiranje je neminovno.) Za jedan mali broj točaka sasvim blizu granice Mandelbrotovog skupa koje prelaze 2 tek nakon više od 50 koraka ovaj test će biti pogrešan, ali to možemo profiniti naknadno, povećavajući broj koraka.

Definiramo odgovarajući predikat

divergiraQ[c_]:=Last[iter[c,50]]>2

divergiraQ[-1+0.2I]

divergiraQ[1+I]

Problem je da je 50 iteracija jako puno i brojevi mogu postati enormni i blokirati računalo (vidi gore 1+i već nakon 5 iteracija). Bolje bi bilo da prekinemo iteriranje čim aposlutna vrijednost pređe 2 i proglasimo točku divergentnom. Definirajmo novi predikat divergiraQ2 u skladu s tim. Također ćemo kao vrijednosti ovog predikata umjesto True i False uzeti 1 i 0 radi lakšeg grafičkog prikaza.

divergiraQ2[c_]:=If[For[z=c;i=1,i<50&&Abs[z]<2,i++,z=(z^2)+c];i==50,0,1]

(Razložite si algoritam ove definicije!)

divergiraQ2[1+i]

DensityPlot[divergiraQ2[x+I y],{x,-2.01,0.7},{y,-1.1,1.1},PlotPoints->200,Mesh->False,Frame->False]

Ovo je ispao zapravo proceduralni program. Ako ipak želimo upražnjavati čisto funkcionalno programiranje, onda možemo definirati funkciju ovako (primijetite da nema pomoćnih funkcija ni varijabli):

mandelbrot[xmin_,xmax_,ymin_,ymax_,n_]:=DensityPlot[-Length[FixedPointList[#^2+(x+I y)&,x+I y,200,SameTest->(Abs[#2]>2.0&)]],{x,xmin,xmax},{y,ymin,ymax},PlotPoints->n,Mesh->False,Frame->False,AspectRatio->Automatic]

mandelbrot[-1.257,-1.247,0.044,0.049,400]

(Prikazan je dio Mandelbrotovog skupa poznat pod imenom "dolina morskog konjića" (seahorse valley).

Zadatak 4-6: Objasnite riječima algoritam ove funkcionalne inačice funkcije mandelbrot.

mandelbrot[0.278,0.2853,-0.007,-0.012,400]

Zadatak 4-7: Logističko preslikavanje je dano iteracijskom formulom (Subscript[x, n+1])= λ (Subscript[x, n])(1-(Subscript[x, n])).  Za male vrijednosti parametra λ, to preslikavanje za veliki n konvergira k jednoj vrijednosti. Kad λ raste, negdje blizu λ=3, dolazi do bifurkacije i iteracije više ne konvergiraju već skaču između dvije vrijednosti. S daljnjim rastom λ dolazi do nove bifurkacije i sustav se za veliki n ponavlja s periodom četiri, itd. Rezultat se može prikazati kao tzv. Feigenbaumov bifurkacijski dijagram

Nacrtajte ga. Mogući način: za dati λ kreirajte listu {(Subscript[x, 0]), (Subscript[x, 1]), ⋯,(Subscript[x, 400])} pa maknite prvih 100-200 članova da dođete u područje u kojem se sustav već stabilizirao. Nacrtajte to kao funkciju od λ, koristeći ListPlot.

4.3 Programiranje transformacijskim pravilima

Ideja ove vrste programiranja u Mathematici je promatrati npr. funkciju f(x)=sin(x), ne kao matematičko pridruživanje koja izvrijednjava izraz x i pridružuje mu sin(x), već kao transformaciju izraza x  u izraz sin(x).  To je u duhu uobičajene procedure računanja rukom — iz reda u red prelazimo tako da matematičke izraze transformiramo u skladu s nekim pravilima.

Globalna i lokalna pravila

Pridruživanje vrijednosti simbolima, koje smo obilato koristili, može se shvatiti kao globalno transformacijsko pravilo. Npr. pridruživanje

x=w;

kaže "kad god se u nekom izrazu pojavi x, transformiraj ga u w".

((Sin[x])^2)+myfun[(x^2)+1]

To pravilo je globalno u smislu da će biti upotrijebljeno u svakom izrazu. Ono je pridruženo simbolu x i to je pridruživanje na snazi do kraja Mathematica sesije tj. komande Quit.

?x

(Ovdje "Global" u Global`x nema veze s globalnošću pravila, već s globalnošću konteksta simbola x u što nećemo ovdje ulaziti. Koga zanima, neka pogleda pojam context  u manualu.)

Clear[x]

Lokalna pravila su ona koja se primjenjuju samo na jedan izraz, putem "/.". I njih smo često koristili:

myfun[(y^2)+3] /.y->w

3+y

?y

Vidimo da je simbol  y  "čist", u smislu da mu nije pridruženo nikakvo pravilo.

Ove dvije vrste pravila imaju, kako znamo, i inačice koje djeluju odgođeno, u smislu da se izrazi u samom pravilu izvrijednjuju tek prilikom primjene pravila. Za globalna pravila to se postiže upotrebom simbola := umjesto =, a za lokalna pravila upotrebom :> umjesto ->.

Ako primjenjujemo listu transformacijskih pravila, ona će se sva primijeniti istovremeno

{a,b}/.{a->b,b->a}

dok ukoliko primjenjuemo niz transformacijskih pravila, ona će se primijeniti konsekutivno, jedno po jedno.

{a,b} /.a->b/.b->a

Uzorci koji se ponavljaju

Ranije smo upoznali uzorke (patterns) koji se u definicijama funkcija i transformacijskih pravila koriste kao zamjena (wildcard) za proizvoljne izraze. Postoje i specijalni uzorci za izraze koji se ponavljaju

({{x_, bilo koji izraz}, {x__, bilo koji izraz koji se ponavlja jednom ili više puta}, {x___, bilo koji izraz koji se ponavlja nula ili više puta }})

To omogućuje definiranje funkcije s proizvoljnim brojem argumenata.

h[x__]:=Power[x]

h[a,b]

h[a,b,c]

h[4,4,4]

Primijetite da ovo odudara od matematičkog pojma funkcije koja uvijek mora imati dobro definiran broj argumenata, ali je sasvim u skladu s idejama transformacijskih pravila i uzoraka.

Ovo možemo iskoristiti npr. za definiranje funkcija s opcionalnim argumentima

data=ReadList["files/fit1.dat",{Number,Number}];

myPlot[fun_,iter_,opts___]:=Show[ListPlot[data,PlotStyle->{PointSize[0.02]},DisplayFunction->Identity],Plot[fun,iter,DisplayFunction->Identity,opts],DisplayFunction->$DisplayFunction]

myPlot[0.4(x^2), {x,0,3}]

Ova funkcija myPlot radi potpuno isto kao i obična funkcija Plot,  jedino što preko svakog dijagrama crta još i gore učitane točke. Tako je pogodna npr. za testiranje kvalitete fita. Treći, opcionalni argument opts___ služi sa prosljeđivanje željenih opcija funkciji Plot.

myPlot[0.5x,{x,0,6}, PlotStyle->{RGBColor[1,0,0],Dashing[{0.05,0.03}]}]

Ako želimo definirati funkciju kod koje opcionalni argumenti kad se izostave imaju neku defaultnu vrijednost, koristimo uzorak oblika:  uzorak_:<default>. Npr:

f4[x_,n_:1]:=(x^n)

f4[a+b,2]

f4[a+b]

Predikati

Već smo ranije upoznali predikate kao funkcije koje poprimaju samo vrijednosti True i False. Koristili smo ih kao argument funkcije Select[<lista>, <predikat>] za izbor elemenata liste koji zadovoljavaju neki uvjet. Predikate koristimo i za ograničavanje područja djelovanja transformacijskih pravila, odnosno, preciznije rečeno, za ograničavanje izraza koji će se spariti s određenim uzorkom. Tako će se npr. uzorak x_ /; IntegerQ[x] spariti samo sa izrazima koji su cijeli brojevi.

2x+3w /. (x_/;EvenQ[x]) ->(x^6)

Primijetite da  uzorak x_ ovdje nema nikakve veze sa simbolom x koji se pojavljuje u izrazu.

Ovo ima i jednostavniji oblik

2x+3w /. (x_?EvenQ) ->(x^6)

Uočite da u prvom obliku, nakon /;, dolazi predikat izvrijednjen za argument x (ime uzorka), a u drugom obliku nakon ? dolazi predikat kao funkcija.

Predikate možemo definirati i sami, kao bilo koju čistu funkciju koja poprima vrijednosti True ili False:

Select[Range[100],IntegerQ[(Sqrt[#])]&]

Ovdje smo pomoću predikata "da li je cijeli broj",  konstruirali predikat "da li je kvadrat cijelog broja".

Koristan predikat je FreeQ[ <izraz>, <forma> ] koji provjerava da li se forma (simbol, uzorak, ...) pojavljuje bilo gdje u izrazu ili njegovim dijelovima (podizrazima).

lst={1,Sin[w],(5^w),a,1,(c^2),(d^w),4,w};

Cases[lst,x_/;FreeQ[x,w]]

Dakle ovdje smo izabrali elemente liste koji su "slobodni od w".

Zadatak 4-8: Napravite operaciju komplementarnu  ovoj gore, tj. načinite listu čiji elementi će biti oni elementi liste lst  u kojima se pojavljuje simbol w.

Primjeri

Koristeći uzorke te svojstvo da funkcije prvo primjenjuju specifičnija pa onda općenitija transformacijska pravila možemo definirati funkcije složenog ponašanja. Npr.  logaritme Mathematica po defaultu ne ekspandira.

Log[a b]

Expand[%]

To je zapravo željeno ponašanje jer pravila ln(a b) = ln(a)+ln(b) i ln((a^n))=n ln(a)  ne vrijede za proizvoljne argumente (npr. prva relacija ne vrijedi ako su i a i b negativni brojevi). Ukoliko ipak želimo koristiti ove relacije, možemo definirati svoju funkciju

logExpand[Log[e1_ e2_]]:=Log[e1]+Log[e2]                      (* 1 *)
logExpand[Log[e_^n_]]:=n Log[e]                                           (* 2 *)

logExpand[Log[a b]]

logExpand[Log[(a^w)]]

No, za kompleksnije izraze logExpand radi loše. Npr. ne zna se distribuirati preko zbrajanja:

logExpand[Log[a b ]+2]

logExpand[Log[(a^2)b]]

U ovom drugom slučaju je problem u tome da  nakon jedne primjene pravila ekspanzija prestaje. Trebamo nekako izvesti da se pravila primjenjuju iznova sve dok ih se više nema na što primjeniti.

Clear[logExpand]
logExpand[Log[e1_ e2_]]:=logExpand[Log[e1]]+logExpand[Log[e2]](* 1 *)
logExpand[Log[e_^n_]]:=n logExpand[Log[e]]       (* 2 *)
logExpand[e1_+e2_]:=logExpand[e1]+logExpand[e2]   (* 3 *)
logExpand[c_ e_/;FreeQ[c,Log]]:=c logExpand[e]        (* 4 *)
logExpand[Log[e_]]:=Log[e]   (* 5 *)

Pravila su modificirana tako da glava logExpand preživljava svaku ekspanziju pa tako nakon jednog djelovanja slijedi novo, sve dok se ijedno pravilo može primijeniti. Na kraju zadnje pravilo (br. 5) miče glavu logExpand. Pravila br. 3 i 4 se brinu za distribuciju preko zbroja i produkta.

logExpand[3Log[a (b^g) c]+2Log[p (w^3)]]

Postoji i ugrađena Mathematica funkcija PowerExpand koja radi isto to, i neke druge stvari povrh toga.

PowerExpand[3Log[a (b^g) c]+2Log[p (w^3)]]

Zadatak 4-9: Korigirajte slijedeću manjkavost funkcije logExpand:

logExpand[Log[a b]+2]

Želimo ponašanje kao za PowerExpand :

PowerExpand[Log[a b]+2]

Zadatak 4-10: Definirajte svoju funkciju der[e_,x_] koja može derivirati proizvoljne polinome, dakle slabiju verziju ugrađene funkcije D[e_,x_].

Zadatak 4-11:   Poboljšajte funkciju der[e_,x_]tako da može izaći na kraj i s izrazima poput (x^2)(b x+c)ili (((x^2) +bx)/(a x)). Naputak: treba implementirati Leibnizovo lančano pravilo deriviranja produkta funkcija.

Kako pronaći prvo pojavljivanje dva ista elementa za redom u nekoj listi? Postoji zgodna tehnika koja koristi "triple-blank" ___ ponavljajuće uzorke:

lst={1,65,13,44,2,7,7,11,0};

lst/.{l___,a_,a_,r___}->a

To nam daje koji je to element. Da bi dobili njegovu poziciju u listi pokušajmo s

lst/.{l___,a_,a_,r___}->{a,Length[{l}]+1}

Problem je da se Length[{l}] izvrijednjava prije sparivanja uzorka l___ pa je to naprosto duljina liste od jednog elementa. Treba upotrijebiti odgođeno pravilo.

lst/.{l___,a_,a_,r___}:>{a,Length[{l}]+1}

Ovo radi dobro jer se sad funkcija Length[{l}] izvrijednjava tek nakon što je uzorak l___ već sparen tj. simbolu l je pridijeljena vrijednost — prvih pet elemenata liste lst.

Zadatak 4-12: Negdje u prvih 1000 decimala broja π pojavljuje se šest istih znamenaka za redom. Pronađite ih. Korisno je upotrijebiti funkcije ToString i Characters  i pomoću njih pretvoriti decimalni zapis broja π  u listu znamenaka.

Usput, funkcija ToString i njen par ToExpression mogu ponekad biti od koristi u konstrukciji raznih izraza.

StringJoin["Sin[","π]"]

ToExpression[%]

Skraćena infiks forma funkcije StringJoin je "<>"

ToExpression["2^"<>"3"]

Zadatak 4-13: Definirajte funkciju traziNtorku[list_, n_] koja će u zadanoj listi pronaći prvo pojavljivanje n istih elemenata zaredom.

Zadatak 4-14: Prim-brojevi blizanci su parovi prim-brojeva koji se razlikuju za dva, poput (3,5) ili (17,19). V. Brun je 1919. dokazao da suma recipročnih vrijednosti prim-brojeva blizanaca konvergira

B≡((1/3)+(1/5))+((1/5)+(1/7))+((1/11)+(1/13))+((1/17)+(1/19))+⋯

Konstanta B se zove Brunova konstanta i iznosi približno B=1.902160583104. Ovako preciznu vrijednost je vrlo teško dobiti, no napišite komandu brun[n_] koja izračunava B koristeći listu od prvih n prim-brojeva. (Koristite ugrađenu funkciju Prime.) Inače, zanimljivo je da je upravo računalno određivanje ove konstante ukazalo na bug u prvoj generaciji Pentium procesora.

4.4 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:

<<"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 4-15: 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`", "VectorAnalysis`"]

Zadatak 4-16: 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 4-17: Modificirajte gornji paket tako da radi striktno u 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 VectorFieldPlots`

Zadatak 4-18: 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}.

Naputak: Postupajte kao u zadatku s kraja biljeznice 2-MatematikaS.nb. Za pocetak kod iscrtavanja uzmite neki fiksirani dio x-y ravnine, a kasnije mozete dodati svojstvo da program crta samo onaj dio ravnine u kojem se cestica giba. Rezultat za centralni potencijal treba biti npr.

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