Pierwsze kroki po środowisku Mathematica

autor : G. Siudem

Notatnik niniejszy jest ilustracją najciekawszych (i początkowo najbardziej użytecznych) funkcjonalności środowiska Mathematica. Stanowi też przykład tego jak estetycznie i czytelnie zredagować swój notatnik, tworząc w Mathematice sprawozdanie ze swojej pracy. Aby sformatować swój notatnik podobnie do tego należy wykorzystać opcje w menu Format (po więcej szczegółów proponuję doczytać hasło MouseOver  w pomocy lub samemu pobawić się możliwościami menu Format). Zagadnienia będą prezentowane tematycznie, zgodnie z dziedziną matematyki.

Uwaga! Pierwszym i najważniejszym źródłem informacji w pracy ze środowiskiem Mathematica jest jego pomoc. Aby ją wywołać wciśnij [F1].

Analiza matematyczna

Potrzebujesz zróżniczkować bardzo złożoną funkcję? Mathematica zrobi to za Ciebie!
Najpierw zdefiniujmy potrzebną funkcję

F[x_]:=Sin[x Log[2+Cos[Exp[x]]]]

F[x]

D[F[x],x]

Można obliczać też pochodne wyższych rzędów.

D[F[x],{x,4}]

Skoro znamy już pochodne rozważanej funkcji, przekonajmy się jak wygląda jej wykres

Plot[F[x],{x,-10,5},PlotRange->{-1.1,1.1},Frame->True,FrameStyle->Directive[18,Bold,FontFamily->"Courier",Background->White],Background->White,PlotStyle->{Thick,Red},FrameLabel->{x,"F[x]"}]

Rozważmy teraz funkcję G zmiennej x parametryzowaną przez β

G[x_,β_]:=(1/(2β))Log[((Cosh[β(x+1)])/(Cosh[β(x-1)]))]

G[x,β]

Uwaga! aby szybko wpisywać greckie litery alfabetu można posłużyć się skrótami - np. α można szybko uzyskać piszac [Esc]a[Esc] lub [Esc]alpha[Esc]. W podobny sposób szybko można uzyskać postać (a/b) wciskając [Ctrl]+/, (a^2) natomiast [Ctrl]+6.

Wyznaczmy funkcję pierwotną funkcji G

Integrate[G[x,β],x]

ten sam efekt osiągniemy stosując poniższą notację [Esc]intt[Esc]

∫G[x,β]\[DifferentialD]x

Analogicznie obliczamy całki oznaczone

Integrate[Sin[x k],{x,0,2π}]

można też w ten sposób [Esc]intt[Esc], a następnie dodać granice dolną [Ctrl]+[Shift]+- i górną [Ctrl]+[Shift]+5.

(∫_0^(2π))Sin[k x]\[DifferentialD]x

Mathematica potrafi też obliczać granice

Limit[Tanh[x],x->∞]

rozwijać w szereg

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

Należy zwrócić szczególną uwagę na fakt, że wszystkie te operacje (poza rysowaniem wykresu) wykonywane były symbolicznie, nie numerycznie. Oczywiście Mathematica umożliwia także obliczenia numeryczne. Proszę zapoznać się na początek z pomocą MouseOver oraz MouseOver.

Algebra liniowa

Czy Mathematica potrafi dokonywać operacji na macierzach? Oczywiście. Najpierw stwórzmy macierz

macierz={{1,2},{3,4}}

W środowisku Mathematica macierze są listami list, jest to wygodne programistycznie, ale utrudnia oglądanie wyników. Zapiszmy więc macierz w wygodniejszej postaci

MatrixForm[macierz]

Obliczmy wyznacznik takiej macierzy

Det[macierz]

zbadajmy jej widmo

Eigenvalues[macierz]

Ile tak naprawdę wynoszą te liczby? Sprawdźmy.

N[Eigenvalues[macierz]]

Rozwiążmy teraz równanie algebraiczne

Solve[{x+a y+1==b,y-x==2},{x,y}]

Zwróćmy uwagę, że wykorzystano tutaj podwójny znak równości! W środowisku Mathematica = jest czym innym niż == i czym innym niż ===. Po szczegóły różnic odsyłam do pomocy.

