Podstawy Wolfram Mathematica

Korzystam tutaj z wersji 14 Mathematiki. W przypadku korzystania z wersji starszej lub nowszej, pewne tresci teoretycznie mogą wyglądać inaczej.
Studenci UW mają prawo do uzyskania licencji: https://www.fuw.edu.pl/oprogramowanie-licencje.html#mathematica

Domyślnym typem komórki jest komórka z kodem Wolfram Language, w którym operuje Mathematica. Wykonujemy ją przez Shift+Enter. Możemy dodawać też inne typy – wybierając z paska na górze ekranu bądż używając skrótów (np. Alt+7 dla zmiany typu na komórkę tekstową). Cały notatnik możemy wykonać przez wybranie z paska menu “Evaluation->Evaluate Notebook”. Może to być przydatne gdy komórki nawzajem do siebie się odwołują (czyli dość często).

(*To jest komentarz.
Może mieć wiele linii (ale musi być w jednej komórce). *)

Możemy też wywoływać Wolfram|Alpha – zaczynając linię od “=” bądź “==”. Pierwsza opcja jest przydatna ze względu na to, że zwraca także poprawny input w języku Wolfram Language – stąd jest doskonała do poznawania nowych komend. Druga opcja zwraca pełny output, jaki uzyskalibyśmy na stronie Wolfram|Alpha, także z możliwością uzyskania rozwiązania krok po kroku.

Przejdźmy jednak do samej Mathematiki

Zmienne w Mathematice są dynamicznie typowane – oznacza to, że nie mają na twardo przypisanego typu (np. liczba całkowita, liczba zmiennoprzecinkowa, tekst) i możemy zapisywać do nich cokolwiek chcemy. Możemy też kazać Mathematice “zapomnieć” o zmiennej (wtedy traktowana jest jako symbol i zaznaczona jest na niebiesko).

x=5

x

x="A teraz zmienna x zawiera tekst."

x=.

x

Jak widać, Mathematica domyślnie wyświetla wynik każdego wykonanego polecenia. Możemy to “wyciszyć” kończąc polecenie średnikiem.

q=Pi;

q=Pi

Domyślnie Mathematica będzie starała się wyświetlać wynik w postaci zawierającej ułamki oraz stałe. Wartość numeryczną uzyskamy przez polecenie “N” – obejmując nim wyrażenie bądź wykonując je na wyniku korzystając z notacji “//”:

2*Pi/3

N[2*Pi/3]

2*Pi/3//N

Liczbę cyfr możemy kontrolować dodatkowym argumentem:

N[2*Pi/3,50]

2*Pi/3//N[#,50]&

Ta druga opcja wykorzystuje tzw. “czyste funkcje” (“pure functions”), którymi nie będziemy się w tym miejscu zajmować. Jedyne, co musimy wiedzieć, to to, że “#” oznacza parametr funkcji (czyli w tym wypadku całe wyrażenie), a “&” na końcu umożliwia skorzystanie z funkcjonalności wspomnianych czystych funkcji. Dla bardziej dociekliwych – czysta funkcja to taka, która zawsze zwraca ten sam wynik dla tych samych parametrów (czyli nie ma jakichś ukrytych zależności od np. zmiennych globalnych) oraz jej działaniem jest jedynie podanie wyniku – nie prowadzi do żadnych efektów ubocznych (jak zmiana zmiennych globalnych, wypisanie czegoś itp.).

Precyzja liczb

Mathematica, oprócz wartości liczbowych, zapamiętuje także precyzję z jaką są one wyrażone – czyli miarę błędu względnego danej wartości (po szczegóły odsyłam do dokumentacji Mathematiki).  Domyślnie jest to tak zwana precyzja maszynowa, czyli taka w której potrafi natywnie operować komputer (wykorzystywana jest tzw. podwójna precyzja, w której mantysa zapisana jest w 53 bitach – stąd w postaci dziesiętnej MachinePrecision=53/Log2[10]). Liczby całkowite oraz stałe takie jak Pi (które Mathematica jest sobie w stanie wyznaczyć dowolnie dokładnie w razie potrzeby) mają nieskończoną precyzję.

Precision[1]

Precision[1.0]

N[MachinePrecision]

