Matlab/Octave – wprowadzenie
Tomasz Sobiech,
Politechnika Warszawska, Wydział Fizyki2 marca 2015
S
krypt ten ma na celu zapoznanie państwa z działaniem i podstawową pracą z Matlab/Octave, czyli obliczeniach na macierzach i wizualizacji wyników. Działania tablicowe omówione zostały już w poprzedniej części, jednak ze względu na ich duże znaczenie i problemy w ich stosowaniu przez po- czątkujących, wracamy do nich raz jeszcze. W skrypcie wprowadzone zostało ważne pojęcie (bardzo ważne) wektoryzacji, pomimo początkowych trudności, musi zostać przyswojone przez każdego użytkow- nika wszystkich języków skryptowych Matlab/Octave, Python, Julia itd.1 Operatory (raz jeszcze)
1.1 Operator zakresu
Składnia operatora : wygląda następująco start : krok : nie większe niż. Przykład:
>> a = 3:0.7:7 a =
3.0000 3.7000 4.4000 5.1000 5.8000 6.5000
Domyślna wartość kroku to jeden, składnia upraszcza się wtedy do postaci:
>> a = 3:7 a =
3 4 5 6 7
W kontekście macierzy operator ten może przybrać jeszcze jedną formę:
>> A = zeros(3,4);
>> A(:,2) = 1 A =
0 1 0 0
0 1 0 0
0 1 0 0
Zapis A(:,2) weź wszystkie wartości z kolumny 2. Operatora zakresu możemy też używać z typem danych char:
>> s = ’a’:’f’
s = abcdef
1.2 Działania arytmetyczne
1.2.1 Operatory macierzowe
Możliwość wykonywania działań macierzowych jest absolutnie niezbędna dla programisty-fizyka, środowisko
Matlab/Octave pozwala na wykonywanie takich obliczeń w szybki sposób przez zdefiniowane operatory (tab. 1), bez pisania zbędnych pętli for.
Tablica 1: Operatory macierzowe
Operacja Operator
Dodawanie +
Odejmowanie -
Mnożenie *
Potęgowanie ^
Lewostronne dzielenie \ Prawostronne dzielenie /
Jeszcze słowo komentarza á propos dzielenia, oprócz dzielenia prawostronnego (/) mamy do dyspozycji dzielenie lewostronne (\). Dzielenie takie stosowane jest w równaniach macierzowych. W szczególności polecenie x = A\b rozwiązuje równanie Ax = b. Jest ono odpowiednikiem inv(A)*b, ale działa szybciej i numerycznie stabilniej.
1.2.2 Operatory tablicowe
Prawdopodobnie dużo częściej wykorzystywane są operatory tablicowe (tab. 2). Za ich pomocą możemy uzyskać takie struktury jak v.*w = [v1w1, v2w2, ..., vnwn], czy v.^w = [v1w1, vw22, ..., vwnn]. Te same zasady stosują się do macierzy. Dla dwóch macierzy A i B polecenie C = A.*B, w wyniku daje macierz z elementami Cij = AijBij.
Tablica 2: Działania tablicowe
Operacja Operator
Dodawanie +
Odejmowanie -
Mnożenie .*
Potęgowanie .^
Lewostronne dzielenie .\
Prawostronne dzielenie ./
1.3 Operatory relacji
Działania z użyciem operatorów relacji dają wynik w postaci macierzy o tej samej wielkości co argumenty, gdzie 1 oznacza prawdziwość relacji, a 0 oznacza fałsz.
Tablica 3: Operatory relacji
Operacja Operator
Mniejszy niż <
Mniejszy lub równy <=
Większy niż >
Większy lub równy >=
Równy ==
Różny od =
Przykłady:
>> x = [1, 5, 3, 7];
>> y = [0, 2, 8, 7];
>> k = x<y k =
0 0 1 0
>> k = x>=y k =
1 1 0 1
>> k = x~=y k =
1 1 1 0
1.4 Operatory logiczne
Operatory logicznie działają analogicznie do operatorów relacji, zwracają macierz o tej samej wielkości co argumenty:
Tablica 4: Operatory logiczne
Operacja Operator
Koniunkcja &
Alternatywa |
Negacja ~
Alternatywa wykluczająca xor
2 Praca z macierzami
2.1 Indeksowanie macierzy
W środowisku Matlab/Octave istnieją dwa sposoby indeksowania tablic (patrz rys. 1) tzw. index oraz subscripts, jak to działa sprawdzimy na prostym przykładzie:
>> A = rand(4, 5) A =
0.6557 0.6787 0.6555 0.2769 0.6948 0.0357 0.7577 0.1712 0.0462 0.3171 0.8491 0.7431 0.7060 0.0971 0.9502 0.9340 0.3922 0.0318 0.8235 0.0344
>> A(4) %w notacji index ans =
0.9340
>> A(4,1) %to samo w notacji subscripts ans =
0.9340
>> A(5) %w notacji index ans =
0.6787
>> A(1, 2) %to samo w notacji subscripts ans =
0.6787
Rysunek 1: Dwa sposoby indeksowania macierzy
Oczywiście tyczy się to też macierzy o większej ilości wymiarów. Do konwersji pomiędzy tymi dwoma stylami indeksowania macierzy można wykorzystać funkcję ind2sub oraz sub2ind. Przykład (Macierz A z poprzedniego przykładu):
>> idx = [sub2ind(size(A), 3, 2), sub2ind(size(A), 4, 4), sub2ind(size(A), 4, 5)] ...
%tworzymy wektor 3 indeksów które chcemy wybrać z macierzy idx =
7 16 20
>> A(idx) %wypisujemy podane elementy ans =
0.7431 0.8235 0.0344
>> A(idx) = 0; %nadpisujemy wartość podanych elementów
>> A A =
0.6557 0.6787 0.6555 0.2769 0.6948 0.0357 0.7577 0.1712 0.0462 0.3171
0.8491 0 0.7060 0.0971 0.9502
0.9340 0.3922 0.0318 0 0
Możemy też podawać w obu konwencjach wektory. W konwencji index przykład jest już wyżej, z macierzy wybierane są wszystkie elementy, których index jest podany w wektorze. Trochę więcej przykładów:
>> A([1, 3, 5, 6, 7]) ans =
0.6557 0.8491 0.6787 0.7577 0
>> A(1:2:20) %co drugi element ans =
Columns 1 through 9
0.6557 0.8491 0.6787 0 0.6555 0.7060 0.2769 0.0971 0.6948
Column 10 0.9502
>> i = 5:5:20 %co piąty indeks od piątego do dwudziestego i =
5 10 15 20
>> A(i) ans =
0.6787 0.1712 0.0971 0
W konwencji subscripts zwracany jest produkt kartezjański podanych indeksów:
>> A([1, 2], [1, 2]) %produkt kartezjański to elementy (1, 1), (1, 2), (2, 1), (2, 2) ans =
0.6557 0.6787 0.0357 0.7577
>> A(1:3,1:3) %wykorzystując operator zakresu ans =
0.6557 0.6787 0.6555 0.0357 0.7577 0.1712
0.8491 0 0.7060
>> A(1:2:3,1:3) %co drugi wiersz od pierwszego do trzeciego i od pierwszej kolumny do trzeciej ans =
0.6557 0.6787 0.6555
0.8491 0 0.7060
Należy tu jeszcze wspomnieć o kolejnym ważnym elemencie dotyczącym indeksowania, mianowicie słowo kluczowe end. Standardowo słowo end kończy blok kodu po for, if itd., w tym kontekście indeksowania ma ono inne znaczenie, tj.
oznacza ostatni element macierzy.
>> A(end) ans =
0
>> A(end-4:end) %ostatnie pięć elementów ans =
0 0.6948 0.3171 0.9502 0
>> A(end-1:end, 1) %korzystając z subscripts: przedostatni i ostatni element pierwszej kolumny ans =
0.8491 0.9340
>> A(2:end, 1:end-3) ans =
0.0357 0.7577
0.8491 0
0.9340 0.3922
Typowe zastosowanie wybierania indeksów z którym na pewno się spotkacie:
>> x = 0:0.01:10*pi;
>> y = sin(x);
>> idx = find(abs(y) < 5e-3); ...
%funkcja find zwraca indeksy elementów wektora równych 1 (działanie operatora relacji już znacie)
>> plot(x, y, ’b-’, x(idx), y(idx), ’ro’)
Na rysunku 2 przedstawiony jest wynik powyższych operacji.
2.2 Wektoryzacja
Przez wektoryzację rozumiemy nadanie obliczeniom takiej struktury, w której można zastosować zmienne wektorowe lub tablicowe z operatorami tablicowymi zamiast seryjnych obliczeń skalarnych. Jako przykład rozważ aproksymację funkcji wykładniczej za pomocą pierwszych dziesięciu wyrazów rozwinięcia w szereg
ex=X
k
xk−1 (k − 1)!
Można to liczyć w pętli, co jest czasochłonne, gdyż interpreter musi dziesięć razy wykonać tą czynność (interpreter jest wolny). W formie zwektoryzowanej interpreter parsuje jedną linijkę, reszta liczy się w wewnętrznych funkcjach
Rysunek 2: Miejsca zerowe funkcji f (x) = sin(x)
Matlab/Octave napisanych w szybkich językach typu C/C++
>> x = 1;
>> k=1:10;
>> e = sum(x.^(k-1)./factorial(k-1)) e =
2.7183
>>%lub w inny sposób po wykonaniu drobnych przekształceń w głowie
>> format long
>> sum(1./factorial(0:10))-exp(1) ans =
-2.731266102173890e-08
>> sum(1./factorial(0:10))-exp(1) %policzmy błąd takiego przybliżenia ans =
-2.731266102173890e-08
>> sum(1./factorial(0:20))-exp(1) %zwiększmy nieco ilość wyrazów szeregu ans =
0
W ten oto prosty sposób policzyliśmy liczbę Nepera. Stosować wektoryzację należy zawsze (o ile to jest możliwe)! Inny przykład, odtworzymy sobie działanie funkcji diff, i policzmy pochodną funkcji:
>> plot(x, y, x, Dy)
>> dx = 0.001;
>> x = 0:dx:2*pi;
>> y = sin(x);
>> Dy = (y(2:end)-y(1:end-1))./dx; ...
%diff => y(2:end)-y(1:end-1), dzieląc przez dx otrzymujemy pochodną
>> Dy = [Dy, Dy(end)]; %powielenie ostatniego elementu, żeby zgadzały się długości wektorów
>> plot(x, y, x, Dy)
Rysunek 3: Funkcja f (x) = sin(x) oraz jej pochodna f0(x) = cos(x)
2.3 Trudniejsze przykłady
Rozwiążmy teraz typowy fizyczny problem tj. rozwiążemy dwuwymiarowe równanie Poissona (∇u(x, y) = f (x, y) ∀(x,y)∈ h0, 100i × h0, 100i = Ω
u(x, y)
∂Ω= 0
gdzie
f (x, y) =
(100 ∀(x,y)∈ h49, 51i × h10, 90i = ω
0 ∀(x,y)∈ Ω\ω
Przedstawione tu rozumowanie jest nieco uproszczone, więcej o rozwiązywaniu równań różniczkowych można dowiedzieć się na wykładnie z Metod Numerycznych lub w internecie. Dyskretyzując równanie dla prostokątnej siatki otrzymujemy:
ui−1,j+ ui,j−1+ ui+1,j+ ui,j+1
4 − ui,j= fi,j (1)
u0,j = u100,j = ui,0= ui,100= 0 (2)
Obliczenie tego wyrażenia w C wymagałoby iteracji po całej macierzy w dwóch pętlach for, możemy jednak to zwektoryzować. Graficznie pierwszy składnik tego równania można przedstawić jak poniżej:
1 6 11 16 21 2 7 12 17 22 3 8 13 18 23 4 9 14 19 24 5 10 15 20 25
1 6 11 16 21 2 7 12 17 22 3 8 13 18 23 4 9 14 19 24 5 10 15 20 25
+
1 6 11 16 21 2 7 12 17 22 3 8 13 18 23 4 9 14 19 24 5 10 15 20 25
1 6 11 16 21 2 7 12 17 22 3 8 13 18 23 4 9 14 19 24 5 10 15 20 25
+
1 6 11 16 21 2 7 12 17 22 3 8 13 18 23 4 9 14 19 24 5 10 15 20 25
1 6 11 16 21 2 7 12 17 22 3 8 13 18 23 4 9 14 19 24 5 10 15 20 25
+
1 6 11 16 21 2 7 12 17 22 3 8 13 18 23 4 9 14 19 24 5 10 15 20 25
1 6 11 16 21 2 7 12 17 22 3 8 13 18 23 4 9 14 19 24 5 10 15 20 25
+
Tak powstałą nową macierz należy jeszcze podzielić przez cztery i odjąć od niej macierz początkową.
1 close all, clear all, clc %zamknij wszystkie wykresy, wyczysc zmienne, wyczysc linie komend
2
3 %definicja geometrii
4 N = 100;
5 u = zeros(N, N);
6 %rysuje prostokat
7 %wybieramy podmacierz i przypisujemy ladunek
8 %floor (podloga) zaokraglenie, aby miec pewnosc ze bedzie liczba calkowita
9 w1 = floor((N/2))-1;
10 w2 = floor((N/2))+1;
11 u(10:N-10, w1:w2) = 100;
12
13 for i = 0:1e4
14 %liczymy laplasjan (to tu cala magia)
15 nabla = zeros(N, N);
16 nabla(1:end-1,:) = nabla(1:end-1,:) + u(2:end,:);
17 nabla(2:end,:) = nabla(2:end,:) + u(1:end-1,:);
18 nabla(:,1:end-1) = nabla(:,1:end-1) + u(:,2:end);
19 nabla(:,2:end) = nabla(:,2:end) + u(:,1:end-1);
20 nabla = nabla./4.0;
21 nabla = nabla - u;
22
23 u = u + nabla;
24
25 %warunki brzegowe
26 u(1,:) = zeros(1, N);
27 u(end,:) = zeros(1, N);
28 u(:,1) = zeros(1, N);
u(:,end) = zeros(1, N);
30
31 u(10:N-10, w1:w2) = 100;
32 end
33
34 contourf(u); hold on;
35 rectangle('Position',[w1, 10, w2-w1, N-20],'FaceColor','k')
Rysunek 4: Rozwiązanie równania
Na koniec jeszcze przykład rozwiązania jednego z zadań z KADD.
Zadanie 1. Wygenerować macierz 3D zgodnie ze wzorem v = x ∗ sin x2− y2− z2 w przedziale [−2, 2] z krokiem 0.05.
Wyświetlić ją na ekranie.
1 close all, clear all, clc %zamknij wszystkie wykresy, wyczysc zmienne, wyczysc linie komend
2
3 [X,Y,Z] = meshgrid(-2:0.05:2,-2:0.05:2,-2:2); %tworzymy 3 macierze 3D
4 V = X.*sin(X.^2-Y.^2-Z.^2); %tworzymy nowa macierz 3D z wartosciami ze wzoru
5
6 for i = 1:5
7 % conturf rysuje wykres w formie mapy
8 [¬,h] = contourf(X(:,:,3), Y(:,:,3), V(:,:,3)); hold on;
9
10 %tu troche magicznych funkcji w sumie nieistotne
11 hh = get(h,'Children');
12 set(hh, {'ZData'}, cellfun(@(x) (i-3)*ones(size(x)), get(hh,{'XData'}), 'UniformOutput',false))
13 end
Uwaga 1. Rysowanie kilku map na jednym wykresie nie zadziała w Octave, przynajmniej nie w ten sposób. Jednak to nieduża strata bo ten sposób jest bardzo nie czytelny.
3 Wizualizacja
W tym rozdziale przedstawione zostały różne przykłady wizualizacji danych w Matlab/Octave. Najbardziej użyteczną funkcją generującą wykresy jest plot, jej składnia wygląda następująco:
plot(vec_x1, vec_y1, ’opcje stylu’, vec_x1, vec_y1, ’opcje stylu’, ...)
Rysunek 5: Mapy kostki
Przykładowo:
plot(y) %tworzy domyślny wykres na osi x odłożone są ideksy wektora y plot(x,y,’--’) %zamiast linii ciągłej jest przerywana
plot(x,y,’ro’) %rysuje czerwone kółka
Tablica 5: Opcje stylów
Opcje koloru Opcje stylu linii Opcje stylu znacznika
yżółty -linia ciągła +symbol plusa
mpurpurowy –linia przerywana okółko cgranatowy : linia kropkowana *gwiazdka rczerwony -. linia kreskowo-kropkowa x znak x
gzielony . kropka
bniebieski ^daszek
wbiały skwadrat
kczarny drąb
3.1 Etykiety, tytuły, legendy i inne
Wykresy można opisywać za pomocą poleceń:
xlabel(’napis x’) % tytuł osi x ylabel(’napis y’) % tytuł osi y title(’tytul’) % tytuł wykresu
text(x, y, ’napis’) % umieszcza na wykresie napis w pozycji (x,,y) Legendę można utworzyć, jak nie trudno się zgadnąć, za pomocą funkcji legend
legend(’linia 1’, ’linia 2’, ...) % tworzy legendę zawierającą etykiety ’linia 1’, ...
legend(’StylLinii1’, ’linia 1’, ...) % przypisuje każdej etykiecie styl linii
legend(..., pos) % pos = 1 (lub inna wartość) ustawia pozycję legendy
legend off % wyłącza legendę
Ustawiać zakresy osi możemy za pomocą polecenia axis
axis([x_min, x_max]) % ustawia zakres osi x axis([x_min, x_max, y_min, y_max]) % ustawia zakres osi x i y
axis(’equal’) % ustawia jednakową skalę na obu osiach itd.
3.2 Wykresy nakładane
Istnieje kilka sposobów narysowania wielu linii na jednym wykresie, oto przykład:
1 t = linspace(0, 2*pi, 100);
2 y1 = sin(t);
3 y2 = t;
4 y3 = t - (t.^3)/6 + (t.^5)/120;
5
6 %rysuje y1 jako ciagla linie (domyslnie)
7 plot(t,y1)
8 % %dodaje y2 jako linie przerywana
9 line(t,y2, 'linestyle', '--')
10 % %dodaje y3 jako serie kolek
11 line(t,y3, 'linestyle', 'o')
12
13 %inaczej ale robi to samo
14 plot(t,y1)
15 hold on
16 plot(t,y2, 'linestyle', '--')
17 plot(t,y3, 'linestyle', 'o')
18 hold off
19
20 %w ten sposob podobnie ale dodatkowo automatycznie zmienia kolory linii
21 plot(t,y1,t,y2, '--',t,y3, 'o')
22
23 axis([0 5 -1 5]) %nowe zakresy osi
24 xlabel('t')
25 ylabel('aproksymacja sin(t)')
26 title('Wykresy nakladane')
27 legend('sin(t)', 'aproksymacja liniowa', 'aproksymacja piatego rzedu')
3.3 Tworzenie wykres ´ ow r ´ ownoległych przy u˙zyciu subplot
Do rysowania kilku wykresów w jednym oknie służy polecenie subplot. Funkcja ta dzieli okno na siatkę (n, m) podwykresów, na przykład:
subplot(3,2,1) % trzeci argument oznacza na którym wykresie aktualnie rysujemy plot(x)
subplot(3,2,2) plot(y)
...
1 t = linspace(0, 8*pi, 200);
2 y = t.*sin(t);
3
4 figure(1) %tworzy okno
5 subplot(1, 2, 1)% dzieli okno na macierz (1,2) wykresow rysuje na 1
6 area(t, y);
7 subplot(1, 2, 2)% teraz rysuje na 2
8 t = linspace(0, 2*pi, 200);
9 y = sqrt(abs(2.*sin(5.*t)));
10 polar(t, y)
11
12 figure(2) %tworzy drugie okno
13 subplot(1, 2, 1)
14 hist(randn(1, 1000));
15 subplot(1, 2, 2)
16 t = linspace(0, 2*pi, 500);
17 r = sqrt(abs(2.*sin(5.*t)));
18 x = r.*cos(t);
19 y = r.*sin(t);
20 fill(x, y, 'r')
