2008-10-01

2. Elementarna i viša  matematika na Mathematici

2.1 Algebarske manipulacije

Moć paketa za simboličku matematiku, poput Mathematice, leži u sposobnosti manipulacije simboličkim izrazima. Kao prvo, svakom izrazu možemo radi lakšeg snalaženja pridijeliti ime.

izraz1=(((b+a))^3)

Kod davanja imena mogu se upotrebljavati i grčka slova, indeksi itd. kako bi se postigla sličnost sa standardnom matematičkom notacijom ((Subscript[β, 1]), (OverHat[ℱ]), …), ali takva imena su čest izvor grešaka i ne preporučaju se početnicima.

Mathematica ne provodi skoro nikakve operacije na izrazima dok je eksplicitno ne instruiramo. Recimo da želimo raspisati ("ekspandirati") izraz1 koristeći binomni teorem

Expand[izraz1]

To bi mogli i na ruke. Međutim:

izraz2=(((a+2b+3c))^3)(((x+y))^3)

Expand[izraz2]

Ukoliko ne želimo vidjeti ispis rezultata neke komande, stavimo na kraj ";".

izraz3=Expand[(izraz2^3)];

Length[izraz3]

Komanda Short ispisuje samo početak i kraj izraza te naznačuje broj ispuštenih članova pa tako možemo dobiti neku ideju o rezultatu bez da zakrčimo ekran:

Short[izraz3]

Inverzna operacija ekspanziji je faktorizacija:

Factor[izraz3]

Ovo faktorizira cijeli izraz u najmanji broj faktora. Faktoriziranje konkretne varijable tj. prikazivanje izraza u obliku polinoma u toj varijabli:

Collect[izraz2, y]

Traženje koeficijenta uz određeni izraz:

Coefficient[izraz2,(y^2)]

Stavljanje na zajednički nazivnik:

Together[(a/(1-a))+(a/(1+a))]

Jedna od najčešće korištenih funkcija je Simplify  koja obično daje najjednostavniji mogući zapis izraza

Simplify[%]

Zadatak 2-1: Uzmite izraz3=(a+b)((c+d x)x+f (x^2))  i korištenjem gornjih komandi natjerajte Mathematicu da ga prikaže u slijedeća tri oblika:
(a) (a+b) c x+(a+b) (d+f) (x^2)
(b)  a c x+b c x+a d (x^2)+b d (x^2)+a f (x^2)+b f (x^2)
(c)  (a+b) x (c+d x+f x)

Funkcija Simplify zna i trigonometrijske identitete

Simplify[((Cos[x])^2)+((Sin[x])^2)]

(1/(Cot[b]))//Simplify

Zadatak 2-2: Uvjerite se da vrijedi identitet
 (((cos^3)x + (sin^3)x)/(cosx + sin x))=1-cos x sin x
Pazite na sintaksu funkcija:  (sin^n)x  je  Sin[x]^n

2.2 Jednadžbe

U Mathematici kao znak jednakosti u jednadžbama stoji "==",  jer je uobičajeni znak "=" , kako smo gore već vidjeli, rezerviran za pridjeljivanje vrijednosti simbolima (imenima). Rješavanje jednadžbi se izvodi ugrađenom funkcijom Solve:

rjes=Solve[2(x^2)-1==0, x]

Valja primijetiti:
1. Prilikom poziva funkcije Solve treba eksplicitno naznačiti po kojoj varijabli se traži rješavanje.
2. Mathematica je pronašla oba rješenja kvadratne jednadžbe.
3. Rezultat je ispisan u obliku liste {...}takozvanih transfromacijskih pravila ("transformation rules"). Ta se pravila mogu onda primijeniti na proizvoljni izraz pomoću simbola /.  koji znači "primijeni na prethodni izraz transformacijska pravila koja slijede":

x+a /.rjes

Tako provjeru dobivenih rješenja možemo napraviti primjenom transformacijskih pravila  na zadanu jednadžbu

2(x^2)-1==0 /. rjes

(Uvrštavanjem rješenja u jednadžbu dobije se jednadžba poput 0==0 koji je uvijek zadovoljena pa odgovara logičkoj istini i Mathematica je automatski transformira u logičku Booleovu konstantu True.)