Manipulacje

Od wersji 6.0 środowisko Mathematica posiada niesamowicie potężne narzędzie - funkcję Manipulate. Zilustrujmy jej zastosowanie badając dynamikę tzw. rodziny logistycznej.

Załóżmy, że chcemy opisać populację drożdży. Niech (Subscript[x, n])oznacza liczność populacji tych istot w n-tej chwili czasu. Drożdże mogą się rozmnażać (w tempie α) , załóżmy zatem na początek, że nasz model będzie wyglądał następująco

(Subscript[x, n+1])=α (Subscript[x, n])
    (Subscript[x, 0])=z

Jest to jednak zły model, gdyż drobny rachunek lub wykorzystanie poniższego kodu pozwala znaleźć zwarte rozwiązanie powyższego równania rekurencyjnego

RSolve[x[n+1]==α x[n]&&x[0]==z,x[n], n]

Widzimy zatem, że dla α>1 osiągamy populacje o przerażająco dużym rozmiarze, a dla 0<α<1 populacja wyginie. Poprawmy więc nasz model dodając do niego czynnik śmierci

(Subscript[x, n+1])=α (Subscript[x, n])(m-(Subscript[x, n]))
    (Subscript[x, 0])=z

gdzie m jest pojemnością środowiska, jeśli zostanie przekroczona nasze drożdże uprzykszają sobie życie będąc stłoczone zbyt ciasno, co skutkuje ograniczeniem tempa rozwoju. 

Można oczywiście przyczynę śmierci tłumaczyć bardziej biologicznie - drożdże w wyniku procesów życiowych produkują alkohol etylowy, jeśli jest ich zbyt dużo, stężęnie alkoholu skutkuje śmiercią tych istot.

Jak jednak, niezależnie od przyczyny śmierci drożdży, badać taki bardziej złożony układ? Z pomocą przyjdzie nam funkcja MouseOver, która umożliwia wielokrotne iterowanie rozważanej funkcji, zaczynając od warunku początkowego.

H[x_]:=α x(m-x);
Nest[H,z,1]
Nest[H,z,2]

Ustalmy, dla ustalenia uwagi, że m=1, można to traktować jako wyrażenie x w jednostkach pojemności środowiska. Stoimy zatem przed pytaniem jak zmiany α wpływają na zachowanie się naszych drożdży

m=1;
α=1;
z=0.5;
ListPlot[NestList[H,z,20],PlotRange->{0,1},Frame->True,FrameStyle->Directive[18,Bold,FontFamily->"Courier",Background->White],Background->White,PlotStyle->PointSize[0.01],FrameLabel->{n,"(Subscript[x, n])"}]

Widzimy, że dla ustalonych wartości parametrów populacja stabilizuje się (można zgadywać, że na zerze, co jest dla drożdży smutne), co jednak jeśli zmienimy wartość parametrów α albo z? Jak to wpłynie na dynamikę. Tutaj z pomocą przyjdzie nam wspomniana funkcja Manipulate.

H[x_,α_]:=α x(1-x);
Manipulate[
ListPlot[NestList[H[#,α]&,z,80],PlotRange->{0,1},Frame->True,FrameStyle->Directive[18,Bold,FontFamily->"Courier",Background->White],Background->White,Joined->True,Mesh->All,PlotStyle->PointSize[0.01],FrameLabel->{n,"(Subscript[x, n])"}],{α,1,3.999,0.001},{{z,0.5},0,1}]

Teraz można wykryć, że jeśli α przekroczy wartość 3 pojawiają się oscylacje w liczebności populacji, a jeśli dopuścimy jeszcze większe α, wówczas zachowania całego układu stają się całkowicie niezrozumiałe! [proszę sprawdzić samemu i dla α>3.8 zmienić nieznacznie wartość z]

Co się dzieje z drożdżami podczas zmiany parametru α? Dlaczego tak prosta funkcja H dała tak nieregularne zachowanie? Czy to jest chaos? 

Po odpowiedzi na te i inne pasjonujące pytania zapraszam na Układy Dynamiczne :)