Programowanie Środ Matlab

Z Skrypty dla studentów Ekonofizyki UPGOW

(Różnice między wersjami)
(Wykresy standardowe)
(Macierze)
 
(Nie pokazano 3 wersji pomiędzy niniejszymi.)
Linia 380: Linia 380:
Różnica czasów wykonania powyższych operacji bez- oraz z wektoryzacją jest także imponująca. Jeżeli szybkość działania przygotowywanego przez Ciebie programu w środowisku MATLAB ma istotne, wektoryzacja powinna pomóc. Naturalnie, o ile wykorzystywane w programie algorytmy poddają się wektoryzacji.
Różnica czasów wykonania powyższych operacji bez- oraz z wektoryzacją jest także imponująca. Jeżeli szybkość działania przygotowywanego przez Ciebie programu w środowisku MATLAB ma istotne, wektoryzacja powinna pomóc. Naturalnie, o ile wykorzystywane w programie algorytmy poddają się wektoryzacji.
Możesz także przyspieszyć działanie pętli '''for''' poprzez wstępną inicjalizację (prealokację) macierzy lub wektorów, w których będą skladowane obliczenia wykonywane w "długich" pętlach. Zobacz kod poniżej,
Możesz także przyspieszyć działanie pętli '''for''' poprzez wstępną inicjalizację (prealokację) macierzy lub wektorów, w których będą skladowane obliczenia wykonywane w "długich" pętlach. Zobacz kod poniżej,
-
ten kod używa funkcji zeros, aby dokonać prealokacji wektora tworzonego w pętli for. Dzięki temu pętla będzie się niekiedy wykonywać szybciej, zobaczmy przyklad  
+
ten kod używa funkcji '''ones'''(), aby dokonać prealokacji wektora tworzonego w pętli for. Dzięki temu pętla będzie się niekiedy wykonywać szybciej, zobaczmy przyklad  
<source lang="matlab">
<source lang="matlab">
Linia 990: Linia 990:
[[Plik:Hist2.png]]
[[Plik:Hist2.png]]
-
Gdy argumentem '''hist'''() jest macierz o wymiarach ''m''*''p'', '''hist'''() traktuje kolumny tej macierzy jako wektory i zwróci histogram/macierz o rozmiarach 10*p. Dalsze przykłady użycia '''hist'''() znajdziesz np. w skrypcie * [[MKZR:Dodatek#Histogramy]] .
+
Drugi parametr '''hist'''() może być wektorem, nazwę go ''x''. Wówczas będzie tyle kontenerów, ile elementów wektora ''x'' czyli ''nk'' = '''length'''(x) kontenerów. Pokrywać będą przedział od najmniejszej do największej wartości wśród jego elementów. Typowy przyklad,
 +
 
 +
<source lang="matlab">
 +
dane = rand(100, 1);     
 +
xp = min(dane);
 +
xk = max(dane);
 +
krok = 0.1;
 +
x = xp: krok: xk
 +
hist(dane, x)
 +
</source>
 +
 +
Gdy pierwszym argumentem '''hist'''() jest macierz o wymiarach ''m''*''p'', '''hist'''() traktuje kolumny tej macierzy jako wektory i zwróci histogram/macierz o rozmiarach 10*p. Dalsze przykłady użycia '''hist'''() znajdziesz np. w skrypcie * [[MKZR:Dodatek#Histogramy]] .
Na posiedzeniach rad nadzorczych, zarządów korporacji, kierownictw przedsiębiorstw, spółek, ...itp. bardzo często spotykaną formą prezentacji rozmaitych danych ekonomicznych jest tzw. wykres kołowy/procentowy (ang. ''pie chart''). Pokazuje on, jak podzielony jest przysłowiowy tort (''pie'' - placek owocowy). Ile procent przychodu firmy poszło np. na: podatki, ubezpieczenia, inwestycje, materiały do produkcji, media, koszty osobowe, ... no i najbardziej interesujące, jaki procent przychodu stanowi dochód/zysk firmy. Taki wykres można sporządzić np. w arkuszu kalkulacyjnym Excel, teraz pokażę, jak to zrobić w MATLAB (zob. ''help'' pie ). W najprostszej wersji, jako argument funkcji '''pie'''() wystarczy podać wektor reprezentujący jakieś dane. Funkcja '''pie'''() oblicza procentowy udział każdego skladnika w ich sumie. W przykładzie poniżej ze wzgledów dydaktycznych przygotowałem wektor danych tak, aby suma jego składników wynosila 100. Wobec tego wartość każdego składnika to procentowy wkład do sumy wszystkich skladników, sprawdźmy to na wykresie
Na posiedzeniach rad nadzorczych, zarządów korporacji, kierownictw przedsiębiorstw, spółek, ...itp. bardzo często spotykaną formą prezentacji rozmaitych danych ekonomicznych jest tzw. wykres kołowy/procentowy (ang. ''pie chart''). Pokazuje on, jak podzielony jest przysłowiowy tort (''pie'' - placek owocowy). Ile procent przychodu firmy poszło np. na: podatki, ubezpieczenia, inwestycje, materiały do produkcji, media, koszty osobowe, ... no i najbardziej interesujące, jaki procent przychodu stanowi dochód/zysk firmy. Taki wykres można sporządzić np. w arkuszu kalkulacyjnym Excel, teraz pokażę, jak to zrobić w MATLAB (zob. ''help'' pie ). W najprostszej wersji, jako argument funkcji '''pie'''() wystarczy podać wektor reprezentujący jakieś dane. Funkcja '''pie'''() oblicza procentowy udział każdego skladnika w ich sumie. W przykładzie poniżej ze wzgledów dydaktycznych przygotowałem wektor danych tak, aby suma jego składników wynosila 100. Wobec tego wartość każdego składnika to procentowy wkład do sumy wszystkich skladników, sprawdźmy to na wykresie
Linia 1101: Linia 1112:
[[Plik:bes.gif]]
[[Plik:bes.gif]]
-
Do tej pory, nie rozpraszać uwagi pomijałem takie "detale" jak opis osi wykresu, dobór koloru lub znaków "piszących" linie wykresu, ...itp. Oczywiście, jeżeli chcemy zrobić porządny, czytelny wykres/rysunek, musimy go staranie opisać.
+
Do tej pory, aby nie rozpraszać uwagi pomijałem takie "detale" jak opis osi wykresu, dobór koloru lub znaków "piszących" linie wykresu, ...itp. Oczywiście, jeżeli chcemy zrobić porządny, czytelny wykres/rysunek, musimy go staranie opisać.
Do opisu wykresu służy paleta funkcji, których parametrami są najczęściej łańcuchy znaków, najczęściej używane to;
Do opisu wykresu służy paleta funkcji, których parametrami są najczęściej łańcuchy znaków, najczęściej używane to;

Aktualna wersja na dzień 12:32, 26 kwi 2011

Środowisko Matlab

MATLAB jest obecnie bardzo rozbudowanym narzędziem obliczeniowo-programistycznym, w tym rozdziale przedstawię tylko wybrane, podstawowe elementy środowiska MATLAB.

Spis treści


Informacje podstawowe

MATLAB jest obecnie bardzo rozbudowanym narzędziem obliczeniowo-programistycznym, w tym podrozdziale przedstawię tylko wybrane, podstawowe elementy środowiska MATLAB. Omawiam zagadnienia i przyklady, które są wykonywalne również w niekomercyjnym środowisku GNU Octave, ma to znaczenie, zwłaszcza dla studentów ekonofizyki. W szczególności nie omawiam ścieżek dochodzenia i "wyklikiwania" odpowiednich komend/funkcji w okienkowym interfejsie MATLAB.

Żródłowe informacje "on line" o środowisku MATLAB oraz GNU Octave znajdziesz pod adresami:

http://www.mathworks.com/products/matlab

http://www.gnu.org/software/octave


Po uruchomieniu MATLAB w oknie pojawia się znak

>>

(ang. prompt) zachęcający do pisania instrukcji. Napiszmy

>> 3 + 4 + 3

Po wpisaniu instrukcji należy nacisnąć klawisz Enter/Return, co oznacza polecenie jej wykonania przez interpreter MATLAB. Wówczas na ekranie pojawi się wynik działania, przypisany do standardowej nazwy ans (ang. answer czyli odpowiedź)

ans = 10

Zakończenie instrukcji średnikiem powoduje, iż rezultat jej działania nie pojawi się na ekranie, jest to wygodne, gdy np. nie interesują nas pośrednie wyniki sekwencji obliczeń.

Instrukcję składającą sie z bardzo dużej liczby wyrażeń, znaków, itp - przez to niewygodną do manipulowania, nieczytelną - możemy dowolnie, według naszego uznania, podzielić na większa liczbę linii. Po napisaniu trzech kropek ... w dogodnym dla nas miejscu, przechodzimy do nowej linii i piszemy dalej. Interpreter potraktuje wyrażenie zapisane w dwóch lub więcej liniach, jako jedną całość. Z drugiej strony, w jednej linii możemy napisać kilka instrukcji, wystarczy je oddzielać przecinkami np.

>> 2 + 4 + 10, 100 – 30, 1000 – 1, 4 + 4 – 8

Na ekranie pojawi się sekwencja wyników,

ans = 16

ans = 70

ans = 999

ans = 0

Do nawigacji pomiędzy poprzednimi i następnymi instrukcjami służą klawisze oznaczone strzałkami (góra, dół). Zadania wykonywane przez MATLAB można w każdej chwili zatrzymać wywołując wbudowaną funkcję pause() - do momentu naciśnięcia klawisza dowolnego znaku - ,pause(n) zatrzymać na n sekund, definitywnie przerwać naciskając kombinację klawiszy ctrl c. Instrukcja cls czyści okno, quit lub exit powoduje wyjście z MATLAB.

Znak % (procent) rozpoczyna krótki komentarz, czyli wszystkie znaki napisane po nim, aż do końca linii, są pomijane przez interpreter, służą tylko programiście do opisu. Obszerniejszy komentarz tworzymy umieszczając wielolinijkowy tekst pomiędzy znakami %{ oraz %} ,

%{

aaaaaaaa

bbbbbbbb

cccccccc

%}

Zmienne w MATLAB nie wymagają wcześniejszego deklarowania, tworzone są „w locie”, w momencie ich pierwszego użycia (inicjalizacji, nadania wartości), np. instrukcja

>> x = 1.4

powoduje utworzenie zmiennej o nazwie x, jej wartość wynosi 1.4. Dokładniej mówiąc, interpreter rezerwuje odpowiednią ilość komórek pamięci (bajtów) dla zmiennej o nazwie x, do tych komórek pamięci wprowadza wartość, którą w systemie dziesiętnym przedstawiamy jako 1.4. Nazwa zmiennej może zawierać tylko litery alfabetu łacińskiego (duże lub małe), cyfry od 0 do 9, znak podkreślenia, musi przy tym zaczynać się od litery. Naturalnie, nazwa zmiennej nie może się pokrywać ze słowem kluczowym MATLAB. Lista słów kluczowych MATLAB pojawi się na ekranie monitora, po komendzie iskeyword. Wartość danej zmiennej można zobaczyć podając jej nazwę i naciskając klawisz Enter,

>> a = 10;

>> a

>> a = 10

Przeczytałeś już jedna stronę o MATLAB, a jeszcze nie było żadnego programu, najwyższy czas to nadrobić. Zgadnij, od czego zaczniemy, .......? Tak, tradycji musi stać się zadość, czysta klasyka, wg niepisanego zwyczaju od tego zaczyna się nauczanie najrozmaitszych języków progamowania. Tradycja jest tradycją, nie będziemy się wyłamywać, oto prawdziwy "evergreen", osławiony, kultowy program "hello, world", w wersji dla środowiska MATLAB

disp(' hello, world')

To już cały program, jedna linijka kodu. Wpisz go w lini komend i naciśnij Enter, oznacza to rozkaz wykonania kodu, w rezultacie na ekranie pojawi się napis

>> Hello world

>>

Ten najprostszy program to po prostu wywołanie dostępnej w MATLAB (wbudowanej) funkcji disp(). Wypisuje ona na ekran jej argument (to, co wewnątrz nawiasów półokrągłych) i przechodzi do nowej linii. Wypiszmy coś innego

>> disp(' Wynik testu ')

>> Wynik testu

>>

>> disp(' Karol 4.0 Basia 4.5 Andrzej 3.5')

>> Karol 4.0 Basia 4.5 Andrzej 3.5

>>

W MATLAB wszystkie zmienne są traktowane jako tablice, zmienna skalarna traktowana jest jako tablica jednoelementowa. MATLAB operuje na liczbach zmiennoprzecinkowych podwójnej precyzji. Domyślnie wyświetla jednak tylko 4 cyfry po kropce, aby zapewnić czytelność - taki format nazywany jest w MATLAB format short. Użytkownik może ustalić inny format wyświetlania liczb, pokażę to na przykładzie predefiniowanej w MATLAB zmiennej o nazwie pi:

>> format long  % pokazywać liczby z precyzją 14 miejsc po kropce dziesiętnej

>> pi

>> 3.14159265358979


>> format longE  % precyzja 14 miejsc po kropce, reprezentacja potęgowa

>> pi

>> 3.141592653589793e+000 . Przy definiowaniu zmiennej łańcuchowej (tekstowej) korzysta się z pojedynczych apostrofów, np.

>> napis = ‘Hej !! ‘

>> napis

>> Hej !!

Zmienna typu łańcuchowego może zawierać nazwę komendy/funkcji, taką komendę uruchamia funkcja eval( ), np.

>> napis = ‘cos(sin(0)) ‘

>> eval(napis)

>> ans = 1

Funkcja input( ) służy do komunikacji z użytkownikiem programu, wypisuje ona na ekran treść swojego argumentu (łańcuch znaków) i MATLAB oczekuje na reakcję użytkownika. Oto jej typowe zastosowanie

>> n = input( ‘Wpisz n: ‘ )

>> Wpisz n:  % system czeka na wpisanie wartości n

Po wpisaniu za dwukropkiem np. liczby 10 i jak zwykle naciśnięciu klawisza Enter, zobaczymy na ekranie wynik podstawienia

>> n = 10

Drugi przykład,

>> pierwszy_atom = input (‘ Nazwa atomu numer jeden: ‘ , ‘s’ )

>> Nazwa atomu numer jeden:  % system czeka na wpisanie nazwy

Wpisanie drugiego parametru funkcji input( ) w postaci ‘s’ oznacza, że interpreter ma potraktować wprowadzoną daną jako ciąg znaków. Pisząc za dwukropkiem np. tlen dostajemy

>> pierwszy_atom = tlen

Funkcja num2str( ) zamienia liczbę na ciąg znaków, str2num( ) zamienia ciąg znaków na liczbę.

Komenda who powoduje wyświetlenie nazw wszystkich zmiennych zdefiniowanych w bieżącym programie (przestrzeni roboczej), whos podaje listę zmiennych wraz z dodatkowymi szczegółami (precyzja, rozmiar macierzy). Usunięcie wszystkich zmiennych powoduje komenda clear all, wybraną zmienną, np. o nazwie a, usuwamy instrukcją clear a.

Użycie komendy diary podczas sesji MATLAB powoduje, iż od tego momentu cały przebieg aktualnej sesji MATLAB jest zapisywany w pliku o nazwie diary. Komenda,

>> diary moja_sesja

zapisuje przebieg sesji w pliku, który nazwałem moja_sesja. Rejestrację sesji zawiesza komenda diary off, diary on nakazuje wznowienie rejestracji sesji.

MATLAB posiada liczny arsenał funkcji standardowych, np. polecenie

>> exp(1)

kreuje odpowiedź

ans = 2.7183

Kolejny przykład,

>> exp(sin(0)) + cos(0)

ans = 2

Wywołując funkcję rand() otrzymujemy zmiennoprzecinkową liczbę losową o rozkładzie równomiernym z przedziału < 0, 1 >, natomiast rand(n, m) tworzy macierz o wymiarach n * m takich liczb losowych. Po każdym uruchomieniu MATLAB algorytm generowania liczby losowej w funkcji rand() startuje od tego samego, domyślnego ziarna (ang. seed). Aby tego uniknąć można wywołać rand() jak następuje

>> rand('state',sum(100*clock));

Tym sposobem „seed” jest wybierany na podstawie wskazań zegara systemowego, a więc za każdym wywołaniem inny.

Użyteczną funkcją jest polyval( ), która ułatwia posługiwanie się wielomianami. Należy podać dwa parametry wejściowe; wektor współczynników wielomianu oraz wektor wartości, dla których ten wielomian ma być obliczony. Jeżeli interesują mnie wartości wielomianu, np. f(x) = x3 + 4x2 – 2x -3 w punktach x = {0, 1, 2, ….9, 10}, wystarczą trzy instrukcje

>> x = 0 : 10; c = [1 4 -2 3]; f = polyval(c, x);

Wykres tego wielomianu f(x) w przedziale <0, 10> otrzymamy po wywołaniu plot( )

>> plot(x, f)

Dotychczasowe przykłady były wypisywane bezpośrednio, w wierszu poleceń po znaku zachęty. Takie postępowanie jest czytelne, niezbyt uciążliwe i efektywne tylko przy niewielkiej sekwencji instrukcji, zajmujących kilka linii. W praktyce kod programu wygodniej przygotowywać, modyfikować bądź poprawiać w oddzielnym pliku. Jeżeli zapiszemy ten plik z rozszerzeniem *.m, na przykład moj_kod.m, to MATLAB potraktuje go jako skrypt, tak zwany m – plik. Wpisanie w linii komend nazwy m-pliku (bez rozszerzenia)

>> moj_kod

powoduje wykonywanie wszystkich instrukcji zapisanych w tym pliku, jak gdyby były one wpisywane po kolei w linii poleceń. Komenda what podaje listę m-plików zapisanych w bieżącym katalogu.

Dostęp do dokumentacji MATLAB zapewnia polecenie doc. Pojawia się lista zagadnień (menu). Wybieramy interesujący nas temat, zagłębiając się w bardziej szczegółowe poziomy dokumentacji. Warto także skorzystać z komendy help, na przykład

>> help sin

wyświetla informację o funkcji sinus. Komenda

>> help plot

prowadzi do szczegółowego opisu opcji funkcji plot() i pokrewnych. Komenda lookfor, np.

>> lookfor sqrt

powoduje przeszukiwanie plików pomocy szukając łańcucha wpisanego po lookfor, tutaj sqrt. Nawiasem mówiąc, jest to nazwa wbudowanej funkcji pierwiastek kwadratowy.

Macierze

W MATLAB zmienne domyślnie są traktowane jako macierze (tablice), niezbędne jest poznanie podstawowych zasad posługiwania się takimi zmiennymi. Macierz składa się z elementów ułożonych w wierszach (poziomo) i kolumnach (pionowo). Macierz o m wierszach i n kolumnach zawiera m * n elementów. Szczególnym przypadkiem jest wektor, czyli macierz o tylko jednym wierszu i n kolumnach (1 * n) lub tylko jednej kolumnie i m wierszach (m * 1). Skalar to po prostu tablica jednoelementowa, dla której m = n = 1. Istnieje kilka sposobów definiowania tablicy, jak w poniższych przykładach:

>> M = [ 2 5; 6 10]  % wyrazy w wierszu rozdziela spacja lub przecinek, po średniku lub

% przejściu do nowej linii rozpoczyna się następny wiersz

M =

2     5
6    10

>> M = [5 10 15 20]  % definicja 4-elementowego wektora wierszowego (poziomego)

M =

5   10   15   20

>> M = [10; 11; 12]  % definicja 3-elementowego wektora kolumnowego (pionowego)

M =

10
11
12

M = [min: krok: max]  % tworzy wektor o elementach od min do max, z krokiem postępu

% krok, pominięcie parametru krok domyślnie oznacza krok = 1

>> M = [0: 1: 10]

M =

0  1  2  3  4  5  6  7  8  9  10 

>> M = [1: 5; 10: 10: 50]

M =

 1    2    3    4    5 
10  20  30  40  50 

Wywołanie funcji linspace(a, b, n) tworzy wektor o n elementach równomiernie (liniowo) rozmieszczonych, poczynając od pierwszego a, kończąc na ostatnim b,

>> M = linespace(0, 20, 3)

M =

 0   10   20

Nic nie stoi na przeszkodzie aby nową macierz (tablicę) ”poskładać” z już istniejących, trzeba przy tym pamiętać o zgodności wymiarów, np.


>> M1 = [1 2 3; 4 5 6];

>> M2 = [7 8; 9 10];

>> M3 = [0 0 0 0 0; 0 0 0 0 0];

>> M = [M1 M2; M3]

M =

 1    2    3    7    8
 4    5    6    9   10
 0    0    0    0    0
 0    0    0    0    0

MATLAB wykonuje operacje rachunku macierzowego; dodawanie, odejmowanie, na szerszy komentarz zasługuje mnożenie. Przy dwóch danych tablicach A i B, z których A ma tyle kolumn co B wierszy, zapis

>> C = A * B

oznacza mnożenie zgodnie z regułą rachunku macierzowego (mnożenie Cauchy’ego), tj. element Cij macierzy C jest sumą Cij = Aik Bkj ( sumowanie po powtarzających się wskaźnikach). Zapis C = A^2 jest równoważny C = A * A, z kolei C = A/B jest równoważny C = A B-1 (B-1 jest macierzą odwrotną do B).

Druga opcja to tzw. mnożenie tablicowe. Komenda A.*B ( mnożenie z kropką) powoduje mnożenie każdego elementu tablicy A przez jego odpowiednik w tablicy B. Tak mnożone tablice muszą mieć ten sam wymiar, przykład

>> A = [1 2; 3 4];

>> B = [5 6; 7 8];

>> C = A. * B

C =

 5   12
21   32

Dzielenie i potęgowanie tablicowe zapisujemy A./B, A.^2, odpowiednio. W celu wywołania pojedynczego elementu macierzy, np. Cij piszemy C(i, j). Pierwszy wskaźnik to numer wiersza, drugi kolumny, przy czym wiersze i kolumny numerowane są od jedynki. MATLAB dysponuje funkcjami, wykonującymi typowe operacje na macierzach, zwłaszcza kwadratowych.

MATLAB posiada funkcje, umożliwiające szybkie tworzenie typowych, często używanych macierzy oraz niektórych tzw. macierzy specjalnych, oto kilka z nich:

zeros() - macierz składająca się z samych zer

eye() - macierz jednostkowa, jedynki na głównej przekątnej, pozostale elementy zerowe

ones() - macierz składająca się z samych jedynek

Parametrem tych funkcji może być liczba naturalna np. n, wówczas powstaje macierz kwadratowa o wymiarach n * n, lub para liczb naturalnych np. (n, m) co prowadzi do utworzenia macierzy o wymiarach n * m.

hadamard() - macierz kwadratowa, elementami są wartości ±1 (tzn. +1 lub -1)

Kolumny macierzy Hadamarda są parami ortogonalne (zob. help hadamard).

magic() - macierz kwadratowa o wymiarze >2, każdy składnik macierzy jest liczbą naturalną (całkowitą dodatnią) różną od wszystkich pozostałych. Taka macierz jest też zwana kwadratem magicznym. Suma liczb w każdej kolumnie, na każdej przekątnej i w każdym wierszu jest taka sama, wartość największego elementu takiej macierzy o wymiarze n wynosi n2 - doprawdy magiczna macierz.

Parametrem magic() może być tylko jedna liczba naturalna np. n ( n > 2), wtedy macierz ma wymiar n * n. Zobaczmy jak wygląda najmniejsza macierz magiczna ( n = 3)

>> M = magic(3)

M =

 8    1    6
 3    5    7
 4    9    2


W celu zwiększenia szybkości pracy programu wykonywanego w środowisku MATLAB, można przeprowadzać tzw. wektoryzację algorytmów. Dotyczy to zwłaszcza algorytmów intensywnie używających instrukcji sterujących (np. pętli for). Instrukcje MATLAB wpisane przez użytkownika są odczytywane i przetwarzane potokowo, za pośrednictwem interpretera kodu. Wprowadza to pewne "organiczne" opóźnienie wykonywania w porównaniu do takich języków jak asembler czy C/C++. MATLAB jest silnie zoptymalizowany pod kątem wykonywania działań na macierzach. Wektoryzacja polega zatem na zastępowaniu klasycznych instrukcji sterujących (pętli) przez takie specyficzne operacje na macierzach czy wektorach, które wykonywane/interpretowane są bardzo szybko w środowisku MATLAB. Posłużę się bardzo prostym przykładem - utworzenie wektora, którego składowymi są wartości funkcji sinus. Przy okazji wspomnę o parze współpracujących ze sobą matlabowskich funkcji tic() i toc(), które w duecie działają jak stoper do pomiaru czasu wykonywania bloku instrukcji w programie. Funkcja tic() włącza stoper, toc() go zatrzymuje i wypisuje na ekran czas (elapsed time) jaki upłynął od ostatniego wywołania tic().

x = 0.0;
disp('Bez wektoryzacji:')
tic
for m = 1:10001
y(m) = sin(x); 
x = x + 1; 
end
toc

Po wykonaniu programu (środowisko Octave 3.2.3 gcc 4.4.0 ) na ekranie mojego komputera zobaczyłem komunikat

Elapsed time is 0.42 seconds.

Wektoryzacja tego samego programu wygląda następująco,

disp('Po vectoryzacji:')
tic
x = 0:10000;
y = sin(x);
toc

Tym razem pojawił się napis,

Elapsed time is 0.0039 seconds.

Kod wektoryzowany został wykonany ponad sto (!) razy szybciej, chyba warto wektoryzować, prawda ? Takie "przebicie" w transakcjach ekonomicznych, rynkowych to już raczej marzenie inwestora/menadżera. Zobaczmy jeszcze drugi przyklad, mnożenie macierzowe i tablicowe dwóch macierzy

N = 50;              / wykonaj oblicznia dla różnych N, np N = 3, 10, 20, 50
A = ones(N);
B = zeros(N);      / zamień sobie zeros(N) na magic(N), eye(N) lub rand(N)
disp('Bez wektoryzacji')
tic
for n = 1:N; 
for m = 1:N; 
Tablica(n,m) = A(n,m)*B(n,m);                  %  mnożenie tablicowe
Macierz(n, m) = 0.0;
for p = 1:N
Macierz(n,m) = Macierz(n,m)+ A(n,p)*B(p,m);            %  mnożenie macierzowe
end 
end 
end
toc
Tablica;                  
Macierz;
 
disp('Po wektoryzacji')
tic
Tablica_wekt = A.*B;
Macierz_wekt = A*B;
toc
Tablica_wekt;
Macierz_wekt;


Różnica czasów wykonania powyższych operacji bez- oraz z wektoryzacją jest także imponująca. Jeżeli szybkość działania przygotowywanego przez Ciebie programu w środowisku MATLAB ma istotne, wektoryzacja powinna pomóc. Naturalnie, o ile wykorzystywane w programie algorytmy poddają się wektoryzacji. Możesz także przyspieszyć działanie pętli for poprzez wstępną inicjalizację (prealokację) macierzy lub wektorów, w których będą skladowane obliczenia wykonywane w "długich" pętlach. Zobacz kod poniżej, ten kod używa funkcji ones(), aby dokonać prealokacji wektora tworzonego w pętli for. Dzięki temu pętla będzie się niekiedy wykonywać szybciej, zobaczmy przyklad

disp(' Bez prealokacji wektora ')
clear all
tic;
for k = 1:10000;
w(k) = cos(k) ;
end
toc;

Bez prealokacji, interpreter MATLAB-a tworzy i inicjalizuje (nadaje wartość) kolejną skladową wektora w przy każdym przebiegu pętli. Nadajmy wstępne wartości (prealokacja) wszystkim składowym wektora w "hurtowo", jeszcze przed jego użyciem w pętli,

disp(' Z prealokacja wektora ')
clear all
w = ones(1,10000);       % alokacja wektora w 
tic;
for k = 1:10000;
w(k) =  cos(k) ;
end
toc;

Do wstępnej inicjalizacji wektorów lub macierzy, rutynowo wykorzystuje się funkcje zeros() lub ones(), z uwagi na wygodę. Prealokacja niekiedy zwiększa szybkość wykonywania pętli, ale na ogół już nie w tak spektakularnym, raczej znacznie mniejszym stopniu, niż wektoryzacja.

Instrukcje sterujące

W MATLAB nie występują zmienne typu logicznego. W zamian obowiązuje konwencja; jeżeli wyrażenie lub zmienna ma liczbową wartość zero, to jego logiczna wartość (stan) jest false. Gdy wartość liczbowa wyrażenia różni się od zera, wówczas jego stan logiczny jest true. Naturalnie, stan logiczny true albo false otrzymamy również, gdy w danym wyrażeniu występują operatory logiczne (np. > , < , = =, ...). Omawiam niżej grupę instrukcji w MATLAB, które służą do podejmowania decyzji w programie, są to tzw. instrukcje sterujące.

Instrukcja warunkowa if

Wymienię warianty instrukcji warunkowej if:

a)

if warunek
instrukcja
end
Przykład,
if ( z < 0 )
disp(‘ z ujemna’)
end

b)

if warunek
instrukcja_1
else
instrukcja_2
end
Przykład,
if  z  < 0
disp(‘ Uwaga ! liczba ujemna’)
else
disp(‘OK')
end

Konstrukcja ze słowem kluczowym elseif umożliwia kolejne sprawdzanie kilku warunków. Warunek występujący przy elseif jest sprawdzany tylko wtedy, gdy warunki pojawiające się wcześniej nie są prawdziwe:

c)

if warunek_1
instrukcja_1
elseif warunek_A
instrukcja_A
elseif warunek_B
instrukcja_B
elseif warunek_Z_
instrukcja_Z
else
instrukcja_2
end

Przykład:

if      z < 0
disp(‘z mniejsza od zera ’) 
elseif  z > 0
disp(‘z większa od zera’)
else
disp(‘Uwaga, zmienna z ma wartość zero !!’)
end

Pętla while

while
warunek
instrukcja
end

Jeżeli warunek jest prawdziwy, to instrukcja zostaje wykonywana i ponownie przystępuje się do testowania prawdziwości warunku. Instrukcja wykonywana jest tak długo, jak długo jest prawdziwy warunek. Przykład,

k  =  10;
while k  >  0
k  =  k - 1;
disp(k)
end

Na ekranie pojawi sie kolumna liczb poczynając od 9 aż do 0.

Pętla for

Pętla for umożliwia wykonanie instrukcji lub sekwencji instrukcji określoną ilość razy. Schemat pętli

for zmienna_iterowana = wyrażenie
instrukcja
end

Zamiast jednej instrukcji, może być ich wiele, oddzielonych np. przecinkami. Wyrażenie występujące w pętli for ma często typową postać: indeks = pocz: krok: koniec. Przykład, obliczanie wartości średniej losowo wybranych n liczb z przedziału <0, 1>

n = 10;
los  = rand(1,n);
suma = 0;
for k = 1 : n
suma =suma + los(k);
end
srednia = suma/n

>> srednia = 0.62738  % ta wartość zależy od wylosowanych wartości wektora los.

Następny przykład,

for k  = 1 : 2
for n = 1 : 2 
wynik =  k + n;
disp( wynik)
end
end

Na ekranie pojawi się 4 – elementowa kolumna

2
3
3
4

Instrukcja wyboru switch

Umożliwia wybieranie jednej z wielu możliwych ścieżek postępowania. Mówiąc w języku kodu programowania, umożliwia wybieranie jednej instrukcji (lub bloku instrukcji) z szerszej palety. Składnia instrukcji wyboru switch:

switch WYBOR  % tylko skalar lub łańcuch
case wyraz_1
instrukcja  % lub blok instrukcji oddzielonych przecinkami
case wyraz_2
instrukcja  % lub blok instrukcji oddzielonych przecinkami
...................
case wyraz_n
instrukcja
otherwise  % nie musi koniecznie wystąpić
instrukcja  % lub blok instrukcji, nie musi wystąpić
end

Porównuje się wyrażenie na wejściu switch (nazwałem je WYBOR) po kolei z wyrażeniami występującymi przy słowach kluczowych case. Jeżeli np. WYBOR = = wyraz_2, to będzie wykonana instrukcja (blok) zapisana poniżej drugiego case, pozostałe instrukcje będą pominięte. Jeżeli wyrażenie WYBOR nie zgadza się z żadną z wartości podanych przy etykietach case, to wykonywane są instrukcje po otherwise. Uwaga, przy etykiecie case można w nawiasach klamrowych wpisać grupę wyrażeń, np. case {wyraz1, wyraz2, wyraz3, ...}. Instrukcja (blok instrukcji) występująca po tej etykiecie będzie wykonana, gdy WYBOR zgadza się z którymkolwiek z wyrażeń w nawiasie klamrowym. Przykład,

switch (sign(x))
case -1
disp(‘ujemna’)
case 1
disp(‘dodatnia’)
otherwise
disp(‘zero’)
end

Następny krótki przykład:

switch k
case 1
n =  10
case 2
n =  50
case 3
n = 100
otherwise
disp(’dopuszczalne tylko k = 1, 2 lub 3)
end

Popatrz jeszcze na poniższy programik o charakterze dydaktycznym, wspomagający nauczanie mnożenia liczb naturalnych (tabliczka mnożenia "do stu"). Innymi słowy, instrukcja wyboru switch na usługach naszych milusińskich.

disp('ILOCZYN LICZB NATURALNYCH - SAMOUCZEK ' )
disp(' ')
a = 't' ;
 
while a == 't' | a == 'T'
 
m = floor(10*rand) + 1;
n = floor(10*rand) + 1;
iloczyn = m * n;
disp(['Ile wynosi iloczyn: ' num2str(m) ' * ' num2str(n) '        '])
disp('  ')
odpowiedz = input( '  ');
if odpowiedz == iloczyn
los = 1 + fix(6 * rand());
disp('  ')
switch los 
case 1
disp(' Wspaniale ')
case 2
disp(' Doskonale sobie radzisz ' )
case 3
disp(' Bardzo dobrze, zadziwiasz mnie ')
case 4
disp(' Doskonale ')
case 5
disp(' Gratulacje, no, no ')
case 6
disp(' Bardzo fajnie, tak trzymaj ')
end
 
else
los = 1 + fix(6 * rand());
disp('  ')
 
switch los 
case 1
disp(' Kiepsko ')
case 2
disp(' No niestety, .... ' )
case 3
disp(' Niedobrze ')
case 4
disp(' Przykro mi, ale musisz jeszcze popracowac ')
case 5
disp(' Wynik niepoprawny ')
case 6
disp(' Nie, nie, policz jeszcze raz ')
end
end
 
disp('  ')
a = input('Jeszcze raz, T/N  ? ', 's');
if a == 'n' | a == 'N'
disp('    ')
disp('Bye, bye, bye, ...')
disp('  ')
end
end

Instrukcja break

Przerywa instrukcję switch, powoduje także natychmiastowe zakończenie/przerwanie wykonywania pętli; for, while, ... .

Jeżeli mamy do czynienia z kilkoma pętlami, zagnieżdżonymi jedna wewnątrz drugiej, to instrukcja break powoduje przerwanie tylko tej najbardziej wewnętrznej pętli, w której występuje. To przerwanie pętli z wyjściem o jeden poziom wyżej.

Instrukcja continue

Instrukcja continue występuje zwykle wewnątrz pętli while, for, ... . Powoduje ona zaniechanie wykonywania instrukcji będących treścią pętli, jednak sama pętla nie zostaje przerwana. Instrukcja continue przerywa tylko bieżący obieg pętli i zaczyna następny obieg, kontynuując pracę pętli. Innymi słowy, napotkanie instrukcji continue odpowiada natychmiastowemu przejściu/skokowi na sam koniec pętli, komputer "pomyśli", że wykonał treść pętli i przystąpi do następnego obiegu. W ten sposób instrukcje pętli leżące poniżej continue, aż do końca pętli, nie będą realizowane.

Funkcje

MATLAB oferuje nie tylko szeroką paletę funkcji standardowych, ale umożliwia tworzenie przez programistę jego własnych funkcji. Definicję takiej funkcji należy koniecznie umieścić w pliku o rozszerzeniu *.m, podobnie jak w przypadku skryptu. Różnica polega na tym, iż nazwa *.m pliku skryptowego może być dowolna, natomiast nazwa *.m pliku zawierającego funkcję definiowaną samodzielnie przez programistę musi by taka sama, jak nazwa funkcji. Schemat definicji funkcji wygląda następująco,

function [zmienne wyjściowe] = nazwa_funkcji(parametry wejściowe);

...........................  % instrukcje do wykonania

end

Definicja własnej funkcji musi zaczynać się od słowa kluczowego function i dalej zawiera:

- nazwę zmiennej zwracającej wartość obliczoną przez funkcję. Jeżeli takich zmiennych jest więcej, ich listę umieszcza się w nawiasie kwadratowym [ ]
- znak podstawienia ( = )
- wybraną przez programistę nazwę funkcji z umieszczoną w nawiasach półokrągłych listą parametrów wejściowych, od których funkcja zależy
- dalej zaleca się zamieszczenie komentarza informującego, co ta funkcja robi, jakie znaczenie mają parametry wejściowe, ..itp.,
- potem jest ciąg instrukcji do wykonywania przez funkcję czyli tzw. ciało funkcji.
- definicję funkcji kończy słowo kluczowe end.

Treść komentarza, który jest zalecany ale nie obowiążkowy, zostanie na nasze życzenie wyświetlona na ekranie komputera, wystarczy w tym celu wydać polecenie

>> help nazwa_funkcji

Przykład bardzo prostej funkcji, deklarującej dwa argumenty wejściowe (nazwane x oraz y) i jedną zmienną wyjściową o nazwie wynik,

function  wynik   =  moja(x, y)
% oblicza kombinację funkcji trygonometrycznych 
wynik = sin(x) + cos(y)sin(y) * cos (x);
end

Po zapisaniu tego kodu w pliku moja.m, korzystamy z tej funkcji wywołując ją z linii komend

z  =  moja(0, pi)
z   =  - 1

Po wywołaniu funkcji, wykonywane są instrukcje będące w jej ciele, po ich wykonaniu funkcja zwraca wartości parametrów wyjściowych (jeżeli występują), kończy działanie i zwraca sterowanie do tego miejsca m-skryptu (linii komend), z ktorego zastała wywołana. Dokładniej, do linijki kodu następnej po instrukcji wywołania funkcji, która właśnie zakończyła pracę. Jednakże, jeżeli w ciele funkcji pojawi się słowo kluczowe return , wówczas w tym miejscu funkcja kończy działanie i zwraca sterowanie. Instrukje ciała funkcji poniżej return nie będą wówczas wykonywane. Tak się często dzieje, gdy w ciele funkcji pojawiają się rozmaite warunki kogiczne, w zależności od ich spełnienia lub nie, dalsza część instrukji funcjie jest wykonywana lub nie.przerywane jest wykonywanie instrukcji. Wszystkie zmienne kreowane w ciele funkcji są lokalne tzn. niewidoczne spoza ciała funkcji, spoza *.m pliku w którym jest definicja tej funkcjiJeżeli chcemy, aby ta sama zmienna/macierz mogła być używana przez kilka różnych funkcji, była "widoczna" w całej przestrzeni roboczej, należy ją poprzedzić słowem global.

W MATLAB spotykamy rozmaite funkcji, ale niekiedy wywoływane z jednym, innym razem z dwoma, trzema lub większą liczbą parametrów. Odpowiednio, liczba zmiennych wyjściowych, zwracanych przez funkcję, bywa także różna. Przypomnij sobie, iż poprawne jest wywołanie np. rand(10) jak i rand(10, 30). Nie tylko wbudowane, matlabowskie funkcje mają taką cechę, dość często spotykaną praktyką przy programowaniu w MATLAB jest konstruowanie przez użytkownika/programistę swoich własnych, które też tak się będą zachowywać. W tym kontekście należy wspomnieć o użytecznych funkcjach nargin() (number of input arguments) i nargout() (number of output arguments). Umieszczone we wnętrzu (ciele) jakiejś funkcji zwracają, odpowiednio, ilość argumentów wejściowych i wyjściowych tej funkcji, w której są zanurzone. Korzystając z tego, za pośrednictwem instrukcji warunkowych można np. tak zaprojektować swoją funkcję, aby wykonywała ona różne zadania (warianty) w zależności od liczby parametrów wejściowych, z którymi funkcja jest wywoływana. Wywołane spoza ciała funkcji, np. z linii komend, nargin()/nargout() podają liczbę zadeklarowanych parametrów wejściowych/wyjściowych tej funkcji, której nazwę podamy jako ich argument. Załóżmy, że w kodzie programu mamy zdefiniowaną funkcję nazwaną "moja_funkcja" i jej deklaracja zawiera trzy parametry wejściowe, wówczas

>> nargin(moja_funkcja)

>> ans = 3

Oto typowy przykład zastosowania nargin

function wykres(p1, p2)
x = inspace(0, 30);
if nargin == 1
disp(' Podano jeden parametr: p1')
disp(' Wykres f(x) = x^(p1) ')
y = x.^(p1);
figure
plot(x, y);
elseif nargin == 2
disp(' Podano dwa parametry; p1 oraz p2 ');
disp(' Wykres f(x) = x.^(p1.*p2) ')
y = x.^(p1.*p2);
figure
plot(x, y)
end

Uruchom kilka razy tą funkcję podając jej zarówno jeden, jak i dwa parametry, zobacz otrzymane wykresy, sprawdź, czy wszystko się zgadza.

W *.m pliku zawierającym definicję funkcji można poniżej umieścić definicje następnych funkcji. Wówczas ta pierwsza funkcja, która nazywa sie tak samo jak *.m plik (dokładniej, jak jego sygnatura tj. nazwa bez rozszerzenia), jest funkcją główną/nadrzędną/pierwszoplanową (primary) w przestrzeni roboczej tego *.m pliku. Inne funkcje zdefiniowane w tym pliku są nazywane w języku angielskim subfunction. Spotykane określenie podfunkcja jest dość automatycznym tłumaczeniem, może nazwa "podrzędna" byłaby alternatywą... ?. Najpierw w *.m pliku należy umieścić definicję funkcji głównej, poniżej sukcesywnie definicje funkcji podrzędnych (podfunkcji). Kolejność funkcji podrzędnych nie ma znaczenia, ale główna musi być na początku. Funkcja nadrzędna i wszystkie inne funkcje (subfunctions) z jej pliku, jej przestrzeni roboczej, widzą się wzajemnie. Inaczej mówiąc, np. z ciała funkcji nadrzędnej mogę wywołać inną z tych funkcji. Zwracam jednak uwagę, iż zmienne, wprowadzone w ciele jednej z tych funkcji, są lokalne w tej funkcji, tj. nie są widoczne, nie można ich używać w innej funkcji z tego samego pliku. Jeżeli chcemy aby były widoczne, należy postąpić standardowo czyli zadeklarować tą zmienną jako globalną. Naturalnie, poza przestrzenia roboczą tego *.m pliku (np. z linii komend), widoczna jest tylko funkcja nazywająca się tak, jak ten plik. Funkcje podrzędne (subfunctions) tego pliku nie są widoczne z linii komend lub z innego m-pliku. Funkcje nargin oraz nargout występujące w kodzie funkcji podrzędnej zwracają dane dotyczące danej funkcji, a nie funkcji nadrzędnej/głównej. Schemat *.m pliku z kilkoma funkcjami,

function [ ] = nadrzedna( );

…............  % instrukcje funkcji nadrzedna

end

function [ ] = sub_1( )

…...........  % instrukcje funkcji sub_1

end

function [ ] = sub_2( )

.............  % instrukcje funkcji sub_2

end

W schemacie powyżej wpisałem tylko dwie podfunkcje, możesz zdefinować tyle pofunkcji, ile potrzebujesz. Przykład,

function [w1, w2, w3] = f(p1, p2, p3)           % początek funkcji głównej/nadrzędnej
% tylko skalarne parametry wejściowe
 
if (nargin < 1)
disp('Podaj przynajmniej jeden argument funkcji !')
return;
endif
 
if (nargin > 3)
disp('Przepraszam, maksymalnie trzy argumenty !')
return;
endif
 
d = length(p1(:))       % liczymy, ile elementów ma p1
if d ~= 1                   
disp('Tylko skalarne parametry wejściowe')
return;
endif
 
global x;
x = 0:0.1:2.*pi;
w1 = p1;
y = sin(w1.*x);
disp('Jestem w f(), wykonuje pierwszy wykres')
disp(' ')
y = sin(w1.*x);
plot(x, y)
pause(8)
close()
if nargin == 1
disp('Bye, bye, ...')
return;
endif
 
d = length(p2(:))
if d ~= 1
disp('Drugi parametr wejsciowy nie jest skalarem ')
return;
endif
 
w2 = w1.*p2;
disp('Uruchomiam a()');
disp(' ')
a(w2)                            % wywołanie podfunkcji a() z ciała funkcji f()
pause(8)
close()
if nargin == 2
disp(' Bye, bye, ...')
return;
endif
 
d = length(p3(:))
if d ~= 1
disp('Trzeci parametr wejsciowy nie jest skalarem ')
return;
endif
w3 = w2.*p3;
disp(' Uruchamiam b() ');
disp(' ')
b(w3)                              % wywolanie pofunkcji b() 
pause(8)
close()
disp(' Bye, bye, ...' )
return;
end                                 % koniec funkcji głównej/nadrzędnej
 
function a = a(a1)              % początek podfunkcji a()
global x;
wa = a1;
ya = cos(wa.*x);
disp(' Jestem w a(), wykonuje drugi wykres' );
plot(x, ya)
pause(8)
close()
end                                 % koniec podfunkcji a()
 
function b = b(b1)              % początek pofunkcji b()
global x;
wb = b1;
yb = cos(wb.*x).*x.^3;
disp(' Jestem w b(), wykonuje trzeci wykres ');
plot(x, yb)
pause(8)
close()
end                                  % koniec podfunkcji b()


W plikach *.m umieszcza się zwykle bardziej rozbudowane funkcje. Istnieje możliwość definiowania własnych funkcji bezpośrednio w linii komend, bez zapisywania do *.m pliku – nazywamy ją wówczas funkcją inline. W praktyce tym sposobem definiuje się funkcje o niewielkim kodzie. W przykładzie

>> f = inline( ‘4*x + 7’ , ‘x’ )

tworzę funkcję f(x) = 4 x + 7 i od tego momentu jest ona do dyspozycji

>> f(2)

>> ans = 15

Równolegle do funkcji inline istnieje możliwość definiowania tzw funkcji anonimowej - jest to wygodne, gdy np. kod tworzonej funkcji nie jest rozległy, pracujemy interaktywnie w otwartej sesji Matlab/Octave, nie chcemy tworzyć m-pliku, ...itp. Przy definicji funkcji anonimowej używa się znaku @ , dobrze znanego Ci elementu adresu poczty elektronicznej. Przejdę od razu do przykładu, który unaoczni sposób tworzenia takiej "podręcznej" funkcji

>> potega = @(x) x.^6

W nawiasie na prawo od @ umieszczamy zmienną (lub oddzielone przecinkami zmienne) funkcji, dalej jest jej treść/ciało funkcji. Operator @ tworzy uchwyt do zbioru instrukcji występujących po jego prawej stronie, w moim przykladzie uchwytowi przypisałem nazwę "potega".

Wywołując funkcję anonimową dostajemy

>> potega(2)

>> ans = 64

Lista parametrów funkcji anonimowej może być pusta, jak w przykładzie niżej

>> data = @() datestr(clock);

gdzie funkcja anonimowa jest skonstruowana z dwóch wbudowanych funkcji MATLABA (zob. help clock, help datestr). Zwracam uwagę, iż tą funkję należy wywołać z pustym nawiasem, wywołanie tylko nazwy uchwytu, bez pustego nawiasu, nie spowoduje uruchomienia funkcji. Wywołajmy tą funkcję,

>> data()

>> ans = 04-Feb-2011 15:16:31  % Ty zobaczysz inną datę, datę Twojego kliknięcia

Tworzę teraz funkcję anonimową, którą następnie całkuję numerycznie w przedziale <0, 1> przy pomocy matlabowskiej funkcji quad() (zob. help quad)

a = 1; b = 1; c = 1; d =1;  
wielomian = @(x) (a*x.^3 + b*x.^2 + c*x + d);
quad(wielomian, 0, 1);
ans = 2.0833

Całka nieoznaczona z funkcji f(x) = ax3 + bx2 + cx + d wynosi ax4/4 + bx3/3 + cx2/2 + dx , obliczona dla a = b = c = d = 1 w granicach od 0 do 1 daje wartość 2 + 1/12 = 2.0833.

Grafika

MATLAB umożliwia tworzenie grafiki, rozmaitych wykresów, prezentację przebiegów funkcji,...itp. Bardzo szczegółowe omówienie tych zagadnień wykraczałoby znacznie poza ramy tego krótkiego wprowadzenia do MATLAB. Podam tutaj esencjonalne informacje dotyczące zwłaszcza tworzenia wykresów funkcji wbudowanych lub definowanych przez programistę, prezentacji zgromadznych danych (eksperymentalnych, giełdowych, ankietowych, ...itp.).

Wykresy standardowe

Jednym z najprostszych jest wykres słupkowy uruchmiany przy pomocy funkcji bar(), jak w przykładzie

x = 0: 0.1: 2*pi;
y = sin(x);
bar(x, y)

Oto rezultat,

Bar.png

Podanie funkcji bar() argumentu wektorowego (macierz o jednym wierszu, w naszym przykładzie wektora y), powoduje wykreślenie słupków reprezentujących wartości tego wektora w funkcji jego indeksów. Inaczej mówiąc, na osi x bedą indeksy wektora, od 1 do n, gdzie n jest liczbą elementów wektora. Gdy argumentem jest macierz o wymiarach n * m, wówczas bar() tworzy m grup (tyle, ile kolumn macierzy). W każdej z tych grup jest n słupków, rerezentujących wartości elementów danego wiersza macierzy. Dla ilustracji utwórzę najpierw trójwymiarową macierz,

>> M = magic(3)

M =

 8    1    6
 3    5    7
 4    9    2

i następniej jej wykres słupkowy,

bar(M)

Bar1.png

Funkcja barh() działa podobnie jak bar(), ale rysuje w poziomie,

x = -2.8:0.1:2.8;
y = exp(-x.*x);
barh(x, y)
grid on                   % tworzy siatkę na wykresie

Barh1.png

Jak prawie każda funkcja wbudowana, bar() posiada wiele opcji, np. można manipulować szerokością rysowanego słupką, podając funkcji bar() trzeci parametr wejściowy. Polecam komendy: help bar oraz help bar3 ("słupki" trójwymiarowe) w celu poznania pełnej palety opcji. Ta sama uwaga dotyczy innych funkcji graficznych pokrótce omawianych poniżej, nie będę jej za każdym razem powtarzał.

Wykres pałeczkowy (ang. stem plot) pokazuje wszystkie rzędne badanej funkcji,

x = 0:0.1:pi;  
y = sin(x.*x).*exp(-x)
stem(x, y)
grid on

Stem.png

Kolejna funkcja graficzna to errorbar() (zob. help errorbar), umożliwia szybkie narysowanie wykresu z zaznaczeniem błędów (pomiarowych, statystycznych, ...itp). Ograniczę się do bodaj najczęściej wykorzystywanego wariantu errorbar() - wykres danej funkcji jednej zmiennej z zaznaczonymi błędami wartości funkcji. Wówczas parametrami wejściowymi są trzy wektory, jak w przykładzie poniżej

x = 0:0.1:2*pi;                 % wektor odciętych
y = sin(x).*cos(x);            % wektor rzędnych (wartości) wykreślanej funkcji y(x) = sin(x)cos(x)
e = rand(size(x))/4;           % wektor błędów wartości y(x), tutaj jego elementy generowane przypadkowo
errorbar(x, y, e);               % parametry wejściowe (trzy wektory)

Error.png

Przy analizie rozmaitych zbiorów danych liczbowych zachodzi potrzeba unaocznienia rozkładu wartości tego zbioru, czyli jego histogramu. MATLAB dysponuje odpowiednią funkcją hist() (zob. help hist ).Jeżeli podamy hist() tylko jeden parametr, który jest wektorem a wartości jego elementów leżą w przedziale < min, max >, to ten przedział zostanie podzielony na 10 równych "'plastrów"/kontenerów. Dalej hist() policzy, ile elementów wektora ma wartości leżące w pierwszym, drugim, ... ostatnim kontenerze i zwróci te liczby elementów w poszczególnych kontenerach jako wektor wyjściowy, który wobec tego będzie zawierał tylko liczby naturalne (10 liczb). Dla przykładu utworzę macierz a,

>> a = [1 2 3 4 5 6 7 8 9 10];

Zawiera ona elementy o wartościach z przedziału <1, 10>. Każdy element/składnik macierzy a ma inną wartość. Spodziewamy się, iż każda wartość "wpadnie" do innego plastra/kontenera, i będzie w danym kontenerze osamotniona. Wobec tego funkcja hist() powinna zwrócić tylko jedynki. Sprawdźmy to,

>> n = hist(a)

>> n =

1   1   1   1   1   1   1  1  1   1

Naturalnie, zechcesz przetestować konstrukję histogramu także na innych, mniej regularnych i bardziej licznych zbiorach danych. Wówczas ich "ręczne" wpisywanie, jak w przykładzie powyżej, nie będzie dobrym pomysłem. Oto wersja poprzedniego kodu nieco bardziej rozwinięta i wygodniejsza do wielokrotnego testowania - losowanie n liczb przyjmujących wartości całkowite dodatnie z określonego przedziału i obliczanie histogramu tego rozkładu.

a = 1;           
b = 11;,        
n = 10;                                     % liczba losowań
for k = 1:n;
dane(k) = floor((b-a)*rand +a);    % losowanie liczby naturalnej z przedziału < a, b - 1>
end
disp(' Wektor danych')
dane
disp(' Histogram danych')
h = hist(dane)

Wynik pracy mojego komputera,

>> dane =

4   2   8   1   6   5   7  4  6   10

W zbiorze dane widzimy: jedynkę, dwójkę, dwie czwórki, jedną piątkę, dwie szóstki, po jednej siódemce, ósemce i dziesiątce, brak trójek i dziewiątek. Bezbłędnie informuje o tym histogram

>> h =

1   1   0   2   1   2   1  1  0   1

Wykonaj kilka obliczeń, zmieniając n np. n = 50, 100, 500, 1000 itd. oraz ewentualnie liczby a i b. Licząc dla większych wartości n zrezygnuj z wypisywania na ekran wektora dane, tj. dane zamień na dane;. Dzieląc wszystkie elementy histogramu (wektora h z przykładu) przez całkowitą liczbę elementów wektora danych (dane) dostajemy zbiór liczb, których suma wynosi 1. Liczby te obrazują prawdopodobieństwa pojawienia się wartości z danego kontenera/plastra, wśród wszystkich wartości wektora dane.

Wykres histogramu dostaniesz "od ręki" wypisując

a = 1;, b = 11;, n = 10;      
for k = 1:n;
dane(k) = floor((b-a)*rand +a);    
end
disp(' Wektor danych')
dane
disp(' Wykres histogramu')
hist(dane)                        %  zamiast instrukcji podstawienia h = hist(dane) z poprzedniego przykładu

Hist.png

Na osi poziomej zaznaczone są kolejne kontenery, na pionowej liczby elementów których wartości mieszczą się w danym kontenerze/plastrze. Domyślna liczba kontenerów wynosi 10, można ją samodzielnie ustalić podając funkcji hist() drugi, skalarny parametr. Podkreślam, iż ten drugi parametr musi być skalarem, aby hist() zinterpretowała go jako wymaganą liczbę kontenerów. Poniżej losowanie liczb z przedziału <0, 1> i histogram o dużej rozdzielczości (w = 300, tj. odcinek < 0, 1> podzielony na 300 części).

n = 900;
dane = rand(1, n);        %  ''n'' liczb losowych z przedziału  < 0, 1 >
hist(dane, 300)

Hist2.png

Drugi parametr hist() może być wektorem, nazwę go x. Wówczas będzie tyle kontenerów, ile elementów wektora x czyli nk = length(x) kontenerów. Pokrywać będą przedział od najmniejszej do największej wartości wśród jego elementów. Typowy przyklad,

dane = rand(100, 1);       
xp = min(dane);
xk = max(dane);
krok = 0.1;
x = xp: krok: xk
hist(dane, x)

Gdy pierwszym argumentem hist() jest macierz o wymiarach m*p, hist() traktuje kolumny tej macierzy jako wektory i zwróci histogram/macierz o rozmiarach 10*p. Dalsze przykłady użycia hist() znajdziesz np. w skrypcie * MKZR:Dodatek#Histogramy .

Na posiedzeniach rad nadzorczych, zarządów korporacji, kierownictw przedsiębiorstw, spółek, ...itp. bardzo często spotykaną formą prezentacji rozmaitych danych ekonomicznych jest tzw. wykres kołowy/procentowy (ang. pie chart). Pokazuje on, jak podzielony jest przysłowiowy tort (pie - placek owocowy). Ile procent przychodu firmy poszło np. na: podatki, ubezpieczenia, inwestycje, materiały do produkcji, media, koszty osobowe, ... no i najbardziej interesujące, jaki procent przychodu stanowi dochód/zysk firmy. Taki wykres można sporządzić np. w arkuszu kalkulacyjnym Excel, teraz pokażę, jak to zrobić w MATLAB (zob. help pie ). W najprostszej wersji, jako argument funkcji pie() wystarczy podać wektor reprezentujący jakieś dane. Funkcja pie() oblicza procentowy udział każdego skladnika w ich sumie. W przykładzie poniżej ze wzgledów dydaktycznych przygotowałem wektor danych tak, aby suma jego składników wynosila 100. Wobec tego wartość każdego składnika to procentowy wkład do sumy wszystkich skladników, sprawdźmy to na wykresie

dane = [16  4  5  5  10  55  5];
pie(dane)


Pie.png

Utworzysz ten sam wykres, ale z opisem składników prezentowanego "tortu", jeżeli w nawiasie klamrowym -po wektorze z danymi - umieścisz tyle 'opisów', ile elementów wetora,

dane = [16  4  5  5  10  55  5];
pie(dane, {'podatek', 'ubezpieczenie', 'materiały', 'media', 'wynagodzenia', 'premie zarządu', 'zysk firmy' } )

Pie2.png

No cóż, teraz wszystko klarownie przedstawione. Jeżeli jesteś członkiem zarządu tej firmy to domniemam, iż przed prezentacją wykresu na zebraniu Rady Nadzorczej, znalazłaś/eś już nową pracę w zupełnie innej branży. Gdyby jednak zależało Ci na powtórzeniu w tej samej firmie tak dobrych dla Ciebie wyników również w przyszłym roku, możesz podjąć techniczną próbę ratowania sytuacji. Pomoże Ci w tym MATLAB, chociaż został stworzony w zasadzie do innych celów. Możesz szczególnie uwypuklić wybrany składnik(składniki) wykresu, tutaj wysoki zysk (!), przy okazji optycznie tonując nastrój i "rozmazując" ewentualne dysproporcje, których nie chcesz eksponować

dane = [16, 4, 5, 5, 10, 55, 5];
eksponuj = [0, 0, 0, 0, 0, 0, 1];
pie(dane, eksponuj, {'PODATEK', 'UBEZPIECZENIE', 'MATERIAŁY', 'MEDIA', 'WYNAGRODZENIA', 'Premie zarządu', 'ZYSK FIRMY' } )
title (' "PROSPERITY" SA, sprawozdanie roczne') 
colormap summer

Pie3.png

Niby to samo, a jednak inaczej, prawda ? Zauważ, iż w celu wyeksponowania danego składnika wykresu kołowego ("tortu"), tworzymy pomocniczy wektor (ja nazwałem go eksponuj ale nazwa dowolna) o tej samej liczbie elementów, co wektor z danymi liczbowymi (dane). Elementami wektora pomocniczego są zera lub jedynki. Jeżeli k-ty element wektora pomocniczego jest jedynką, wówczas k-ty element wektora z danym zostanie wyeksponowany. Wektor pomocniczy należy podać jako drugi argument funkcji pie(), po wektorze z danymi.

Funkcja plot()

Wspomnę o najprostszych opcjach popularnej funkcji plot(),służącej do tworzenia dwuwymiarowych wykresów. Funkcja plot() samoczynnie otwiera nowe okno graficzne, o ile nie istnieje już okno wcześniej otwarte. Jeśli już istnieje otwarte okno graficzne, plot() skorzysta z tego okna i nie otworzy nowego. W każdej chwili możesz otworzyć nowe, aktywne okno graficzne, wpisując komendę

>> figure

Bywa tak, iż masz kilka otwartych okien graficznych, aby jedno z nich uczynić bieżącym (aktywnym), użyj komendy

>> figure(n)

gdzie n jest numerem okna. W MATLAB każde okno graficzne ma swój numer, który jest widoczny na na pasku tytułowym okna. Wyniki działania używanych w programie funkcji graficznych będą pokazywane w tym oknie, wskazanym jako bieżące. Chcesz wyczyścić aktualnie aktywne okno, bez jego zamykania, piszesz clf, clf(n) czyści okno numer n. Komenda close zamyka aktywne okno, close(n) zamyka okno numer n, close all zamyka wszystkie okna.

Komenda hold umożliwia dodawanie krzywych do istniejących już wykresów. Gdy wpiszesz

>> hold on

MATLAB nie usunie istniejącego wykresu, lecz doda do niego nowe dane i jeśli będzie to konieczne, przeskaluje wykres. Komenda hold off przywraca stan naturalny tzn. aktualne okno jest czyszczone przed narysowaniem w nim nowej grafiki. Argumentami funkcji plot() są dwa ciągi liczb, podane np. jako nazwy wektorów, lub wpisane bezpośrednio. Wartości pierwszego wektora są pokazywane na osi X (poziomej), drugiego na osi Y (pionowej). Argumentem może być tylko jeden wektor, na osi X wówczas pojawiają się kolejne indeksy wektora, natomiast jego wartości na osi Y. Dla wektora o wartościach zespolonych np. o nazwie Z, instrukcja plot(Z) jest ekwiwalentna następującej,

>> plot(real(Z),imag(Z)).

Przykład tworzenia wykresu funkcji trygonometrycznej cos(x) w przedziale - 2 pi < x < 2 pi

>> x = -2*pi: 0.1 : 2*pi; >> y = cos(x); >> plot(x, y)


Pm1.jpg

Wyczyść to okno komendą clf, następnie je zamknij poleceniem close. Funkcja subplot(), która posiada trzy parametry (liczby naturalne), pozwala podzielić aktualne okno graficzne na wiele samodzielnych podokien. Wywołanie

>> subplot(n, m, k)

dzieli okno graficzne na m * n podokien (m numeruje wiersze, n kolumny) i ustala jako aktywne podokno o numerze k. Każde podokno ma bowiem swój własny układ współrzędnych i nadany numer wynikający z numerowania podokien wierszami od góry do dołu, natomiast w wierszu od lewej do prawej). Zobaczmy konkretny przykład, kilka niżej wypisanych linii kodu pozwala w jednym oknie przedstawić cztery wykresy różnych funkcji matematycznych.

t = 0:pi/15: 4*pi;
[x,y] = meshgrid(t);
subplot(2,2,1)
plot(sin(t),cos(t)) 
axis equal
subplot(2,2,2)
z = sin(x).*sin(x)-cos(y).*cos(y);
plot(t,z)
axis([0 4*pi -2 2])
subplot(2,2,3)
% z = sin(x).*cos(y);
z = (sin(x).^2).*(cos(y).^2);
plot(t,z)
axis([0 4*pi -1 1])
subplot(2,2,4)
z = (sin(x).^3)-(cos(y).^2);
plot(t,z)
axis([0 4*pi -1 1])
end

Oto wynik interpretacji powyższego kodu przez MATLAB,

Sub3.gif

Wykorzystajmy jeszcze MATLAB do obliczenia wartości czterech funkcji Bessela J(n, x) (n = 0, 1, 2, 3) i spójrzmy na ich wykresy

x = 0:0.1:30;
for k = 1:4,
subplot(2,2,j)
plot(x, besselj(k*ones(size(x)), x))   
end

Bes.gif

Do tej pory, aby nie rozpraszać uwagi pomijałem takie "detale" jak opis osi wykresu, dobór koloru lub znaków "piszących" linie wykresu, ...itp. Oczywiście, jeżeli chcemy zrobić porządny, czytelny wykres/rysunek, musimy go staranie opisać. Do opisu wykresu służy paleta funkcji, których parametrami są najczęściej łańcuchy znaków, najczęściej używane to;

xlabel() – umieszcza łańcuch znaków pod osią x (odciętych)

ylabel() – umieszcza łańcuch znaków wzdłuż osi y (rzędnych)

title() – umieszcza łańcuch znaków centralnie nad wykresem, to tytuł wykresu

legend() – umieszcza legendę (znaczenie różnych symboli lub linii, krótkie info, ... itp.) na wykresie lub obok niego

Funkcja text() służy do umieszczania łańcucha znaków w dowolnym miejscu wykresu. Jeżeli x oraz y są zmiennymi skalarnymi, komenda text(x, y, ’napis’) umieszcza dolny lewy kraniec łańcucha (tutaj napis) w punkcie (x, y) wykresu, przy czym x oraz y są mierzone w jednostkach bieżącego wykresu, zob. help text. Funkcja gtext( ) spełnia analogiczne zadanie jak text( ), jednakże łańcuch będący jej argumentem jest umieszczany na wykresie w miejscu wskazanym myszką, po kliknięciu jej klawisza lub naciśnięciu jakiegokolwiek znaku na klawiaturze (ta funkcja może nie działać, gdy używasz Octave i nowszej wersji Windows). Przykłady i opis wielu opcji funkcji plot() znajdziesz w skrypcie LK http://el.us.edu.pl/ekonofizyka/index.php/Analiza_Szereg%C3%B3w_Czasowych/GNU_Octave#Wizualizacja_danych

Zadania

  1. Napisz program, który oblicza sumę kwadratów liczb naturalnych, od 1 do n.
  2. Dane są dwie macierze M = [3 1; 5 2] oraz N = [1 2; 3 7]. Użyj MATLAB do obliczenia:
    sumy, różnicy, iloczynu macierzowego, iloczynu tablicowego tych macierzy. Oblicz także ich iloraz macierzowy i tablicowy. Sprawdź wyrywkowo wyniki wykonując rachunki bez MATLAB.
  3. Wygeneruj 50 liczb losowych z przedziału < 0, 1 > i przedstaw ten zbiór liczby na wykresie.
  4. Utwórz funkcję inline f(x) = x + cos(x2 ) i przedstaw ją na wykresie dla -5 < x < 5.
  5. Utwórz wektor o n elementach i oblicz średnią arytmetyczną elementów tego wektora.
  6. Napisz program, który pokazuje przebiegi następujących funkcji: cos(x), cos(2x), cos(4x), cos(6x), cos(8x) na jednym, wspólnym wykresie, dla 0 ≤ x ≤ 2 π. Skorzystaj z systemu pomocy MATLAB (komenda help plot), aby poznać więcej opcji funkcji plot( ).
  7. Narysuj wykres wielomianu x5 + 4x3 − 2x + 4, używając do tego celu 50 punktów z przedziału −4 < x < 4.
  8. Wypisz wszystkie kwadraty n2 tych liczb naturalnych n, które spełniają warunek n3 < 2200.