Transformacija izraza pomoću simbola /. je vrlo korisna operacija. Više riječi o njoj će biti u odjeljku 3.5, a sada ćemo samo pokazati njenu najčešću upotrebu: pridjeljivanje vrijednosti simbolima u nekom izrazu:

(((a+b-c))^4) //Expand

%/.b->c

%/. a->2

Funkcija Solve može rješavati i sustave jednadžbi, ako se kao argumenti zadaju liste jednadžbi odnosno varijabli:

Solve[{(x^2)+(y^2)==1,x-2y==0},{x,y}]

Zadatak 2-3: Riješite sustav jednadžbi:
(x^4)+(a^4)=1
(x^2)+(a^2)=1
Koliko ima rješenja?

Funkcija Solve nastoji naći simboličko rješenje, no to je ponekad nemoguće. Npr. jednadžbe viših stupnjeva uopće ne moraju imati sva rješenja u zatvorenom obliku.

Solve[9(x^6)+4(x^4)+3(x^3)+x-17==0, x]

x=1 je trivijalno rješenje, ali ostalih pet rješenja Mathematica ne može napisati u zatvorenom obliku već ih zapisuje na poseban način u koji sad nećemo ulaziti. U takvim slučajevima možemo pribjeći numeričkom rješavanju koje radi  funkcija NSolve:

NSolve[9(x^6)+4(x^4)+3(x^3)+x-17==0, x]

(Često simboličke funkcije u Mathematici imaju svoj numerički pandan koji ima isto ime do na početno slovo N: Integrate ↔ NIntegrate, Sum ↔ NSum, ... )

S transcendentalnim jednadžbama situacija je nešto teža. Naime, i Solve i NSolve su predviđene samo za rješavanje polinomnih jednadžbi.

NSolve[2 ArcTan[x]==(x^2), x]

Funkcija kojom zaista numerički rješavamo jednadžbe je FindRoot

?FindRoot

Vidimo da FindRoot traži da mu zadamo početnu točku od koje kreće tražiti nultočku nekom od numeričkih metoda poput Newtonove.

FindRoot[2 ArcTan[x]==(x^2), {x,0}]

FindRoot[2 ArcTan[x]==(x^2), {x,1}]

Pronašli smo dva rješenja. Daljnji pokušaji daju uvijek jedno od ova dva:

FindRoot[2 ArcTan[x]==(x^2), {x,5}]

FindRoot[2 ArcTan[x]==(x^2), {x,1000}]

FindRoot[2 ArcTan[x]==(x^2), {x,-5}]

Ovo zadnje rješenje je zapravo 0 jer Mathematica, ako joj se ne kaže drugačije, radi numeriku s preciznošću od otprilike (10^(-15)).

Options[FindRoot]

N[MachinePrecision]

$MachinePrecision

Ove varijable s "$" na početku su tzv. environment  varijable koje definiraju okolinu rada Mathematice.

$Version

Važna varijabla je $Path koja specificira direktorije po kojima će Mathematica tražiti datoteke koje želimo učitati.

$Path

FindRoot[2 ArcTan[x]==(x^2), {x,-5}, WorkingPrecision->50]

Vidimo da se s povećanjem radne preciznosti rješenje još više približilo nuli. Funkcija Chop stavlja u egzaktnu nulu sve brojeve koji su joj bliži od (10^(-10)) i korisna je za povećanje preglednosti računa.

Chop[%]

Možemo se uvjeriti da su ova gore dva rješenja zaista jedina, tako da skiciramo grafove funkcija 2 arctg(x) i (x^2) i vidimo da se sijeku na samo dva mjesta. Koristimo funkciju Plot o kojoj će kasnije biti više riječi.

Plot[{2ArcTan[x],(x^2)}, {x,-5,5}]

Zadatak 2-4:

Riješite jednadžbu:
tg(x)-(x/10)=0

Zadatak 2-5:

Nađite (kompleksno!) rješenje jednadžbe:
sin(x)=2
i provjerite uvrštavanjem. Pokušajte pronaći još neko rješenje

Zadatak 2-6:

Riješite nejednadžbu:
(x^2)+x-12<0

Zadatak 2-7:

Pronađite pozicije minimuma i maksimuma Besselove funkcije (Subscript[J, 1])(x) koji su najbliži točki x=0.

2.3 Matematička analiza

Mathematica sadrži sve standardne operacije koje se uče u matematičkoj analizi, poput limesa, nizova, redova i diferencijalnog računa. Limesi se izvrijednjavaju funkcijom Limit

Limit[((Sin[x])/x),x->0]

Razvoj neke funkcije po nekoj varijabli oko neke točke do nekog reda radi Series:

Series[Sin[x],{x,0,5}]

Zadatak 2-8: Odredite prva četiri člana u razvoju funkcije  arctg(x) oko točke x=∞.

Zadatak 2-9

Kolika je relativna pogreška koju radimo ako za izvrijednjavanje arctg(5) koristimo četveročlani red dobiven u gornjem zadatku? Da li bi greška bila manja da smo upotrebljavali četveročlani red oko x=0?

Deriviranje se radi funkcijom D:

D[(E^(2x))Cos[3x],x]

Višestruko deriviranje:

D[(E^(2x))Cos[3x],{x,4}]

Deriviranje (nepoznate) opće funkcije f(x), korištenjem Leibnitzovog lančanog pravila:

D[(x^2)f[x],x]

Simboličko neodređeno integriranje:

Integrate[2 (E^(2 x)) Cos[3 x]-3 (E^(2 x)) Sin[3 x],x]

Primijetite da se konstanta integracije podrazumijeva bez navođenja. Alternativno, postoji i zapis koji grafički bolje izgleda:

∫x Log[x]\[DifferentialD]x

Simboličko rješavanje određenih integrala:

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

Neki integrali su preteški ili se naprosto ne daju prikazati u zatvorenoj formi pa Mathematica, kao i uvijek kad nešto ne zna izračunati, vraća zadani izraz:

Integrate[ArcTan[BesselJ[1,x]],{x,0,1}]

Tada nam ostaje numerička integracija pomoću NIntegrate

NIntegrate[ArcTan[BesselJ[1,x]],{x,0,1}]

Zadatak 2-10: Izračunajte neodređeni integral
∫(((x^2)+3)/((x^5)+(x^4)-x-1))dx
i onda provjerite dobiveni rezultat deriviranjem i algebarskim manipulacijama.

Zadatak 2-11: Izračunajte određene integrale
((∫_0)^1) sin(sin(x))dx            ((∫_0)^π) sin(sin(x))dx
simbolički i numerički.

Zadatak 2-12: Izračunajte dvostruki određeni integral
((∫_0)^1)dx((∫_0)^(1-x))dy (x+y)(Sqrt[x (y^3)])

Zadatak 2-13: Izračunajte 
∫dx(∫dy (x+y)(Sqrt[x (y^3)])) 
i onda provjerite dobiveni rezultat deriviranjem i algebarskim manipulacijama.

2.4 Linearna algebra

Vektori su u Mathematici reprezentirani kao liste

vec1={1,1,2}
vec2={2,2,4}

Skalarni produkt vektora se dobije točkom "." (ili funkcijom Dot),

vec1.vec2

a vektorski produkt funkcijom Cross

Cross[vec1,vec2]

Matrice su liste odgovarajućih red-vektora

mat={{1,2,1},{4,3,3},{9,1,7}}

Za prirodniji ispis koristimo funkciju MatrixForm

MatrixForm[%]

Za upis matrice možemo koristiti i izbornik Insert->Table/Matrix

mat=(({{1, 2, 1}, {4, 3, 3}, {9, 1, 7}}))

S matricama i vektorima možemo raditi sve ono na što smo navikli. Važno je samo znati da se množenje izvodi točkom "." (koja označava nekomutativno množenje), a ne zvjezdicom "*" ili razmakom kao kod običnog množenja koje je komutativno. Zaboravljanje točke je najčešći uzrok pogrešaka kod rada s matricama.

Množenje matrica (uočite točku!):

mat.mat

MatrixForm[%]

Množenje matrice i vektora (uočite točku!):

mat.vec1

Inverz matrice:

Inverse[mat]

Provjera

mat.% //MatrixForm

Primijetite da smo upotrijebili postfiksni oblik funkcije MatrixForm.

Moguće je i mijenjati pojedine elemente matrice. Sjetite se da se element liste dobiva pomoću dvostrukih uglatih zagrada. (Njihov elegantan zapis se postiže sekvencama Esc-[[-Esc  i  Esc-]]-Esc

mat[[1,1]]=0

mat //MatrixForm

Zadatak 2-14: Odredite svojstvene vrijednosti (eigenvalues) i svojstvene vektore (eigenvectors) matrice
(({{2.3`, 4.5`}, {6.7`, -1.2`}}))
i provjerite za jedan od tih vektora da je zaista svojstveni eksplicitnim množenjem.

Zadatak 2-15:

(1)Koristeći funkcije Table i RandomReal (tj. Random u verzijama Mathematice<6.0) kreirajte 3x3 matricu sa slučajnim realnim brojevima između 0 i 10. (2)Zamijenite središnji element matrice simbolom x.  
(3)Invertirajte novonastalu matricu te je 
(4)pomnožite s neinvertiranom i provjerite da je rezultat jedinična matrica.

Zadatak 2-16:

Kakvo množenje matrica dobivamo ako izostavimo točku u A.B?

2.5 Crtanje grafova

2D

Glavna komanda za dvodimenzionalno crtanje je Plot.

Plot[Sin[x],{x,0,2π}]

Poput većine Mathematica komandi i Plot ima svoje opcije pomoću kojih možemo podešavati njegovo ponašanje. Evo ih zajedno s njihovim defaultnim vrijednostima:

Options[Plot]

Za detaljno značenje pojedinih opcija pogledajte dokumentaciju. Evo nekoliko primjera najzanimljivijih opcija:

Plot[Sin[x],{x,0,2π},AxesLabel->{"x","sin(x)"},PlotRange->{{0,2},{0,1.5}}]

Redosljed kojim specificiramo opcije je nebitan.

Za specificiranje svojstava samih grafova funkcija koristi se opcija PlotStyle prema donjem primjeru koji pokazuje i kako  na istoj slici istovremeno prikazati više grafova.

Plot[{Sin[x],Cos[x]},{x,0,2π},PlotStyle->{{Thick,Dashed,Red},{Thickness[0.01],RGBColor[0,0,1]}}]

Za ponovno prikazivanje grafa možemo koristiti funkciju Show, koja omogućava i modificiranje opcija grafa.

Show[%,AspectRatio->1.5,PlotRange->{0,1},Frame->True,FrameLabel->{"x","f(x)","naslov","f(x)"}]

Grafove često želimo upotrijebiti i izvan Mathematica bilježnice, npr. u nekom članku ili prezentaciji. Eksportiranje grafova i drugih objekata postiže se uporabom funkcije Export.Slijedeći primjer bi trebao rezultirati datotekom sinus.gif  u radnom direktoriju. (Za popis mogućih grafičkih formata vidi u manualu funkciju Export.)

Export["sinus.gif",%%]

Zadatak 2-17: Pronađite ovu datoteku na hard disku, prikažite je na ekranu pa onda izbrišite. (Radni direktorij možete pronaći izdavanjem komande Directory[]).

Zadatak 2-18: Nacrtajte graf funkcije ln(x) između x=0.5 i x=1.5 bez ikakvih oznaka, okvira i osi, dakle samo liniju.

Zadatak 2-19:

Koristeći funkciju ParametricPlot nacrtajte kružnicu (na ekranu treba izgledati baš kao kružnica, a ne kao elipsa).

3D

Glavna komanda za trodimenzionalno crtanje je Plot3D.

Plot3D[(x^4)+(y^4),{x,-2,2},{y,-2,2}]

Za promjenu kuta gledanja u Mathematica verzijama≥ 6.0 naprosto "uhvatimo" graf mišem i okrećemo,  a u ranijim verzijama koristimo opciju ViewPoint, čiju vrijednost možemo odrediti pomoću dijaloga Input->3D ViewPoint Selector ...

Plot3D  ima većinu opcija kao i Plot plus još dosta novih

Show[%,Boxed->False,AxesLabel->{"x-os","y-os","z-os"},Lighting->{{"Directional",RGBColor[1,.7,.1],{{5,5,4},{5,5,0}}}}]

2.6 Vektorska analiza

Osim velikog broja ugrađenih ("built-in") funkcija, s Mathematicom dolazi i određen broj paketa s dodatnim funkcijama. Za njihovu upotrebu prvo je potrebno učitati odgovarajući paket. Npr. paket VectorAnalysis definira funkcije poput gradijenta, divergencije, rotacije i laplasijana

<<"VectorAnalysis`"

Defaultni koordinatni sistem s kojim ovaj paket radi je kartezijev s koordinatama Xx, Yy i Zz:

Coordinates[]

Neka se česticu u ravnini nalazi u potencijalu U(x,y)= -4(x^2)(y^3).  Odredimo odgovarajuću silu (Overscript[F, ->])=-∇U:

sila=-Grad[-4 (x^2) (y^3)/.{x->Xx,y->Yy}]

Dakle, polje sile je (Overscript[F, ->])=8 x (y^3) (Overscript[i, ->])+ 12 (x^2) (y^2) (Overscript[j, ->]).

Vektorsko polje možemo nacrtati korištenjem slijedećeg paketa:

<<"VectorFieldPlots`"

VectorFieldPlot[Take[sila,2], {Xx,-2,2},{Yy,-2,2},ColorFunction->Hue]

(Inače, dosta specijaliziranih dodatnih paketa se može skinuti s  http://library.wolfram.com/infocenter/MathSource/ )

2.7 Grafički prikaz i fitanje podataka

Često je potrebno neke podatke dobivene npr. mjerenjima grafički prikazati i približno opisati nekom funkcijom (tzv. fitanje). Prvo je potrebno učitati podatke u Mathematica listu oblika

{{(Subscript[x, 1]),(Subscript[y, 1])},{(Subscript[x, 2]),(Subscript[y, 2])},...}

Neka u radnom direktoriju bude poddirektorij files s datotekom fit1.dat  s podacima složenim u dva stupca. Izravan ispis datoteke je moguć ovako (na MS Windows platformi je možda nužno zamijeniti "/" -> "\")

FilePrint["files/fit1.dat"]

Importiranje ovih podataka u listu za daljnju obradu najelegantnije se obavlja funkcijom Import,

Import["files/fit1.dat"]

dok funkcija ReadList omogućuje nešto precizniju kontrolu

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

Za grafički prikaz možemo koristiti funkcije ListPlot i ScatterPlot3D.

ListPlot[podaci1,PlotStyle->{PointSize[0.015]}]

(Default veličina točaka je obično presitna za ugodno gledanje.)

Fitanje se izvodi funkcijom

Fit[<lista s podacima>, <lista funkcija kojima fitamo>, <lista varijabli>]

koja izvodi fitanje metodom najmanjih kvadrata na linearnu kombinaciju specificiranih funkcija. Npr. linearni fit ovih podataka na funkciju  f(x)=a x+b   je

Fit[podaci1,{1,x},{x}]

Grafičku provjeru kvalitete fita možemo provesti zajedničkim iscrtavanjem ovog pravca i podataka. (Primjetite kako cut'n'paste gornjeg rezultata u donju input čeliju pokazuje da Mathematica interno reprezentira brojeve preciznije nego što ih ispisuje. Takav cut'n'paste je redovito nepoželjan jer lako dovodi do grešaka i nepreglednih izraza! Bolje je izrazima davati imena!)

Show[Plot[0.020159523556311256`+0.5130561168421428` x,{x,0,5}],ListPlot[podaci1,PlotStyle->{PointSize[0.015]}]]

Zadatak 2-20: Podatke iz datoteke fit2.dat nafitajte na polinom trećeg reda, te na f(x)=a sin(x) te nacrtajte sve na jednom grafu. Neka sinusna funkcija bude u boji radi lakšeg snalaženja.

Fit radi s linearnom kombinacijom zadanih funkcija. Jedno moguće poboljšanje je korištenje funkcije FindFit (pojavljuje se tek u Mathematici 5.0) koja može raditi nelinearni fit.

fitfun=a Sin[b x];
nelfit=FindFit[podaci2,fitfun,{a,b},{x}]

Show[Plot[fitfun/.nelfit,{x,0,5}],
ListPlot[podaci2,PlotStyle->{PointSize[0.02]}]]

Koji je fit najbolji? Trebala bi pažljiva analiza koju ćemo ostaviti za neki kasniji zadatak.

Zadatak 2-21:

Pronađite neka svoja mjerenja koja ste napravili na Fizičkom praktikumu i načinite fitanje nekom razumnom funkcijom. Nacrtajte grafički prikaz i stavite na njega legendu.

2.8 Diferencijalne jednadžbe

Osnovna Mathematica funkcija koja traži analitička rješenja običnih diferencijalnih jednadžbi je DSolve. Riješimo pomoću nje jednadžbu gušenog harmoničkog oscilatora

(((d^2)y)/(dt^2))+2 (dy/dt)+4y=0

Važno je voditi računa o tome da je nužno specificirati (u uglatim zagradama) nezavisne varijable o kojima funkcije ovise, u ovom slučaju vrijeme t

DSolve[y''[t]+2y'[t]+4y[t]==0,y[t],t]

Primijetite da smo dobili opće rješenje s dvije nepoznate konstante C[1] i C[2] koje treba odrediti iz početnih uvjeta. Ukoliko želimo partikularno rješenje za neke konkretne početne uvjete, zadamo te uvjete kao dodatne jednadžbe:

sol=DSolve[{y''[t]+2y'[t]+4y[t]==0,y'[0]== 1,y[0]==1},y[t],t]

Plot[(1/3) (E^(-t)) (3 Cos[(Sqrt[3]) t]+2 (Sqrt[3]) Sin[(Sqrt[3]) t]),{t,0,5}]

Neke jednadžbe se ne mogu riješiti analitički, pa moramo pribjeći numeričkom integriranju funkcijom NDSolve. Npr. na kolegiju Diferencijalne jednadžbe i dinamički sustavi numerički se analizira tzv. Van der Polova jednadžba.

sol=NDSolve[{x''[t]-(1-((x[t])^2))x'[t]+x[t]==0,x'[0]==1,x[0]==1},x[t],{t,0,15}]

Rezultirajući objekt je specijalni oblik funkcije, InterpolatingFunction, koji nije analitički, ali je možemo izvrijedniti u nekoj konkretnoj točki ili nacrtati.

Plot[Evaluate[x[t]/.sol],{t,0,15}]

InterpolatingFunction objekt možemo i derivirati

dersol=D[Evaluate[x[t]/.sol],t]

Plot[Evaluate[dersol],{t,0,15}]

NDSolve može riješiti i neke jednostavnije parcijalne diferencijalne jednadžbe.

Zadatak 2-22: Pomoću ParametricPlot nacrtajte graf derivacije gornjeg rješenja Van der Polove jednadžbe x'(t) u ovisnosti o rješenju x(t)  (tzv. fazni dijagram).

Zadatak 2-23: Neka na česticu u ravnini djeluje sila (Overscript[F, ->])=8 x (y^3) (Overscript[i, ->])+ 12 (x^2) (y^2) (Overscript[j, ->])   tj. potencijal U(x,y)= -4(x^2)(y^3) kao u odjeljku 2.6. Nacrtajte potencijal U(x,y) i pomoću Plot3D i pomoću ContourPlot. Nadalje, rješavajući Newtonovu diferencijalnu jednadžbu gibanja, pronađite putanju čestice. Nacrtajte tu putanju, a onda je pokušajte superponirati na gornje dijagrame potencijala. Igrajte se malo s početnim uvjetima i uvjerite se da je sve OK. Izvrijednite ukupnu energiju za razna vremena i uvjerite se da je sustav zaista konzervativan tj. da je energija konstantna.

Zadatak 2-24: Isto kao gornji zadatak, ali neka putanje ne budu kontinuirane linije već točke koje su ekvidistantne u vremenu tako da njihov razmak opisuje trenutnu brzinu tijela. Veći razmak odgovara većoj brzini.