N[53/Log2[10]]

Precyzję podawanych liczb możemy ustawić ręcznie. Należy jednak pamiętać, że wyjście poza precyzję maszynową może spowolnić nasze obliczenia.

0.33`18

Precision[0.33`18]

Odwoływanie się do poprzednich wyników

Jeśli  chcemy  dalej  wykorzystywać wynik jakiegoś polecenia, możemy odwołać się do niego na różne sposoby

w=1/(x^2+1)-x^2/(x-1)

Przez znak “%” – oznacza on poprzedni wynik

2*%

Przez odpowiedni numer komórki – odradzam, bo komórek nie musimy wykonywać po kolei oraz możemy wykonywać je wielokrotnie, przez co może się zmienić numeracja:

2*Out[16]

Najlepiej po prostu zapisać wyrażenie w jakiejś zmiennej (jak zrobiliśmy to wyżej) i odwołać się do niej:

2*w

Możemy też skopiować wynik do dalszego wykorzystania – albo zaznaczyć go i Ctrl+c, Ctrl+v,  albo, jeśli nie chcemy wyszukanego formatowania, użyć “InputForm”:

2*(-((x^2)/(-1+x))+(1/(1+(x^2))))

InputForm[w]

Dość często wykorzystywaną operacją jest podstawienie:

-(x^2/(-1 + x)) + (1 + x^2)^(-1)/.x->y

-(x^2/(-1 + x)) + (1 + x^2)^(-1)/.x->0

Podstawowe operacje matematyczne

Granice:

Limit[1/x,x->0]

Limit[1/x,x->0,Direction->"FromBelow"]

Limit[1/x,x->0,Direction->"FromAbove"]

Całkowanie:

Integrate[Exp[-x^2],x]

Integrate[Exp[-x^2],{x,0,Infinity}]

Różniczkowanie:

D[x^2/(x+1),x]

D[x^2/(x+1),{x,3}]

Rozwinięcie  w  szereg (“Normal” porzuca informację o pominiętych kolejnych wyrazach “O[]”, przez co możemy swobodnie operować szeregiem jak zwykłym wyrażeniem matematycznym):

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

Series[Sin[x],{x,0,10}]//Normal

SeriesCoefficient[Sin[x],{x,0,7}]

Przekształcanie wyrażeń

Wyrażenia zwracane przez Mathematikę bywają dość rozwlekłe. Możemy jej kazać podjąć próbę ich uproszczenia:

w=D[Sin[x]/(x+1),{x,2}]

Simplify[w]

FullSimplify[w]

Ogólne polecenia upraszcające “Simplify” i “FullSimplify” bywają bardzo kosztowne (bo Mathematica próbuje całej długiej listy potencjalnych uproszczeń). Często sporo możemy osiągnąć za pomocą prostszych poleceń.

Rozkład na ułamki proste:

(2*x+1)/(x^2+3*x-5)*(x+1)//Apart

Zebranie do wspólnego mianownika:

(2*x^2)/(1+x)^3-(4*x)/(1+x)^2+2/(1+x)//Together

Rozwinięcie wyrażenia:

(a+b*x)^2//Expand

Tu jeszcze raz przypominam, że omawiane polecenia możemy wykonywać też przez potraktowanie całego wyrażenia jako ich argumentu:

Expand[(a+b*x)^2]

Faktoryzacja – zebranie do postaci iloczynowej:

Factor[a^2+2*a*b*x+b^2*x^2]

Zebranie  członów  stojących  przy  tych  samych  potęgach  podanej  zmiennej  (lub  wyrażenia) :

Collect[(a+b*x)^2*(c+d*x)^2,x]

Tablice i macierze:

Tablice wstawiamy w Mathematice za pomocą nawiasów klamrowych. Odwołujemy się do ich elementów za pomocą nawiasów kwadratowych:

m={a,b}

m[[2]]

Można to robić również za pomocą funkcji “Table”:

m=Table[i^2,{i,1,10}]

Kolejne  wymiary dodajemy przez zagnieżdżanie kolejnych nawiasów klamrowych:

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

Liczbę wymiarów redukujemy z kolei poleceniem “Flatten”:

Flatten[m]

Wygodne przy wizualizacji jest tu polecenie “MatrixForm”:

MatrixForm[m]

Można też posłużyć się interfejsem – na górnym pasku: Insert->Table/Matrix:

m=({{1, 2, 3}, {4, 5, 6}, {7, 8, 9}})

Znów odwołujemy się przez podwójne nawiasy kwadratowe:

m[[2,3]]

Na takich macierzach możemy już wykonywać działania. Np. wyznacznik:

Det[m]

Mnożenie  macierzy :

m.m

Inaczej:

Dot[m,m]

Macierz odwrotna (jeśli istnieje!):

Inverse[m]

m2=({{1, 1, 1}, {1, 2, 1}, {1, 1, 3}});

m2i=Inverse[m2]

m2.m2i//MatrixForm

Iloczyn skalarny macierzy z wektorem:

v={1,2,3}

m.v

Wektora z wektorem:

v.v

Iloczyn wektorowy dwóch wektorów:

Cross[v,v]

v2={1,2,4};

Cross[v,v2]

Zagadnienie własne (podaje wektory własne oraz wartości własne):

m=({{1, 2, 3}, {4, 5, 6}, {7, 8, 9}})

{vals,vecs}=Eigensystem[m]//Simplify

Sprawdźmy czy to działa! Wybierzmy sobie np. drugi wektor:

m.vecs[[2]]-vals[[2]]*vecs[[2]]//Simplify

Równania:

Mathematica umożliwia łatwe rozwiązywanie równań:

s=Solve[a*x^2+b*x+c==0,x]

Polecenie “Solve” zwraca wynik w postaci tablicy z podstawieniami.

x/.s

Po zadziałaniu nim na zmienną x, jak widać dostaliśmy tablicę z wszystkimi możliwymi podstawieniami. Możemy wybrać konkretne:

x/.s[[2]]

Możemy jeszcze sprawdzić czy wyniki są prawidłowe:

a*x^2+b*x+c/.s//Simplify

Oczywiście Mathematica radzi sobie także z trudniejszymi równaniami:

w=a*x^4+b*x^3-c*x^2+d*x+e

s=Solve[w==0,x]

Wzory te są już dość pokaźne. W praktyce przy bardziej złożonych problemach zwykle rozwiązujemy równania dla konkretnych wartości współczynników a,b,c...

w=x^4+2*x^3-3*x^2+5x+8

s=Solve[w==0,x]

Mathematica lubi wyrażać pierwiastki równań wyższych stopni w nieco zakamuflowany sposób. Jeśli chcemy wymusić na niej jawne wypisanie pełnych wzorów, możemy zrobić to w następujący sposób:

ToRadicals[s]

Nadal są to dość skomplikowane kombinacje liczb całkowitych (tego się przecież raczej spodziewaliśmy). Możemy więc posłużyć się znanym nam “N”:

N[s]

Możemy  jeszcze sprawdzić czy rozwiązania są poprawne:

w/.s//Simplify

Skoro i tak w takich wypadkach interesują nas konkretne wartości liczbowe, a nie pełne formuły, możemy też – zamiast robić to naokoło – od razu rozwiązać równanie numerycznie:

NSolve[w==0,x]

Możemy  też  rozwiązywać  równania  macierzowe (m.v==w) :

m={{1,2},{3,4}};
w={{5,6},{7,8}};

v=LinearSolve[m,w]

Łatwo możemy sprawdzić, że rozwiązanie jest poprawne:

m.v-w

Równania różniczkowe również są w naszym zasięgu:

DSolve[y'[x]-y[x]==a Cos[x],y[x],x]

Nie podając parametru przy “y”, dostajemy rozwiązanie w postaci wspomnianej już “czystej funkcji”:

s=DSolve[y'[x]-y[x]==a Cos[x],y,x]

(y'[x]-y[x])/.s//Simplify

Ponownie mamy wariant numeryczny. Zwraca nam on wynik w postaci zinterpolowanej:

NDSolve[{y'[x]-y[x]==Cos[x],y[0]==0},y,{x,0,10}]

Wykresy:

Do tworzenia wykresów w Mathematice służą polecenia o nazwach mówiących same za siebie. W podstawowej wersji są dość wygodne do sprawdzania przebiegu funkcji. Po więcej opcji odsyłam do dokumentacji poszczególnych poleceń.

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

Plot3D[Sin[x]*Cos[y],{x,-5,5},{y,-5,5}]

DensityPlot[Sin[x]*Cos[y],{x,-5,5},{y,-5,5}]

Poniższe polecenie zakomentowane dla bezpieczeństwa – na niektórych słabszych kartach graficznych DensityPlot3D może powodować zawieszenie się całego systemu. Dla świętego spokoju warto więc zapisać przedtem pracę. :-)

(*DensityPlot3D[Cos[x]*Cos[y]*Cos[z],{x,-5,5},{y,-5,5},{z,-5,5}]*)

To dobry moment na wspomnienie o poleceniu “Manipulate” – objęcie nim dowolnego innego polecenia powoduje dodanie interaktywnego interfejsu pozwalającego na... cóż, manipulowanie jego wynikiem. Tu przykład z manipulowaniem dokładnością wykresu (wraz ze wzrostem “a” powinien stawać się “gładszy”):

Manipulate[DensityPlot[Sin[x]*Cos[y],{x,-5,5},{y,-5,5},MaxRecursion->a],{a,0,15,1}]

Dyskretne dane, dopasowywanie i interpolacja

Możemy też rysować dane dyskretne. Wczytajmy je najpierw z pliku. Podajemy nazwę pliku z danymi oraz to jak Mathematika ma te dane interpretować (tu ma się spodziewać w każdej linii trzech liczb rzeczywistych). Jako przykład wziąłem dane z wykładu o Gnuplocie, czyli punkty leżące w przybliżeniu na paraboli, wraz z niepewnościami (które nas w tej chwili nie interesują). Pierwsze polecenie ma na celu ustawienie roboczego katalogu (czyli tego, w którym będą domyślnie poszukiwane i zapisywane pliki) na katalog tego notatnika.

SetDirectory[NotebookDirectory[]]

dat=ReadList["parab.dat",{Real,Real, Real}];

MatrixForm[dat]

Pozbądźmy się najpierw trzeciej kolumny. Możemy to zrobić na więcej niż jeden sposób. Krótko i elegancko – wybrać wszystkie wiersze, ale tylko pierwszą i drugą kolumnę.

data=dat[[All,1;;2]];

Można także użyć znanego nam już polecenia “Table”:

data=Table[{dat[[i,1]],dat[[i,2]]},{i,1,Length[dat]}];

MatrixForm[data]

Do narysowania dyskretnych danych służy polecenie “ListPlot”:

plt1=ListPlot[data]

plt2=ListPlot[data,Joined->True]

Tak jak w Gnuplocie, możemy dopasować jakąś funkcję matematyczną do danych:

fit=Fit[data,{1,x,x^2},x]

plt3=Plot[fit,{x,Min[data[[All,1]]],Max[data[[All,1]]]}]

Aby wyrysować zarówno dane (“ListPlot”) oraz dopasowaną funkcję (“Plot”), możemy użyć polecenia “Show”:

Show[plt1,plt3]

Oprócz dopasowywania funkcji, możemy próbować również interpolacji. Wówczas interpolująca funckcja jest dokładnie równa interpolowanym wartościom  tam gdzie są one podane, a przybliża je jakąś funkcją (zwykle wielomianem Hermita lub funkcją B-Spline) poza nimi. Z jednej strony wyższy stopień wielomianu może lepiej oddać skomplikowaną zmienność funkcji, ale z drugiej przesada może prowadzić do patologicznych oscylacji. Można to zaobserwować zwiększając opcję “InterpolationOrder”. Domyślna wartość “3” wydaje się sensowna do typowych zastosowań.

inter=Interpolation[data,InterpolationOrder->3,Method->"Spline"];

plt4=Plot[inter[x],{x,-4,4}]

Show[plt1,plt4]

Pętle i warunki:

Jak przystało na język programowania, Wolfram Language zawiera pętle i instrukcje warunkowe:

Do[
Print[i]
,{i,1,5}]

x1=0;

x2=If[x1!=0,1,-1]

Własne funkcje:

Własne funkcje możemy definiować na różne sposoby. Pierwszy jest bardzo formalny i tworzy wspomnianą już wcześniej “czystą funkcję”.

f=Function[{x,a},x^a]

f[x,2]

To samo, ale w innej notacji:

f= #1^#2 &

f[x,2]

Czyste funkcje mają choćby tę zaletę,że Mathematica rozumie czym są ich pochodne:

f'[x,2]

f=.

Drugi sposób, nieco bardziej elastyczny, wygląda tak:

f[x_,a_]:=x^a

f[x,2]

Parametry  funkcji  zaznaczone  są  podkreślnikiem  po  nazwie . Możemy  ich  potem  swobodnie  używać . Znak  ":="  oznacza  "odroczone przypisanie" – czyli  za  każdym  razem, gdy  Mathematica  natrafi  na  to  wyrażenie, przelicza  je  na  nowo . W przypadku funkcji jest to najbezpieczniejsze podejście. Oczywiście  z  drugiej  strony – są  przypadki, gdy  nie  będziemy  tego  chcieli  (np . by  zaoszczędzić  na  czasie) .

Korzystając z tej metody definiowania, nie dostajemy automatycznie definicji pochodnej:

f'[x,2]

Nadal możemy jednak swobodnie korzystać z operacji różniczkowania:

D[f[x,2],x]

Jeżeli nasza funkcja jest bardziej skomplikowana, być może będziemy chcieli użyć zmiennych tymczasowych, ale jednocześnie nie chcemy by powodowało to niepotrzebnie  globalne definicje. Wówczas możemy posłużyć się środowiskiem “Module”. W nawiasie klamrowym na początku podajemy zmienne, których zakres chcemy, by ograniczał się do tego jednego środowiska.

f[ran_]:=Module[{u,w},
u=0;w=1;
Do[
tmp=w;
w=u+w;
u=tmp;
,{i,1,ran}];
u]

f[10]

Zmienne “u” i “w” nie są zdefiniowane poza ciałem funkcji. Za to “tmp” pozostaje z ostatnią wartością, jaka była w nim zapisana. Moglibyśmy to oczywiście łatwo naprawić dodając je do nawiasu klamrowego na początku “Module”.

{u,w, tmp}

Eksport danych do pliku

Nierzadko chcemy wyeksportować wyniki naszych obliczeń do pliku tekstowego, na przykład w celu dalszej ich obróbki. Najprościej chyba zrobić to poleceniem “Export”

t=Table[{i,Sin[i]//N},{i,0,100}];
Export["wyniki.txt",t];

W ten sposób w pliku “wyniki.txt” zapisana zostanie tabela z wartościami funkcji sin(x) . Gdy jednak zajrzymy do niego, może nie podobać się nam domyślny format (z nawiasami klamrowymi i przecinkami). W takim wypadku możemy doprecyzować format. Np. “Table” spowoduje zapisanie w prostym formacie gdzie rekordy oddzielone są pustą przestrzenią:

Export["wyniki2.txt",t,"Table"];

Jeśli wolimy spacje (lub inny separator) zamiast tabulatorów, również możemy to doprecyzować:

Export["wyniki3.txt", t, "Table","FieldSeparators"->" "];

Odsyłam do dokumentacji po więcej szczegółów. W tym miejscu warto jednak zwrócić uwagę na szerokie zastosowanie tego polecenia. Możemy w analogiczny sposób również zapisać wykres do obrazka:

Export["wykres.png",ListPlot[t,Joined->True,InterpolationOrder->3]];

Pytania kontrolne

Jak zdefiniować w Mathematice macierz 2x2, zawierającą w pierwszym wierszu wartości 1 i 2, a w drugim 3 i 4?

Jakiego polecenia użyć do narysowania wykresu dyskretnych danych?

Jakiego polecenia użyć by zróżniczkować funkcję f dwukrotnie po q?

Jak przypisać do symbolu x wartość numeryczną 0.3 tak, by Mathematica traktowała ją jako liczbę o precyzji 36 cyfr?

Jakiego polecenia użyć do zinterpolowania danych “dane.dat” wielomianem 5 stopnia?

Jakich poleceń można użyć by stworzyć własną funkcję? Podaj dwa.

Jakiego polecenia użyć, by numerycznie rozwiązać równanie?