pobieranie: Matlab 1
Pobierz dokument pdfPrzeglądaj wersję html pliku:
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Tryb interpretera
2 * 3
ans = 6 2 * 3; pi
ans = 3.1416 sin(pi / 2) 1 [1 4; 6 8] 1 6 6: 2: 20 6 8 10 12 14 16 18 20 4 8
ans =
ans =
ans =
Wprowadzenie do Matlaba, folia 1
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
b
? Undefined function or variable ’b’ A = [1 4; 6 8]
A = 1 6 A - 1 4 8
ans = 0 5 A * A 3 7
ans = 25 54 sin(A)
36 88
ans = 0.8415 -0.2794
-0.7568 0.9894
Wprowadzenie do Matlaba, folia 2
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Programowanie
Zadanie 1. Napisa´ funkcj˛ znajdujaca c e ˛ ˛ pierwiastek równania liniowego. function x = rown1(a, b) if a == 0 if b ~= 0 x = []; else x = NaN; end; else x = -b / a; end; Uwaga: T˛ funkcj˛ zapisuje si˛ w pliku rown1.m ! e e e Efekt działania: rown1(1, 2) ans = -2
Wprowadzenie do Matlaba, folia 3
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Zadanie 2. Napisa´ funkcj˛ rozwiazujaca c e ˛ ˛ ˛ równania liniowe i kwadratowe. function x = rown21(a, b, c) if nargin == 2 x = rown1(a, b); elseif a == 0 x = rown1(b, c); else delta = b * b - 4 * a * c; if delta >= 0 if b > 0 x = (-b-sqrt(delta))/(2*a); else x = (-b+sqrt(delta))/(2*a); end; if x == 0 x = [x 0]; else x = [x (c/a)/x]; end; end; end;
Wprowadzenie do Matlaba, folia 4
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Wskazana jest równie˙ diagnostyka błedów: z if nargin < 2 error(’Za malo parametrow’); end;
P˛ tla for i wektoryzacja e
p˛ tla for e for i = 1: 100 y(i)=sin(i); end; Co wybiera´ ? Zawsze wektoryzacj˛ ! c e Zadanie. Napisa´ funkcj˛ obliczajaca sum˛ c e ˛ ˛ e p-tych pot˛ g n pierwszych liczb naturalnych. e Wersja „klasyczna”: function s = sump(n, p) s = 0; for i = 1: n s = s + i^p; end;
¡¡¡¡¡
wektoryzacja y = sin(1: 100);
Wprowadzenie do Matlaba, folia 5
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Wersja zwektoryzowana: function s = sumpv(n, p) s = sum((1: n).^p); Wektoryzacja jest mo˙ liwa, je˙ eli operacje w p˛ tli z z e nie zale˙ a od porzadku, w ktorym sa wykonywane. z˛ ˛ ˛
Specyfika p˛ tli for e
Zadanie. Co b˛ dzie rezultatem poni˙ szych e z instrukcji: for i=12:-2:1; disp(i), end; for i=12:-2:1; disp(i), ... i = 0; end; for i=12:-2:1; disp(i); ... if i == 8, break; end; end;
Wprowadzenie do Matlaba, folia 6
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Zadanie. Napisa´ funkcj˛ wyznaczajaca n-ty c e ˛ ˛ element ciagu Fibonacciego. ˛ Wersja prostsza: function f = fibs(n) if n == 1 | n == 2 f = 1; else a = 1; b = 1; for i = 3: n c = a + b; a = b; b = c; end; f = c; end; Wersja bardziej efektywna: function f = fib(n) persistent rez if length(rez) < 2 rez = [1 1]; end;
Wprowadzenie do Matlaba, folia 7
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
if length(rez) < max(n) for i = length(rez)+1: max(n) rez(i)=rez(i-1)+rez(i-2); end; end; f = rez(n); Efekt działania: fib([3, 9, 4]) ans = 2 34 3
Zadanie. Napisa´ funkcj˛ obliczajaca n c e ˛ ˛ pierwszych linii trójkata Pascala. ˛ function c = trojkat(n) c = zeros(n, n); c(:,1) = 1; for i = 2:n for j = 2:i c(i,j) = c(i-1,j-1)+c(i-1,j); end; end;
Wprowadzenie do Matlaba, folia 8
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Efekt działania: trojkat(6) ans = 1 1 1 1 1 1 0 1 2 3 4 5 0 0 1 3 6 10 0 0 0 1 4 10 0 0 0 0 1 5 0 0 0 0 0 1
Sposoby wektoryzacji
Prze´led´ my poni˙ sza sesj˛ : s z z ˛ e a = 1: 16; a = reshape(a, 4, 4) a = 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
Wprowadzenie do Matlaba, folia 9
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
a1 = a(:, 4: -1: 1)
a1 = 13 14 15 16 9 10 11 12 5 6 7 8 1 2 3 4
a2 = a(4: -1: 1, 4: -1: 1)
a2 = 16 15 14 13 a3 = a’ 12 11 10 9 8 7 6 5 4 3 2 1
a3 = 1 5 9 13 2 6 10 14 3 7 11 15 4 8 12 16
Wprowadzenie do Matlaba, folia 10
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
To samo wykonuje zreszta predefiniowana funkcja ˛ nchoosek. Zadanie. Utworzy´ wektor, w którym ka˙ demu c z indeksowi parzystemu p b˛ dzie odpowiada´ e c kwadrat, a indeksowi nieparzystemu – sze´cian s p-tej liczby naturalnej. v(1: 2: n) = (1: 2: n).^3; v(2: 2: n) = (2: 2: n).^2; Zadanie. Majac dany wektor x, utworzy´ tablic˛ ˛ c e a o n wierszach b˛ dacych kopiami x. e ˛ x = 1: 6; n = 5; a = x(ones(1, n), :) a = 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4 5 5 5 5 5 6 6 6 6 6
Wprowadzenie do Matlaba, folia 11
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Zadanie. Dla tablicy A = 3 -1 -1 1 1 2 0 3 0 -3 2 5 3 2 0 -4 2 -1 -4 1 3 2 5 5 -4 -1 0 -2 1 -5 1 2 5 3 2 -1
wygenerowa´ tablic˛ zło˙ ona z elementów le˙ acych c e z ˛ z˛ w wierszach od trzeciego do piatego i kolumnach 6, ˛ 1 i 3. B = A(3: 5, [6 1 3]) B = 5 3 2 -1 1 1 0 -4 2
Jak w tej tablicy wyzerowa´ wiersze nieparzyste, a c elementy wi˛ ksze od 3 zamieni´ na 1. e c A(1: 2: 6, :) = 0; A(A > 3) = 1;
Wprowadzenie do Matlaba, folia 12
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Zadanie. Narysowa´ wykres funkcji c
w przedziale nieciagło´ci w punkcie ˛ s
x = linspace(pi-1, pi+1, 50); y = 3*x.^2+log((x-pi).^2)/pi^4+1; plot(x, y)
55
50
45
40
35
30
25
20
15
10
2
¡ ¤ ' ¨ # (! &%$! " ! ¨ ¦ ¤ ©§¤¥£ ¢ ¨ ¦ ¡
2.5 3 3.5
. Dlaczego nie wida´ c ?
4
4.5
Wprowadzenie do Matlaba, folia 13
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Zadanie. Narysowa´ powierzchni˛ c e
dla
x = linspace(-1, 1, 5)
x = -1.00
-0.50
y = x; [X, Y] = meshgrid(x, y)
X = -1.00 -1.00 -1.00 -1.00 -1.00 Y = -1.00 -1.00 -1.00 -1.00 -1.00 -0.50 -0.50 -0.50 -0.50 -0.50 0 0 0 0 0 0.50 0.50 0.50 0.50 0.50 1.00 1.00 1.00 1.00 1.00 -0.50 -0.50 -0.50 -0.50 -0.50 0 0 0 0 0 0.50 0.50 0.50 0.50 0.50 1.00 1.00 1.00 1.00 1.00
¦ §¤ ¡ ¦
. 0 0.50 1.00
Wprowadzenie do Matlaba, folia 14
! £ ¤ ¡¢! ¡ #
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Z = -X.^2 - Y.^2
Z = -2.00 -1.25 -1.00 -1.25 -2.00 -1.25 -0.50 -0.25 -0.50 -1.25 -1.00 -0.25 0 -0.25 -1.00 -1.25 -0.50 -0.25 -0.50 -1.25 -2.00 -1.25 -1.00 -1.25 -2.00
surf(x, y, Z) colormap(gray)
0
−0.5
−1
−1.5
−2 1 0.5 0 0 −0.5 −1 −1 −0.5 0.5 1
Wprowadzenie do Matlaba, folia 15
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Zadanie. Napisa´ funkcj˛ v = vnd(x), która c e na bazie wektora x o składowych utworzy nast˛ pujaca macierz (tzw. macierz e ˛ ˛ Vandermonde’a):
function v = vnd(x) x = x(:); % x powinno byc kolumna v = ones(size(x)); for i = 1: length(x) - 1 v = [v(:, 1).*x v]; end; vnd([2 3 4 5])
ans = 8 27 64 125 4 9 16 25 2 3 4 5 1 1 1 1
Wprowadzenie do Matlaba, folia 16
!
£ ¤
¤ ¦£ ¤
¨
¨
¨
¨
¨
. . .
¨ £ £
. . .
¤¤%# ¡¢¢¡¢# ¤ £ ¡
. . . . . .
! !
©
¤ ©¦ ¤ ¤ ¦ ¤
¨ ¨
¨ ¨ ¨
¨
©£
£
¤ ¤ ¤
¦ ¦ § ¦ ¦ ¦
¥¦
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
P˛ tla while i rekurencja e
P˛ tla niesko´ czona z instrukcjami break: e n while 1 ... if condition_1 break; end; ... if condition_2 break; end end; Zadanie. Napisa´ funkcj˛ szukajaca w wektorze c e ˛ ˛ t ostatniego wystapienia zadanego elementu. ˛ function i = szuk1(t, x) i = 0; for j = length(t): -1: 1 if t(j) == x i = j;
Wprowadzenie do Matlaba, folia 17
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
break; end end; function i = szuk2(t, x) i = find(x == t); if isempty(i) i = 0; else i = i(end); end; function i = szuk3(t, x) i = length(t); while 1 if i == 0 break; end; if t(i) == x break; end; i = i - 1; end;
Wprowadzenie do Matlaba, folia 18
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
t = randperm(6000); tic; szuk1(t, 3000); toc
elapsed_time = 0.0400 tic; szuk2(t, 3000); toc
elapsed_time = 0.0110 tic; szuk3(t, 3000); toc
elapsed_time = 0.1000 Zadanie. Znale´ c pierwiastki równania z´
¡
Wprowadzenie do Matlaba, folia 19
Zastosowa´ algorytm: c 1. Znale´ c pierwiastek rzeczywisty z´ Newtona.
© ©
. (Łatwo 2. Podzieli´ trójmian przez c pokaza´ , ze odpowiedni iloraz ma posta´ c ˙ c , gdzie: , .) 3. Rozwiaza´ otrzymane równanie kwadratowe ˛ c otrzymujac pozostałe dwa pierwiastki. ˛
£
¤
¦
¨©¤ ¡
¡
¥
¤ ¤
¨ ¤ ¤ ¨ §¤ £ ¨ ¢¤ ¦ ¡
§ ¨ ¤
metoda ˛
¨ § ¡ ¤ ¨ § ¨ ¦ ¤ ¤©
¨
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
function x1 = newton1(a, b, c, d) epsi = 1.0e-5; x1 = 0; while 1 fprim = (3*a*x1+2*b)*x1+c; if fprim == 0 x1 = -b/a; else dx = -(((a*x1+b)*x1+c)*x1+d)... /fprim; x1 = x1 + dx; if abs(dx) < epsi break; end end end; function x = rown3(a, b, c, d) if a == 0 x = rown21(b, c, d); else x = newton1(a, b, c, d); beta = b + a * x;
Wprowadzenie do Matlaba, folia 20
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
gamma = c + beta * x; x = [x rown21(a, beta, gamma)]; end; rown3(2, -12, 22, -12)
ans = 1.0000
3.0000
2.0000
Uwaga! W praktyce stosuje si˛ funkcj˛ roots: e e roots([2, -12, 22, -12]) ans = 3.0000 2.0000 1.0000 Ogólny schemat rekurencji: function y = p(n) if n == 1 y = wyrazenie_explicite; else y = zwiazek_z( p(n - 1) ); end;
Wprowadzenie do Matlaba, folia 21
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Zadanie. Napisa´ rekurencyjna wersj˛ funkcji c ˛ e wyznaczajacej warto´c , ˛ s´ . function y = potega(n, p) if p == 0 y = 1; elseif rem(p, 2) == 1 y = n * potega(n*n, floor(p/2)); else y = potega(n*n, floor(p/2)); end; Zadanie. Napisa´ rekurencyjna wersj˛ funkcji c ˛ e szukajacej w wektorze t ostatniego wystapienia ˛ ˛ zadanego elementu. function i = szuk4(t, x) if isempty(t) i = 0; elseif t(end) == x i = length(t); else i = szuk4(t(1: end-1), x); end;
¦
§ ¦§¤# ¢ ¥ £ ¡
Wprowadzenie do Matlaba, folia 22
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Zadanie. Napisa´ funkcj˛ obliczajaca c e ˛ ˛ przybli˙ ona warto´c całek adaptacyjna metoda z ˛ s´ ˛ ˛ Simpsona. Wykorzysta´ ja do obliczenia całki c ˛
Rozwiazanie: Wzór Simpsona jest postaci ˛
mo˙ na obliczy´ na dwa sposoby: z c
gdzie:
´ – srodek
.
Wprowadzenie do Matlaba, folia 23
' © %$#¤%%!"¤ " ! #
Aby obliczy´ całk˛ z dokładno´cia c e s ˛ , dzielimy przedział na podprzedziały w taki sposób, aby na ka˙ dym z nich osiagna´ z ˛ ˛c dokładno´c s´ .
' © %$#¤%#%#¤ " ! ! ¥ © %$#%# ¦ 67© $%! ¤%# ©¨ ¨ ¦ 67© ! ¤
§
§
©£
¤ ¢ ! ¨ ¦ ¤ ¡
¥
¨
§ ¨ ¥ § # #
¨ ¤ ¤ ¦
#¤ ! ' £ # "
©
$ %!
$ %!
¥
©¨ §¤ ¢ ¦ ¡
¦ %$! ¤ 6 7© ¥ ¤ ¤ %# !#%# ©¨ ¡ ¥ ¤%# ! %# © ¨ ¡ ¤ ¥ ¤ ¢ ¤ ¡
©
¨
$! %#¤
¥
¤
¥
¡
! 4 ) ' 231(& 0) ' ! 54
£ ¤
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Mo˙ na udowodni´ , ze z c ˙
gdzie:
– dokładna warto´c odpowiedniej całki. s´
Wynika stad algorytm wyznaczania ˛ : Oblicza si˛ warto´ci e s
oraz „bład” ˛ . Je˙ eli „bład” jest wi˛ kszy z ˛ e ni˙ z , dokonuje si˛ podziału e na podprzedziały i ,a nast˛ pnie powtarza powy˙ sze obliczenia na ka˙ dym e z z z nich z osobna. Implementacja: Funkcj˛ podcałkowa mo˙ na e ˛ z zdefiniowa´ dwoma sposobami: c w pliku fun.m function y = fun(x) y = x.^4.*log(x+sqrt(x.*x+1));
' %$#¤%# ¦ %$! ¤ " ' ¦ 67© ! ' © %$#¤%%! © ¤ " 67© ! #
¥ %$#%# ¦ $%! ¤%# ! ¤ 67© © ¥ ¦ 7© $%! %# #¤%# ¥ ¤ ! 6 © $%#%# #¤%# ! ¤ !
Wprowadzenie do Matlaba, folia 24
! 54
! ! 4 !
£
%! $
¦ ! 54 © ¨ © ¨ © ¨
¤# %%! ¤ " #¤ %#¤ ! ! $! ¦ £ ¨ © £© ¡ ¦£
¡ ¡
© %$#%%#%# ! ¤#! ¤
©
! ¨
£ £
! ¨
¥
¡ ¢©¨
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
lub w linii polece´ : n fun = inline( ’x.^4.*log(x+sqrt(x.*x+1))’); Moduł inicjujacy rekurencj˛ , zapisany w pliku ˛ e simps.m , ma posta´ c function int = simps(f, przedz, epsi) int = []; a = przedz(1); b = przedz(2); if nargin < 3 epsi = 0.5*sqrt(eps)/(b-a); else epsi = 0.5*epsi/(b-a); end; przyb_pocz = (feval(f,a) ... +4*feval(f,0.5*(a+b)) ... + feval(f,b))*(b-a)/6.0; int = simps_rek(f, przedz, ... przyb_pocz, epsi); W tym samym pliku zapisuje si˛ funkcj˛ e e wewn˛ trzna simps_rek realizujaca rekurencj˛ . e ˛ ˛ ˛ e
¡¡¡¡¡ ¡¡¡
Wprowadzenie do Matlaba, folia 25
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
% funkcja wewnetrzna function intr = simps_rek(f, ... przedz, stare_przyb, epsi) a = przedz(1); b = przedz(2); srod = (b+a)*0.5; h = (srod-a)/6.0; fsrod = feval(f,srod); s1 = (feval(f,a) ... + 4.0*feval(f,(a+srod)*0.5) ... + fsrod)*h; s2 = (fsrod ... + 4.0*feval(f,(srod+b)*0.5) ... + feval(f,b))*h; err = abs(s1+s2-stare_przyb); if err>=15.0*epsi*(b-a) s1 = simps_rek(f, [a, srod], ... s1, epsi); s2 = simps_rek(f, [srod, b], ... s2, epsi); end; intr = s1 + s2;
Wprowadzenie do Matlaba, folia 26
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Zauwa˙ my, ze wywołanie simps zale˙ y od z ˙ z sposobu zdefiniowania funkcji podcałkowej: ¥ z zastosowaniem pliku: simps(’fun’, [0 2], 1.0e-4) ans = 8.1534 ¥ za zastosowaniem inline: simps(fun, [0 2], 1.0e-4) ans = 8.1534
Ciagi znaków i pliki tekstowe ˛
Zadanie. Wy´wietli´ wszystkie znaki ASCII o s c kodach od 32 do 255, oprócz znaku o kodzie 127. Dlaczego nie wymaga si˛ wy´wietlenia pozostałych e s znaków? i = char(32:255); i(127) = []; i
Wprowadzenie do Matlaba, folia 27
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Zadanie. Zinterpretowa´ poni˙ sze rezultaty: c z abs(’zorro’) ans = 122
111
114
114
111
sin(’zorro’)
ans = 0.4987 -0.8646 0.7850 0.7850 -0.8646 sin zorro
ans = 0.4987 -0.8646 0.7850 0.7850 -0.8646 Zadanie. Napisa´ funkcj˛ okre´lajaca liczb˛ c e s ˛ ˛ e wystapie´ danego znaku. ˛ n function lwyst = ile(s, x) lwyst = sum(s == x); Zadanie. Napisa´ funkcj˛ okre´lajaca czy dane c e s ˛ ˛ słowo jest palindromem. function r = palindrom(s) r = strcmp(s, s(end: -1: 1));
Wprowadzenie do Matlaba, folia 28
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Zadanie. Napisa´ funkcj˛ zamieniajaca pierwsze c e ˛ ˛ litery słów danego ciagu znaków na du˙ e litery. ˛ z function s1 = duze_lit(s) s = [’ ’ s]; ind = findstr(s, ’ ’); s = [s ’ ’]; s(ind+1) = upper(s(ind+1)); s1 = s(2: end-1); duze_lit(’to jest
ciag’)
ans = To Jest
Ciag
Zadanie. Napisa´ funkcj˛ odwracajaca kolejno´c c e ˛ ˛ s´ słów w zdaniu. function s1 = odwroc(s) [pslowo reszta] = strtok(s); if isempty(pslowo) s1 = s; else s1 = [odwroc(reszta), ’ ’* ... ones(1,(length(reszta)>0)), pslowo]; end;
Wprowadzenie do Matlaba, folia 29
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
odwroc(’to jest jedno zdanie’)
ans = zdanie jedno jest to Zadanie. Napisa´ funkcj˛ zapisujaca w pliku ciag c e ˛ ˛ ˛ znaków oraz odpowiednia funkcj˛ czytajaca ten ˛ e ˛ ˛ ciag. ˛ function zapisz(s, nazwa) f = fopen(nazwa,’w’); fwrite(f, s); fclose(f); function s = czytaj(nazwa) f = fopen(nazwa,’r’); s = fread(f); fclose(f); s = char(s’);
Wprowadzenie do Matlaba, folia 30
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Formatowane wyj´cie: s x = cumprod(1.0e-4*ones(5, 1)); a = [x, exp(x), exp(-x), ... (exp(x)- exp(-x))./(2*x)]; fprintf(’%15s %15s %15s %15s\n’, ... ’x’, ’exp(x)’, ’exp(-x)’, ’poch.’); fprintf( ... ’%15.4e %15.4e %15.4e %15.4e\n’, a’); Rezultat: x 1.0000e-004 1.0000e-008 1.0000e-012 1.0000e-016 1.0000e-020 exp(-x) 9.9990e-001 1.0000e+000 1.0000e+000 1.0000e+000 1.0000e+000 exp(x) 1.0001e+000 1.0000e+000 1.0000e+000 1.0000e+000 1.0000e+000 poch. 1.0000e+000 1.0000e+000 1.0000e+000 5.5511e-001 0.0000e+000
Wprowadzenie do Matlaba, folia 31
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Bł˛ dy zaokraglen e ˛ ´
Zadanie. Wyja´ni´ rezultat wykonania s c poni˙ szego skryptu: z x y z z = = = = 1.0e+29; 1.0e-9; ((y + x) - x) / y % => 0 (y + (x - x)) / y % => 1
Co to jest eps? eps
ans = 2.2204e-016 2^(-52) ans = 2.2204e-016 1 + eps == 1 ans = 0
1 + eps/2 == 1 ans = 1
Wprowadzenie do Matlaba, folia 32
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Zadanie. Napisa´ program obliczajacy warto´ci c ˛ s dla , -20, -10, 0, 10, 20, 30. Wykorzysta´ w tym celu definicj˛ c e
i porówna´ rezultaty z generowanymi przez exp. c function s = expo( x ) t = 1.0; n = 1; s = 0.0; while 1 if ( abs(t) <= 1.0e-14 ) break; end; s = s + t; t = t * (x/n); n = n + 1; end Skrypt: for i = -3:3 fprintf(’%4d %10.3e %10.3e \n’, ... 10*i, exp(10*i), expo(10*i)); end
¢ ¤ ¥ ££
£
¤
¡
¡
'
¦
£ ¡ ¤
'
Wprowadzenie do Matlaba, folia 33
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
-30 -20 -10 0 10 20 30
9.358e-014 2.061e-009 4.540e-005 1.000e+000 2.203e+004 4.852e+008 1.069e+013
6.103e-006 6.148e-009 4.540e-005 1.000e+000 2.203e+004 4.852e+008 1.069e+013
Pytanie: Jak ograniczy´ bł˛ dy zaokragle´ ? c e ˛ n Odpowied´ : Ograniczy´ obliczenia do przedziału z c :
Wprowadzenie do Matlaba, folia 34
¤ © ¤ ¤ ! ¨ ¡ ¥¦ ¦ ¤ ! ¨ ¤ ! ¡ ¤ ¨ ¤
Zadanie. Porówna´ warto´ci funkcji c s w przedziale z warto´ciami tej samej funkcji, ale liczonymi wg s formuły
' (!
" #
¦ ¦ ¦¡
¦ ¦
!#
¨
!!!
¤
¤ #
¡
¥ ©
¦"
'
¦
' ' ' !
© ¥ § ¦ ¤ ¥
#
¢£¡
¤ ! ¡ ¤
¡
' ¥ ¥
' (! # ¦ "
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
function e = bin1(x) x = 1-x; e = x.*x.*x.*x.*x.*x; function e = bin2(x) x2 = x.*x; x4 = x2.*x2; e = 1 - 6.*x + 15.*x2 ... - 20.*x2.*x + 15.*x4 ... - 6.*x4.*x + x4.*x2; Skrypt: poc = 99990; kon = 100010; for i=(poc: kon)/100000 fprintf(’%15.5f %10.3e %10.3e\n’, ... i, bin1(i), bin2(i)); end; i=(poc: kon)/100000; plot(i, bin1(i)); pause; plot(i, bin2(i));
Wprowadzenie do Matlaba, folia 35
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
bin1.m:
1 x 10
−24
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0 0.9999
1
1
1
1.0001
bin2.m:
4 x 10
−15
2
0
−2
−4 0.9999
1
1
1
1.0001
Wprowadzenie do Matlaba, folia 36
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Ró˙ niczkowanie numeryczne z
Dla pierwszych pochodnych stosuje si˛ e
lub
Dla drugich pochodnych mamy:
for x = cumprod(0.1 * ones(1, 10)) fprintf(’%2.0e %12.10f %14.10f \n’, ... x, (exp(x) - exp(-x))/(2 * x), ... (exp(x) - 2 + exp(-x))/x^2); end;
¨
Wprowadzenie do Matlaba, folia 37
¡
¦ ¥ ¨
¡
¥
¨
¡
¥
¡
¡
¡
¥ ! ¦ ! ¡ ¥ ¥ ¨ ¦
¡
¥
© § ¤¥ ¨
¦ ! ¥ ¡ ¡ ' ¨ # " ¦
¥
¡ ¡
! ¡ ¥ ¡ ' ¨ # " ¦
¥ £ ¦ ¤ ¡ ¡ ¤ ¥ £ ¡
¥
¡
¥
¡
¥
£
¢
¥
¥ ¥
¥ ¥
¢ ¢
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Efekty sa intrygujace: ˛ ˛ 1e-001 1e-002 1e-003 1e-004 1e-005 1e-006 1e-007 1e-008 1e-009 1e-010 1.0016675002 1.0008336112 1.0000166667 1.0000083334 1.0000001667 1.0000000834 1.0000000017 1.0000000050 1.0000000000 0.9999989725 1.0000000000 0.9999778783 0.9999999995 0.9992007222 0.9999999939 0.0000000000 1.0000000272 111.0223024625 1.0000000827 0.0000000000
Wyja´nienie: s
¦ ¦
Wprowadzenie do Matlaba, folia 38
¢
¢ ¡
$
¡ ¢ ¡ ¡
¦ ! ¥ ¡ ¨ ¥ ¥ ¨ ¦ ! ¥ £¤¡¢ ¥ ¥ ¨ ¨ ¡
¢ ¡ ¤¥ ¦£ ¤¥ ¦£ ¡ ¡
$
¥
¡
¡
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Całkowanie numeryczne
Zadanie. Z zastosowaniem funkcji quad lub quad8 obliczy´ całk˛ c e
f1 = inline(... ’x.^4.*log(x+sqrt(x.*x+1))’); quad8(f1, 0, 2) ans = 8.1534 Zadanie. Z zastosowaniem funkcji dblquad obliczy´ całk˛ c e
function v = f2(x, y) v = exp(-x.*x-y.*y);
¤ ¢ ! ¨ ¦ ¤ ¡ ¡ ¤
¦
¡ ¢ ' ¡
¨ ¤ ¤ ¦
¦ "
"
Wprowadzenie do Matlaba, folia 39
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
dblquad(’f2’, 0, 5, 0, 5, ... [], [], ’quad8’) ans = 0.7854 Zadanie. Narysowa´ wykres funkcji c
w przedziale
.
f = inline(’cos(x.*sin (u))./pi’, ... ’u’, ’x’); x = linspace(0, 10, 100); for i = 1: 100 y(i) = quad(f, 0, pi, ... [], [], x(i)); end; plot(x, y);
¡¢ ¤ §¥£ ¦ ¤
¨©¦
¢
' ¦ ! # ¦"
!
¡
¤
Wprowadzenie do Matlaba, folia 40
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
1.2
1
0.8
0.6
0.4
0.2
0
−0.2
−0.4
−0.6
0
1
2
3
4
5
6
7
8
9
10
Metoda Monte-Carlo
– ciag niezale˙ nych zmiennych losowych o ˛ z jednakowym rozkładzie, ,
§¨
jednostajny na hiperkostce
¡ ¡ ¡ ¡
Wprowadzenie do Matlaba, folia 41
!
¡ ¡¢¢¡ ¡
§¨
¡ !"(! # ¦ " '
©
¡
¡ ¥
¥ ¡ ¢ ! # ¢¢¡¢# ¡ ¡ ©
¦
¨
¡
¡
¡
¡
¥
¡
©
©
¢
¨ ¨ ¢¢¨
¡
!
©
©
¥
©
¢
¡
£ ¤¡
¡
!
©
¡ ¢
¡
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Bład mo˙ na oszacowa´ za pomoca odchylenia ˛ z c ˛ standardowego z próby (zob. std) pomno˙ onego przez z . Zadanie. Napisa´ funkcj˛ generujaca c e ˛ ˛ macierz, której ka˙ dy wiersz zawiera prób˛ z z e rozkładu jednostajnego na danym przedziale.
function x = random(gran, n) x = rand(size(gran,1), n); dlug = gran(:, 2) - gran(:, 1); x = gran(:, ones(1, n)) ... + x.*dlug(:, ones(1, n));
Wprowadzenie do Matlaba, folia 42
!
¡
¡¢¡¢¡ ©
§ ¨
¤ # ¢¡¢¢# © ¤ ¡ ¡
¡ ¢ !
¥
¡¢¡¢¡¢ ©
©
#
© ¤
#
¡
¥
¦
¡ ¡ ¡
¥
¤ ¥¢
¤ £
¡
o hiperobj˛ to´ci e s
¡
¢!
jednostajny na hiperkostce
§ ¨
' ! £ #%! " ©
!
¡
¨¢¢¨ ¨
2
2
¢
¥
£
©
¢
!
¡
¡
¡
!
©
¡
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Zadanie. Napisa´ funkcj˛ przybli˙ ajaca całk˛ c e z ˛ ˛ e wielokrotna po zadanej hiperkostce. ˛ function [int, st] = mtc( ... fun, gran, npkt); if nargin <= 2 npkt = 10000; end; obj = prod(gran(:,2)-gran(:,1)); z = feval(fun, random(gran, npkt)); if nargout == 2 st = obj * std(z) / sqrt(npkt); end; int = obj * mean(z); Zadanie. Obliczy´ obj˛ to´c przeci˛ cia kul o c e s´ e ´ srodkach w i oraz promieniach odpowiednio 3 i 2. function y = ff(x) y = ((x(1,:)).^2 + x(2,:).^2 ... + x(3,:).^2 < 9)... & ((x(1,:)-2).^2 + x(2,:).^2 ... + x(3,:).^2 < 4);
¦# ¦# ¦# ¦# ¦
Wprowadzenie do Matlaba, folia 43
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
[int, st] = mtc(’ff’, ... [0 3; -2 2; -2 2], 100000) int = 24.6931 st = 0.0759
Równania ró˙ niczkowe z
Metoda Eulera
¡
gdzie:
Wprowadzenie do Matlaba, folia 44
¦
¡
¢
¨ ! %# ! ¤
¡
gdzie:
,
,
£
¡
©
£ %# ¡¢¢¡¢# © ¡¢ ¡ $ £ ¥ ¤ ¡ # ¤ # %¤ ¡
¥
¥
¡
¨ !¤ ¡ ©$ !¤ ¨ ! ¡ © $ !
¡ ¡
¥
¤
"
¡
function yh = euler(f, x, y, h) yh = y + h * feval(f, x, y);
Zadanie. Zaimplementowa´ wymienione metody. c
¡
¢
¡
¡
¡ ¨ ¦ ¨ ! ¡ ¥ ¦© © ! ! ! ¨ %# ! ¨ ¤ ¥ ¡ ! %# ! ¤ ¡ ©
¡
Metoda RK4
Metoda RK2
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Wprowadzenie do Matlaba, folia 45
¢
¡
¡
¡
¡
¡
¨ ! ¥ ¤ ¥ ¤ ¥ ¤
¡
¥
¡ © $ ! ¡ ¡ ¡ ¡ ¦ ¡
©
¨ ¨ ¡ ¨ ¦ ¨ © ¡ ! ¡ ¨ ! %# ¨ ! ¦ ! ¨ ! %# ! ¨ ! ¨ ! %# ! © ! ! ¨ ! %# !
¤
¡
$ !
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
function k1 = h * k2 = h * yh = y + function k1 = h * k2 = h * k3 = h * k4 = h * yh = y +
yh = rk2(f, x, y, h) feval(f, x, y); feval(f, x + h/2, y + k1/2); k2; yh = rk4(f, x, y, feval(f, x, y); feval(f, x + h/2, feval(f, x + h/2, feval(f, x + h, y (k1 + 2*k2 + 2*k3 h) y + k1/2); y + k2/2); + k3 ); + k4)/6.0;
Zadanie. Porówna´ przybli˙ one rozwiazania c z ˛ zagadnienia
f = inline(’x.*(1-y)’, ’x’, ’y’); h = 0.1; x = 0: h: 3; ye = 0;
Wprowadzenie do Matlaba, folia 46
¦
z rozwiazaniem dokładnym ˛ x = 0: 0.1: 3.
¦
6
'
¦
¨ ! ¡ ¤ ¤ ¡ # ¦ ¤ ¡ ¤ ¤ ¨ ¤
¡ ¡
dla
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
for i = x(1: end-1) ye = [ye euler(f, i, ye(end), h)]; end; y2 = 0; for i = x(1: end-1) y2 = [y2 rk2(f, i, y2(end),h)]; end; y4 = 0; for i = x(1: end-1) y4 = [y4 rk4(f, i, y4(end), h)]; end; delete(gcf) plot(x, ye, ’r+’) hold on plot(x, y2, ’b*’) plot(x, y4, ’mo’) sol = inline(’1-exp(-x.*x/2)’); plot(x, sol(x), ’y-’); legend(’Euler’, ’rk2’, ’rk4’, ... ’rozw. analityczne’, 0);
Wprowadzenie do Matlaba, folia 47
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
1
0.9
0.8
0.7
0.6
0.5
0.4 Euler rk2 rk4 0.2 rozw. analityczne
0.3
0.1
0
0
0.5
1
1.5
2
2.5
3
Zadanie. Powtórzy´ poprzednie zadanie dla c problemu
oraz x = 0: 0.3: 10.
Wprowadzenie do Matlaba, folia 48
¦
¡ ¦ ¤
¡
¡
# $! ¡ ¦ # ¦ ¡ ¤ ¤ ¨ ¤ ¤ ¡¡ ¨ ¤ ¦ ¦ ¤ ¡¡ ¤
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
f = inline(... ’[-y(1)/(x+eps)-y(2); y(1)]’, ... ’x’, ’y’); h = 0.3; x = 0: h: 10; ye = [0; 1]; for i = x(1: end-1) ye = [ye euler(f, i, ye(:,end), h)]; end; y2 = [0; 1]; for i = x(1: end-1) y2 = [y2 rk2(f, i, y2(:, end), h)]; end; y4 = [0; 1]; for i = x(1: end-1) y4 = [y4 rk4(f, i, y4(:, end), h)]; end; delete(gcf) plot(x, ye(2, :), ’r+’) hold on plot(x, y2(2, :), ’b*’) plot(x, y4(2, :), ’mo’) plot(x, besselj(0, x), ’y-’)
Wprowadzenie do Matlaba, folia 49
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
y45 = [0; 1]; [x y45]=ode45(f, x, y45); plot(x, y45(:, 2), ’g^’) legend(’Euler’, ’rk2’, ’rk4’,... ’rozw. analityczne’, ’ode45’, 0);
1
0.8
0.6
0.4
0.2
0
−0.2
−0.4 Euler −0.6 rk2 rk4 −0.8 rozw. analityczne ode45
−1
0
1
2
3
4
5
6
7
8
9
10
Wprowadzenie do Matlaba, folia 50
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Zadanie. W zbiorniku umieszczono moli i 1 mol benzenu. Oznaczmy odpowiednio przez , , , , i ilo´ci s , benzenu, chlorobenzenu, dwuchlorobenzenu i trójchlorobenzenu w chwili . Dla obj˛ to´ci produktu obowiazuja równania e s ˛ ˛
gdzie:
,
. Mo˙ na przyja´ z ˛c
.
Zwró´ my uwag˛ , ze ostatnie równanie ró˙ niczkowe c e ˙ z . nie jest istotne, bo zachodzi Zu˙ ywana ilo´c z s´ wynosi . Narysowa´ ewolucj˛ , , i w funkcji c e .
dla
Dla jakiej warto´ci otrzymuje si˛ maksimum s e chlorobenzenu? Jaka warto´c ˛ s´ na mol benzenu nale˙ y z wprowadzi´ do zbiornika, aby osiagna´ to c ˛ ˛c maksimum?
Wprowadzenie do Matlaba, folia 51
¦
! ¡ §
£
¦
¡ ¢
¡ ¡ ¥ § ¡ ££ ¦¡ # £ ¡¡ § ¡ ¨£ ¡ £ ¤
£
¥£ ¨ £ ¨ ¤ ¡ ! ¡ ¥ ¨ £ ¨ ¤ ¨ £
¦
£
¥
£
£
¤ £
©
¡ ¦
£
! ¡ £ ¡ ¦ ¦ © ¡ ©
# #
¦£ ¡¡ £ ¡ ¦
¦
¡¢
¡
¦
¡¢
§ ¡ £ § ¡ £
¦
¡¢ §
£¡
£
©
¡ ¡
¨£ ¦ ¤
¦
¥
£
¡
¤
i
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
%obliczenia opcje=odeset(’abstol’, 1.0e-4, ... ’reltol’, 1.0e-4); [t y]=ode45(’tricl’, ... [0:0.001:0.7], [1 0 0 ], opcje); q = y(:,1); r = y(:,2); s = y(:,3); t = 1-y(:,1)-y(:,2)-y(:,3); p = r+2*s+3*t; [m i] = max(r); %prezentacja delete(gcf); figure(’name’, ... ’Trojchlorowanie benzenu’, ... ’number’, ’off’, ’Menu’, ’none’); subplot(’Position’, ... [0.1 0.4 0.8 0.55]); plot(p, q, ’y-’,p, r, ’r:’,p, ... s, ’b-.’, p, t, ’g--’) legend(’benzen’, ’C_6H_5Cl’, ... ’C_6H_5CL_2’, ’C_6H_5CL_3’);
Wprowadzenie do Matlaba, folia 52
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
subplot(’Position’,... [0.1 0.01 0.8 0.35]); axis off st{1} = sprintf([’Maksimum ’, ... ’C_6H_5Cl wynosi %6.4f’, ... ’ i odpowiada’], m); st{2} = sprintf(... [’koncentracji C_6H_6’, ... ’ rownej %6.4f.’], q(i)); st{3} = sprintf([’Ilosc ’, ... ’chloru na mol C_6H_6’, ... ’ potrzebna do jego’, ... ’ otrzymania’]); st{4} = sprintf(... ’wynosi %6.4f.’, p(i)); text(0, 0.55, st, ’fontsize’, 14)
Wprowadzenie do Matlaba, folia 53
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
function yp = tricl(theta,y) k = [ 240 30 1]; p0 = 2; v = 1; p = p0 -y(2)-2*y(3)... -3*(1-y(1)-y(2)-y(3)); q = y(1); r = y(2); s = y(3); yp = [ -k(1)*p*q; (k(1)*q - k(2)*r)*p; (k(2)*r -k(3)*s)*p; ]/v;
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 benzen C6H5Cl C6H5CL2 C6H5CL3
Maksimum C6H5Cl wynosi 0.7430 i odpowiada koncentracji C6H6 rownej 0.0907. Ilosc chloru na mol C6H6 potrzebna do jego otrzymania wynosi 1.0762.
Wprowadzenie do Matlaba, folia 54
Matlab 1
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Tryb interpretera
2 * 3
ans = 6 2 * 3; pi
ans = 3.1416 sin(pi / 2) 1 [1 4; 6 8] 1 6 6: 2: 20 6 8 10 12 14 16 18 20 4 8
ans =
ans =
ans =
Wprowadzenie do Matlaba, folia 1
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
b
? Undefined function or variable ’b’ A = [1 4; 6 8]
A = 1 6 A - 1 4 8
ans = 0 5 A * A 3 7
ans = 25 54 sin(A)
36 88
ans = 0.8415 -0.2794
-0.7568 0.9894
Wprowadzenie do Matlaba, folia 2
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Programowanie
Zadanie 1. Napisa´ funkcj˛ znajdujaca c e ˛ ˛ pierwiastek równania liniowego. function x = rown1(a, b) if a == 0 if b ~= 0 x = []; else x = NaN; end; else x = -b / a; end; Uwaga: T˛ funkcj˛ zapisuje si˛ w pliku rown1.m ! e e e Efekt działania: rown1(1, 2) ans = -2
Wprowadzenie do Matlaba, folia 3
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Zadanie 2. Napisa´ funkcj˛ rozwiazujaca c e ˛ ˛ ˛ równania liniowe i kwadratowe. function x = rown21(a, b, c) if nargin == 2 x = rown1(a, b); elseif a == 0 x = rown1(b, c); else delta = b * b - 4 * a * c; if delta >= 0 if b > 0 x = (-b-sqrt(delta))/(2*a); else x = (-b+sqrt(delta))/(2*a); end; if x == 0 x = [x 0]; else x = [x (c/a)/x]; end; end; end;
Wprowadzenie do Matlaba, folia 4
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Wskazana jest równie˙ diagnostyka błedów: z if nargin < 2 error(’Za malo parametrow’); end;
P˛ tla for i wektoryzacja e
p˛ tla for e for i = 1: 100 y(i)=sin(i); end; Co wybiera´ ? Zawsze wektoryzacj˛ ! c e Zadanie. Napisa´ funkcj˛ obliczajaca sum˛ c e ˛ ˛ e p-tych pot˛ g n pierwszych liczb naturalnych. e Wersja „klasyczna”: function s = sump(n, p) s = 0; for i = 1: n s = s + i^p; end;
¡¡¡¡¡
wektoryzacja y = sin(1: 100);
Wprowadzenie do Matlaba, folia 5
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Wersja zwektoryzowana: function s = sumpv(n, p) s = sum((1: n).^p); Wektoryzacja jest mo˙ liwa, je˙ eli operacje w p˛ tli z z e nie zale˙ a od porzadku, w ktorym sa wykonywane. z˛ ˛ ˛
Specyfika p˛ tli for e
Zadanie. Co b˛ dzie rezultatem poni˙ szych e z instrukcji: for i=12:-2:1; disp(i), end; for i=12:-2:1; disp(i), ... i = 0; end; for i=12:-2:1; disp(i); ... if i == 8, break; end; end;
Wprowadzenie do Matlaba, folia 6
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Zadanie. Napisa´ funkcj˛ wyznaczajaca n-ty c e ˛ ˛ element ciagu Fibonacciego. ˛ Wersja prostsza: function f = fibs(n) if n == 1 | n == 2 f = 1; else a = 1; b = 1; for i = 3: n c = a + b; a = b; b = c; end; f = c; end; Wersja bardziej efektywna: function f = fib(n) persistent rez if length(rez) < 2 rez = [1 1]; end;
Wprowadzenie do Matlaba, folia 7
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
if length(rez) < max(n) for i = length(rez)+1: max(n) rez(i)=rez(i-1)+rez(i-2); end; end; f = rez(n); Efekt działania: fib([3, 9, 4]) ans = 2 34 3
Zadanie. Napisa´ funkcj˛ obliczajaca n c e ˛ ˛ pierwszych linii trójkata Pascala. ˛ function c = trojkat(n) c = zeros(n, n); c(:,1) = 1; for i = 2:n for j = 2:i c(i,j) = c(i-1,j-1)+c(i-1,j); end; end;
Wprowadzenie do Matlaba, folia 8
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Efekt działania: trojkat(6) ans = 1 1 1 1 1 1 0 1 2 3 4 5 0 0 1 3 6 10 0 0 0 1 4 10 0 0 0 0 1 5 0 0 0 0 0 1
Sposoby wektoryzacji
Prze´led´ my poni˙ sza sesj˛ : s z z ˛ e a = 1: 16; a = reshape(a, 4, 4) a = 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
Wprowadzenie do Matlaba, folia 9
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
a1 = a(:, 4: -1: 1)
a1 = 13 14 15 16 9 10 11 12 5 6 7 8 1 2 3 4
a2 = a(4: -1: 1, 4: -1: 1)
a2 = 16 15 14 13 a3 = a’ 12 11 10 9 8 7 6 5 4 3 2 1
a3 = 1 5 9 13 2 6 10 14 3 7 11 15 4 8 12 16
Wprowadzenie do Matlaba, folia 10
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
To samo wykonuje zreszta predefiniowana funkcja ˛ nchoosek. Zadanie. Utworzy´ wektor, w którym ka˙ demu c z indeksowi parzystemu p b˛ dzie odpowiada´ e c kwadrat, a indeksowi nieparzystemu – sze´cian s p-tej liczby naturalnej. v(1: 2: n) = (1: 2: n).^3; v(2: 2: n) = (2: 2: n).^2; Zadanie. Majac dany wektor x, utworzy´ tablic˛ ˛ c e a o n wierszach b˛ dacych kopiami x. e ˛ x = 1: 6; n = 5; a = x(ones(1, n), :) a = 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4 5 5 5 5 5 6 6 6 6 6
Wprowadzenie do Matlaba, folia 11
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Zadanie. Dla tablicy A = 3 -1 -1 1 1 2 0 3 0 -3 2 5 3 2 0 -4 2 -1 -4 1 3 2 5 5 -4 -1 0 -2 1 -5 1 2 5 3 2 -1
wygenerowa´ tablic˛ zło˙ ona z elementów le˙ acych c e z ˛ z˛ w wierszach od trzeciego do piatego i kolumnach 6, ˛ 1 i 3. B = A(3: 5, [6 1 3]) B = 5 3 2 -1 1 1 0 -4 2
Jak w tej tablicy wyzerowa´ wiersze nieparzyste, a c elementy wi˛ ksze od 3 zamieni´ na 1. e c A(1: 2: 6, :) = 0; A(A > 3) = 1;
Wprowadzenie do Matlaba, folia 12
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Zadanie. Narysowa´ wykres funkcji c
w przedziale nieciagło´ci w punkcie ˛ s
x = linspace(pi-1, pi+1, 50); y = 3*x.^2+log((x-pi).^2)/pi^4+1; plot(x, y)
55
50
45
40
35
30
25
20
15
10
2
¡ ¤ ' ¨ # (! &%$! " ! ¨ ¦ ¤ ©§¤¥£ ¢ ¨ ¦ ¡
2.5 3 3.5
. Dlaczego nie wida´ c ?
4
4.5
Wprowadzenie do Matlaba, folia 13
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Zadanie. Narysowa´ powierzchni˛ c e
dla
x = linspace(-1, 1, 5)
x = -1.00
-0.50
y = x; [X, Y] = meshgrid(x, y)
X = -1.00 -1.00 -1.00 -1.00 -1.00 Y = -1.00 -1.00 -1.00 -1.00 -1.00 -0.50 -0.50 -0.50 -0.50 -0.50 0 0 0 0 0 0.50 0.50 0.50 0.50 0.50 1.00 1.00 1.00 1.00 1.00 -0.50 -0.50 -0.50 -0.50 -0.50 0 0 0 0 0 0.50 0.50 0.50 0.50 0.50 1.00 1.00 1.00 1.00 1.00
¦ §¤ ¡ ¦
. 0 0.50 1.00
Wprowadzenie do Matlaba, folia 14
! £ ¤ ¡¢! ¡ #
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Z = -X.^2 - Y.^2
Z = -2.00 -1.25 -1.00 -1.25 -2.00 -1.25 -0.50 -0.25 -0.50 -1.25 -1.00 -0.25 0 -0.25 -1.00 -1.25 -0.50 -0.25 -0.50 -1.25 -2.00 -1.25 -1.00 -1.25 -2.00
surf(x, y, Z) colormap(gray)
0
−0.5
−1
−1.5
−2 1 0.5 0 0 −0.5 −1 −1 −0.5 0.5 1
Wprowadzenie do Matlaba, folia 15
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Zadanie. Napisa´ funkcj˛ v = vnd(x), która c e na bazie wektora x o składowych utworzy nast˛ pujaca macierz (tzw. macierz e ˛ ˛ Vandermonde’a):
function v = vnd(x) x = x(:); % x powinno byc kolumna v = ones(size(x)); for i = 1: length(x) - 1 v = [v(:, 1).*x v]; end; vnd([2 3 4 5])
ans = 8 27 64 125 4 9 16 25 2 3 4 5 1 1 1 1
Wprowadzenie do Matlaba, folia 16
!
£ ¤
¤ ¦£ ¤
¨
¨
¨
¨
¨
. . .
¨ £ £
. . .
¤¤%# ¡¢¢¡¢# ¤ £ ¡
. . . . . .
! !
©
¤ ©¦ ¤ ¤ ¦ ¤
¨ ¨
¨ ¨ ¨
¨
©£
£
¤ ¤ ¤
¦ ¦ § ¦ ¦ ¦
¥¦
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
P˛ tla while i rekurencja e
P˛ tla niesko´ czona z instrukcjami break: e n while 1 ... if condition_1 break; end; ... if condition_2 break; end end; Zadanie. Napisa´ funkcj˛ szukajaca w wektorze c e ˛ ˛ t ostatniego wystapienia zadanego elementu. ˛ function i = szuk1(t, x) i = 0; for j = length(t): -1: 1 if t(j) == x i = j;
Wprowadzenie do Matlaba, folia 17
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
break; end end; function i = szuk2(t, x) i = find(x == t); if isempty(i) i = 0; else i = i(end); end; function i = szuk3(t, x) i = length(t); while 1 if i == 0 break; end; if t(i) == x break; end; i = i - 1; end;
Wprowadzenie do Matlaba, folia 18
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
t = randperm(6000); tic; szuk1(t, 3000); toc
elapsed_time = 0.0400 tic; szuk2(t, 3000); toc
elapsed_time = 0.0110 tic; szuk3(t, 3000); toc
elapsed_time = 0.1000 Zadanie. Znale´ c pierwiastki równania z´
¡
Wprowadzenie do Matlaba, folia 19
Zastosowa´ algorytm: c 1. Znale´ c pierwiastek rzeczywisty z´ Newtona.
© ©
. (Łatwo 2. Podzieli´ trójmian przez c pokaza´ , ze odpowiedni iloraz ma posta´ c ˙ c , gdzie: , .) 3. Rozwiaza´ otrzymane równanie kwadratowe ˛ c otrzymujac pozostałe dwa pierwiastki. ˛
£
¤
¦
¨©¤ ¡
¡
¥
¤ ¤
¨ ¤ ¤ ¨ §¤ £ ¨ ¢¤ ¦ ¡
§ ¨ ¤
metoda ˛
¨ § ¡ ¤ ¨ § ¨ ¦ ¤ ¤©
¨
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
function x1 = newton1(a, b, c, d) epsi = 1.0e-5; x1 = 0; while 1 fprim = (3*a*x1+2*b)*x1+c; if fprim == 0 x1 = -b/a; else dx = -(((a*x1+b)*x1+c)*x1+d)... /fprim; x1 = x1 + dx; if abs(dx) < epsi break; end end end; function x = rown3(a, b, c, d) if a == 0 x = rown21(b, c, d); else x = newton1(a, b, c, d); beta = b + a * x;
Wprowadzenie do Matlaba, folia 20
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
gamma = c + beta * x; x = [x rown21(a, beta, gamma)]; end; rown3(2, -12, 22, -12)
ans = 1.0000
3.0000
2.0000
Uwaga! W praktyce stosuje si˛ funkcj˛ roots: e e roots([2, -12, 22, -12]) ans = 3.0000 2.0000 1.0000 Ogólny schemat rekurencji: function y = p(n) if n == 1 y = wyrazenie_explicite; else y = zwiazek_z( p(n - 1) ); end;
Wprowadzenie do Matlaba, folia 21
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Zadanie. Napisa´ rekurencyjna wersj˛ funkcji c ˛ e wyznaczajacej warto´c , ˛ s´ . function y = potega(n, p) if p == 0 y = 1; elseif rem(p, 2) == 1 y = n * potega(n*n, floor(p/2)); else y = potega(n*n, floor(p/2)); end; Zadanie. Napisa´ rekurencyjna wersj˛ funkcji c ˛ e szukajacej w wektorze t ostatniego wystapienia ˛ ˛ zadanego elementu. function i = szuk4(t, x) if isempty(t) i = 0; elseif t(end) == x i = length(t); else i = szuk4(t(1: end-1), x); end;
¦
§ ¦§¤# ¢ ¥ £ ¡
Wprowadzenie do Matlaba, folia 22
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Zadanie. Napisa´ funkcj˛ obliczajaca c e ˛ ˛ przybli˙ ona warto´c całek adaptacyjna metoda z ˛ s´ ˛ ˛ Simpsona. Wykorzysta´ ja do obliczenia całki c ˛
Rozwiazanie: Wzór Simpsona jest postaci ˛
mo˙ na obliczy´ na dwa sposoby: z c
gdzie:
´ – srodek
.
Wprowadzenie do Matlaba, folia 23
' © %$#¤%%!"¤ " ! #
Aby obliczy´ całk˛ z dokładno´cia c e s ˛ , dzielimy przedział na podprzedziały w taki sposób, aby na ka˙ dym z nich osiagna´ z ˛ ˛c dokładno´c s´ .
' © %$#¤%#%#¤ " ! ! ¥ © %$#%# ¦ 67© $%! ¤%# ©¨ ¨ ¦ 67© ! ¤
§
§
©£
¤ ¢ ! ¨ ¦ ¤ ¡
¥
¨
§ ¨ ¥ § # #
¨ ¤ ¤ ¦
#¤ ! ' £ # "
©
$ %!
$ %!
¥
©¨ §¤ ¢ ¦ ¡
¦ %$! ¤ 6 7© ¥ ¤ ¤ %# !#%# ©¨ ¡ ¥ ¤%# ! %# © ¨ ¡ ¤ ¥ ¤ ¢ ¤ ¡
©
¨
$! %#¤
¥
¤
¥
¡
! 4 ) ' 231(& 0) ' ! 54
£ ¤
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Mo˙ na udowodni´ , ze z c ˙
gdzie:
– dokładna warto´c odpowiedniej całki. s´
Wynika stad algorytm wyznaczania ˛ : Oblicza si˛ warto´ci e s
oraz „bład” ˛ . Je˙ eli „bład” jest wi˛ kszy z ˛ e ni˙ z , dokonuje si˛ podziału e na podprzedziały i ,a nast˛ pnie powtarza powy˙ sze obliczenia na ka˙ dym e z z z nich z osobna. Implementacja: Funkcj˛ podcałkowa mo˙ na e ˛ z zdefiniowa´ dwoma sposobami: c w pliku fun.m function y = fun(x) y = x.^4.*log(x+sqrt(x.*x+1));
' %$#¤%# ¦ %$! ¤ " ' ¦ 67© ! ' © %$#¤%%! © ¤ " 67© ! #
¥ %$#%# ¦ $%! ¤%# ! ¤ 67© © ¥ ¦ 7© $%! %# #¤%# ¥ ¤ ! 6 © $%#%# #¤%# ! ¤ !
Wprowadzenie do Matlaba, folia 24
! 54
! ! 4 !
£
%! $
¦ ! 54 © ¨ © ¨ © ¨
¤# %%! ¤ " #¤ %#¤ ! ! $! ¦ £ ¨ © £© ¡ ¦£
¡ ¡
© %$#%%#%# ! ¤#! ¤
©
! ¨
£ £
! ¨
¥
¡ ¢©¨
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
lub w linii polece´ : n fun = inline( ’x.^4.*log(x+sqrt(x.*x+1))’); Moduł inicjujacy rekurencj˛ , zapisany w pliku ˛ e simps.m , ma posta´ c function int = simps(f, przedz, epsi) int = []; a = przedz(1); b = przedz(2); if nargin < 3 epsi = 0.5*sqrt(eps)/(b-a); else epsi = 0.5*epsi/(b-a); end; przyb_pocz = (feval(f,a) ... +4*feval(f,0.5*(a+b)) ... + feval(f,b))*(b-a)/6.0; int = simps_rek(f, przedz, ... przyb_pocz, epsi); W tym samym pliku zapisuje si˛ funkcj˛ e e wewn˛ trzna simps_rek realizujaca rekurencj˛ . e ˛ ˛ ˛ e
¡¡¡¡¡ ¡¡¡
Wprowadzenie do Matlaba, folia 25
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
% funkcja wewnetrzna function intr = simps_rek(f, ... przedz, stare_przyb, epsi) a = przedz(1); b = przedz(2); srod = (b+a)*0.5; h = (srod-a)/6.0; fsrod = feval(f,srod); s1 = (feval(f,a) ... + 4.0*feval(f,(a+srod)*0.5) ... + fsrod)*h; s2 = (fsrod ... + 4.0*feval(f,(srod+b)*0.5) ... + feval(f,b))*h; err = abs(s1+s2-stare_przyb); if err>=15.0*epsi*(b-a) s1 = simps_rek(f, [a, srod], ... s1, epsi); s2 = simps_rek(f, [srod, b], ... s2, epsi); end; intr = s1 + s2;
Wprowadzenie do Matlaba, folia 26
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Zauwa˙ my, ze wywołanie simps zale˙ y od z ˙ z sposobu zdefiniowania funkcji podcałkowej: ¥ z zastosowaniem pliku: simps(’fun’, [0 2], 1.0e-4) ans = 8.1534 ¥ za zastosowaniem inline: simps(fun, [0 2], 1.0e-4) ans = 8.1534
Ciagi znaków i pliki tekstowe ˛
Zadanie. Wy´wietli´ wszystkie znaki ASCII o s c kodach od 32 do 255, oprócz znaku o kodzie 127. Dlaczego nie wymaga si˛ wy´wietlenia pozostałych e s znaków? i = char(32:255); i(127) = []; i
Wprowadzenie do Matlaba, folia 27
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Zadanie. Zinterpretowa´ poni˙ sze rezultaty: c z abs(’zorro’) ans = 122
111
114
114
111
sin(’zorro’)
ans = 0.4987 -0.8646 0.7850 0.7850 -0.8646 sin zorro
ans = 0.4987 -0.8646 0.7850 0.7850 -0.8646 Zadanie. Napisa´ funkcj˛ okre´lajaca liczb˛ c e s ˛ ˛ e wystapie´ danego znaku. ˛ n function lwyst = ile(s, x) lwyst = sum(s == x); Zadanie. Napisa´ funkcj˛ okre´lajaca czy dane c e s ˛ ˛ słowo jest palindromem. function r = palindrom(s) r = strcmp(s, s(end: -1: 1));
Wprowadzenie do Matlaba, folia 28
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Zadanie. Napisa´ funkcj˛ zamieniajaca pierwsze c e ˛ ˛ litery słów danego ciagu znaków na du˙ e litery. ˛ z function s1 = duze_lit(s) s = [’ ’ s]; ind = findstr(s, ’ ’); s = [s ’ ’]; s(ind+1) = upper(s(ind+1)); s1 = s(2: end-1); duze_lit(’to jest
ciag’)
ans = To Jest
Ciag
Zadanie. Napisa´ funkcj˛ odwracajaca kolejno´c c e ˛ ˛ s´ słów w zdaniu. function s1 = odwroc(s) [pslowo reszta] = strtok(s); if isempty(pslowo) s1 = s; else s1 = [odwroc(reszta), ’ ’* ... ones(1,(length(reszta)>0)), pslowo]; end;
Wprowadzenie do Matlaba, folia 29
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
odwroc(’to jest jedno zdanie’)
ans = zdanie jedno jest to Zadanie. Napisa´ funkcj˛ zapisujaca w pliku ciag c e ˛ ˛ ˛ znaków oraz odpowiednia funkcj˛ czytajaca ten ˛ e ˛ ˛ ciag. ˛ function zapisz(s, nazwa) f = fopen(nazwa,’w’); fwrite(f, s); fclose(f); function s = czytaj(nazwa) f = fopen(nazwa,’r’); s = fread(f); fclose(f); s = char(s’);
Wprowadzenie do Matlaba, folia 30
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Formatowane wyj´cie: s x = cumprod(1.0e-4*ones(5, 1)); a = [x, exp(x), exp(-x), ... (exp(x)- exp(-x))./(2*x)]; fprintf(’%15s %15s %15s %15s\n’, ... ’x’, ’exp(x)’, ’exp(-x)’, ’poch.’); fprintf( ... ’%15.4e %15.4e %15.4e %15.4e\n’, a’); Rezultat: x 1.0000e-004 1.0000e-008 1.0000e-012 1.0000e-016 1.0000e-020 exp(-x) 9.9990e-001 1.0000e+000 1.0000e+000 1.0000e+000 1.0000e+000 exp(x) 1.0001e+000 1.0000e+000 1.0000e+000 1.0000e+000 1.0000e+000 poch. 1.0000e+000 1.0000e+000 1.0000e+000 5.5511e-001 0.0000e+000
Wprowadzenie do Matlaba, folia 31
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Bł˛ dy zaokraglen e ˛ ´
Zadanie. Wyja´ni´ rezultat wykonania s c poni˙ szego skryptu: z x y z z = = = = 1.0e+29; 1.0e-9; ((y + x) - x) / y % => 0 (y + (x - x)) / y % => 1
Co to jest eps? eps
ans = 2.2204e-016 2^(-52) ans = 2.2204e-016 1 + eps == 1 ans = 0
1 + eps/2 == 1 ans = 1
Wprowadzenie do Matlaba, folia 32
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Zadanie. Napisa´ program obliczajacy warto´ci c ˛ s dla , -20, -10, 0, 10, 20, 30. Wykorzysta´ w tym celu definicj˛ c e
i porówna´ rezultaty z generowanymi przez exp. c function s = expo( x ) t = 1.0; n = 1; s = 0.0; while 1 if ( abs(t) <= 1.0e-14 ) break; end; s = s + t; t = t * (x/n); n = n + 1; end Skrypt: for i = -3:3 fprintf(’%4d %10.3e %10.3e \n’, ... 10*i, exp(10*i), expo(10*i)); end
¢ ¤ ¥ ££
£
¤
¡
¡
'
¦
£ ¡ ¤
'
Wprowadzenie do Matlaba, folia 33
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
-30 -20 -10 0 10 20 30
9.358e-014 2.061e-009 4.540e-005 1.000e+000 2.203e+004 4.852e+008 1.069e+013
6.103e-006 6.148e-009 4.540e-005 1.000e+000 2.203e+004 4.852e+008 1.069e+013
Pytanie: Jak ograniczy´ bł˛ dy zaokragle´ ? c e ˛ n Odpowied´ : Ograniczy´ obliczenia do przedziału z c :
Wprowadzenie do Matlaba, folia 34
¤ © ¤ ¤ ! ¨ ¡ ¥¦ ¦ ¤ ! ¨ ¤ ! ¡ ¤ ¨ ¤
Zadanie. Porówna´ warto´ci funkcji c s w przedziale z warto´ciami tej samej funkcji, ale liczonymi wg s formuły
' (!
" #
¦ ¦ ¦¡
¦ ¦
!#
¨
!!!
¤
¤ #
¡
¥ ©
¦"
'
¦
' ' ' !
© ¥ § ¦ ¤ ¥
#
¢£¡
¤ ! ¡ ¤
¡
' ¥ ¥
' (! # ¦ "
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
function e = bin1(x) x = 1-x; e = x.*x.*x.*x.*x.*x; function e = bin2(x) x2 = x.*x; x4 = x2.*x2; e = 1 - 6.*x + 15.*x2 ... - 20.*x2.*x + 15.*x4 ... - 6.*x4.*x + x4.*x2; Skrypt: poc = 99990; kon = 100010; for i=(poc: kon)/100000 fprintf(’%15.5f %10.3e %10.3e\n’, ... i, bin1(i), bin2(i)); end; i=(poc: kon)/100000; plot(i, bin1(i)); pause; plot(i, bin2(i));
Wprowadzenie do Matlaba, folia 35
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
bin1.m:
1 x 10
−24
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0 0.9999
1
1
1
1.0001
bin2.m:
4 x 10
−15
2
0
−2
−4 0.9999
1
1
1
1.0001
Wprowadzenie do Matlaba, folia 36
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Ró˙ niczkowanie numeryczne z
Dla pierwszych pochodnych stosuje si˛ e
lub
Dla drugich pochodnych mamy:
for x = cumprod(0.1 * ones(1, 10)) fprintf(’%2.0e %12.10f %14.10f \n’, ... x, (exp(x) - exp(-x))/(2 * x), ... (exp(x) - 2 + exp(-x))/x^2); end;
¨
Wprowadzenie do Matlaba, folia 37
¡
¦ ¥ ¨
¡
¥
¨
¡
¥
¡
¡
¡
¥ ! ¦ ! ¡ ¥ ¥ ¨ ¦
¡
¥
© § ¤¥ ¨
¦ ! ¥ ¡ ¡ ' ¨ # " ¦
¥
¡ ¡
! ¡ ¥ ¡ ' ¨ # " ¦
¥ £ ¦ ¤ ¡ ¡ ¤ ¥ £ ¡
¥
¡
¥
¡
¥
£
¢
¥
¥ ¥
¥ ¥
¢ ¢
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Efekty sa intrygujace: ˛ ˛ 1e-001 1e-002 1e-003 1e-004 1e-005 1e-006 1e-007 1e-008 1e-009 1e-010 1.0016675002 1.0008336112 1.0000166667 1.0000083334 1.0000001667 1.0000000834 1.0000000017 1.0000000050 1.0000000000 0.9999989725 1.0000000000 0.9999778783 0.9999999995 0.9992007222 0.9999999939 0.0000000000 1.0000000272 111.0223024625 1.0000000827 0.0000000000
Wyja´nienie: s
¦ ¦
Wprowadzenie do Matlaba, folia 38
¢
¢ ¡
$
¡ ¢ ¡ ¡
¦ ! ¥ ¡ ¨ ¥ ¥ ¨ ¦ ! ¥ £¤¡¢ ¥ ¥ ¨ ¨ ¡
¢ ¡ ¤¥ ¦£ ¤¥ ¦£ ¡ ¡
$
¥
¡
¡
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Całkowanie numeryczne
Zadanie. Z zastosowaniem funkcji quad lub quad8 obliczy´ całk˛ c e
f1 = inline(... ’x.^4.*log(x+sqrt(x.*x+1))’); quad8(f1, 0, 2) ans = 8.1534 Zadanie. Z zastosowaniem funkcji dblquad obliczy´ całk˛ c e
function v = f2(x, y) v = exp(-x.*x-y.*y);
¤ ¢ ! ¨ ¦ ¤ ¡ ¡ ¤
¦
¡ ¢ ' ¡
¨ ¤ ¤ ¦
¦ "
"
Wprowadzenie do Matlaba, folia 39
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
dblquad(’f2’, 0, 5, 0, 5, ... [], [], ’quad8’) ans = 0.7854 Zadanie. Narysowa´ wykres funkcji c
w przedziale
.
f = inline(’cos(x.*sin (u))./pi’, ... ’u’, ’x’); x = linspace(0, 10, 100); for i = 1: 100 y(i) = quad(f, 0, pi, ... [], [], x(i)); end; plot(x, y);
¡¢ ¤ §¥£ ¦ ¤
¨©¦
¢
' ¦ ! # ¦"
!
¡
¤
Wprowadzenie do Matlaba, folia 40
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
1.2
1
0.8
0.6
0.4
0.2
0
−0.2
−0.4
−0.6
0
1
2
3
4
5
6
7
8
9
10
Metoda Monte-Carlo
– ciag niezale˙ nych zmiennych losowych o ˛ z jednakowym rozkładzie, ,
§¨
jednostajny na hiperkostce
¡ ¡ ¡ ¡
Wprowadzenie do Matlaba, folia 41
!
¡ ¡¢¢¡ ¡
§¨
¡ !"(! # ¦ " '
©
¡
¡ ¥
¥ ¡ ¢ ! # ¢¢¡¢# ¡ ¡ ©
¦
¨
¡
¡
¡
¡
¥
¡
©
©
¢
¨ ¨ ¢¢¨
¡
!
©
©
¥
©
¢
¡
£ ¤¡
¡
!
©
¡ ¢
¡
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Bład mo˙ na oszacowa´ za pomoca odchylenia ˛ z c ˛ standardowego z próby (zob. std) pomno˙ onego przez z . Zadanie. Napisa´ funkcj˛ generujaca c e ˛ ˛ macierz, której ka˙ dy wiersz zawiera prób˛ z z e rozkładu jednostajnego na danym przedziale.
function x = random(gran, n) x = rand(size(gran,1), n); dlug = gran(:, 2) - gran(:, 1); x = gran(:, ones(1, n)) ... + x.*dlug(:, ones(1, n));
Wprowadzenie do Matlaba, folia 42
!
¡
¡¢¡¢¡ ©
§ ¨
¤ # ¢¡¢¢# © ¤ ¡ ¡
¡ ¢ !
¥
¡¢¡¢¡¢ ©
©
#
© ¤
#
¡
¥
¦
¡ ¡ ¡
¥
¤ ¥¢
¤ £
¡
o hiperobj˛ to´ci e s
¡
¢!
jednostajny na hiperkostce
§ ¨
' ! £ #%! " ©
!
¡
¨¢¢¨ ¨
2
2
¢
¥
£
©
¢
!
¡
¡
¡
!
©
¡
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Zadanie. Napisa´ funkcj˛ przybli˙ ajaca całk˛ c e z ˛ ˛ e wielokrotna po zadanej hiperkostce. ˛ function [int, st] = mtc( ... fun, gran, npkt); if nargin <= 2 npkt = 10000; end; obj = prod(gran(:,2)-gran(:,1)); z = feval(fun, random(gran, npkt)); if nargout == 2 st = obj * std(z) / sqrt(npkt); end; int = obj * mean(z); Zadanie. Obliczy´ obj˛ to´c przeci˛ cia kul o c e s´ e ´ srodkach w i oraz promieniach odpowiednio 3 i 2. function y = ff(x) y = ((x(1,:)).^2 + x(2,:).^2 ... + x(3,:).^2 < 9)... & ((x(1,:)-2).^2 + x(2,:).^2 ... + x(3,:).^2 < 4);
¦# ¦# ¦# ¦# ¦
Wprowadzenie do Matlaba, folia 43
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
[int, st] = mtc(’ff’, ... [0 3; -2 2; -2 2], 100000) int = 24.6931 st = 0.0759
Równania ró˙ niczkowe z
Metoda Eulera
¡
gdzie:
Wprowadzenie do Matlaba, folia 44
¦
¡
¢
¨ ! %# ! ¤
¡
gdzie:
,
,
£
¡
©
£ %# ¡¢¢¡¢# © ¡¢ ¡ $ £ ¥ ¤ ¡ # ¤ # %¤ ¡
¥
¥
¡
¨ !¤ ¡ ©$ !¤ ¨ ! ¡ © $ !
¡ ¡
¥
¤
"
¡
function yh = euler(f, x, y, h) yh = y + h * feval(f, x, y);
Zadanie. Zaimplementowa´ wymienione metody. c
¡
¢
¡
¡
¡ ¨ ¦ ¨ ! ¡ ¥ ¦© © ! ! ! ¨ %# ! ¨ ¤ ¥ ¡ ! %# ! ¤ ¡ ©
¡
Metoda RK4
Metoda RK2
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Wprowadzenie do Matlaba, folia 45
¢
¡
¡
¡
¡
¡
¨ ! ¥ ¤ ¥ ¤ ¥ ¤
¡
¥
¡ © $ ! ¡ ¡ ¡ ¡ ¦ ¡
©
¨ ¨ ¡ ¨ ¦ ¨ © ¡ ! ¡ ¨ ! %# ¨ ! ¦ ! ¨ ! %# ! ¨ ! ¨ ! %# ! © ! ! ¨ ! %# !
¤
¡
$ !
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
function k1 = h * k2 = h * yh = y + function k1 = h * k2 = h * k3 = h * k4 = h * yh = y +
yh = rk2(f, x, y, h) feval(f, x, y); feval(f, x + h/2, y + k1/2); k2; yh = rk4(f, x, y, feval(f, x, y); feval(f, x + h/2, feval(f, x + h/2, feval(f, x + h, y (k1 + 2*k2 + 2*k3 h) y + k1/2); y + k2/2); + k3 ); + k4)/6.0;
Zadanie. Porówna´ przybli˙ one rozwiazania c z ˛ zagadnienia
f = inline(’x.*(1-y)’, ’x’, ’y’); h = 0.1; x = 0: h: 3; ye = 0;
Wprowadzenie do Matlaba, folia 46
¦
z rozwiazaniem dokładnym ˛ x = 0: 0.1: 3.
¦
6
'
¦
¨ ! ¡ ¤ ¤ ¡ # ¦ ¤ ¡ ¤ ¤ ¨ ¤
¡ ¡
dla
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
for i = x(1: end-1) ye = [ye euler(f, i, ye(end), h)]; end; y2 = 0; for i = x(1: end-1) y2 = [y2 rk2(f, i, y2(end),h)]; end; y4 = 0; for i = x(1: end-1) y4 = [y4 rk4(f, i, y4(end), h)]; end; delete(gcf) plot(x, ye, ’r+’) hold on plot(x, y2, ’b*’) plot(x, y4, ’mo’) sol = inline(’1-exp(-x.*x/2)’); plot(x, sol(x), ’y-’); legend(’Euler’, ’rk2’, ’rk4’, ... ’rozw. analityczne’, 0);
Wprowadzenie do Matlaba, folia 47
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
1
0.9
0.8
0.7
0.6
0.5
0.4 Euler rk2 rk4 0.2 rozw. analityczne
0.3
0.1
0
0
0.5
1
1.5
2
2.5
3
Zadanie. Powtórzy´ poprzednie zadanie dla c problemu
oraz x = 0: 0.3: 10.
Wprowadzenie do Matlaba, folia 48
¦
¡ ¦ ¤
¡
¡
# $! ¡ ¦ # ¦ ¡ ¤ ¤ ¨ ¤ ¤ ¡¡ ¨ ¤ ¦ ¦ ¤ ¡¡ ¤
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
f = inline(... ’[-y(1)/(x+eps)-y(2); y(1)]’, ... ’x’, ’y’); h = 0.3; x = 0: h: 10; ye = [0; 1]; for i = x(1: end-1) ye = [ye euler(f, i, ye(:,end), h)]; end; y2 = [0; 1]; for i = x(1: end-1) y2 = [y2 rk2(f, i, y2(:, end), h)]; end; y4 = [0; 1]; for i = x(1: end-1) y4 = [y4 rk4(f, i, y4(:, end), h)]; end; delete(gcf) plot(x, ye(2, :), ’r+’) hold on plot(x, y2(2, :), ’b*’) plot(x, y4(2, :), ’mo’) plot(x, besselj(0, x), ’y-’)
Wprowadzenie do Matlaba, folia 49
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
y45 = [0; 1]; [x y45]=ode45(f, x, y45); plot(x, y45(:, 2), ’g^’) legend(’Euler’, ’rk2’, ’rk4’,... ’rozw. analityczne’, ’ode45’, 0);
1
0.8
0.6
0.4
0.2
0
−0.2
−0.4 Euler −0.6 rk2 rk4 −0.8 rozw. analityczne ode45
−1
0
1
2
3
4
5
6
7
8
9
10
Wprowadzenie do Matlaba, folia 50
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
Zadanie. W zbiorniku umieszczono moli i 1 mol benzenu. Oznaczmy odpowiednio przez , , , , i ilo´ci s , benzenu, chlorobenzenu, dwuchlorobenzenu i trójchlorobenzenu w chwili . Dla obj˛ to´ci produktu obowiazuja równania e s ˛ ˛
gdzie:
,
. Mo˙ na przyja´ z ˛c
.
Zwró´ my uwag˛ , ze ostatnie równanie ró˙ niczkowe c e ˙ z . nie jest istotne, bo zachodzi Zu˙ ywana ilo´c z s´ wynosi . Narysowa´ ewolucj˛ , , i w funkcji c e .
dla
Dla jakiej warto´ci otrzymuje si˛ maksimum s e chlorobenzenu? Jaka warto´c ˛ s´ na mol benzenu nale˙ y z wprowadzi´ do zbiornika, aby osiagna´ to c ˛ ˛c maksimum?
Wprowadzenie do Matlaba, folia 51
¦
! ¡ §
£
¦
¡ ¢
¡ ¡ ¥ § ¡ ££ ¦¡ # £ ¡¡ § ¡ ¨£ ¡ £ ¤
£
¥£ ¨ £ ¨ ¤ ¡ ! ¡ ¥ ¨ £ ¨ ¤ ¨ £
¦
£
¥
£
£
¤ £
©
¡ ¦
£
! ¡ £ ¡ ¦ ¦ © ¡ ©
# #
¦£ ¡¡ £ ¡ ¦
¦
¡¢
¡
¦
¡¢
§ ¡ £ § ¡ £
¦
¡¢ §
£¡
£
©
¡ ¡
¨£ ¦ ¤
¦
¥
£
¡
¤
i
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
%obliczenia opcje=odeset(’abstol’, 1.0e-4, ... ’reltol’, 1.0e-4); [t y]=ode45(’tricl’, ... [0:0.001:0.7], [1 0 0 ], opcje); q = y(:,1); r = y(:,2); s = y(:,3); t = 1-y(:,1)-y(:,2)-y(:,3); p = r+2*s+3*t; [m i] = max(r); %prezentacja delete(gcf); figure(’name’, ... ’Trojchlorowanie benzenu’, ... ’number’, ’off’, ’Menu’, ’none’); subplot(’Position’, ... [0.1 0.4 0.8 0.55]); plot(p, q, ’y-’,p, r, ’r:’,p, ... s, ’b-.’, p, t, ’g--’) legend(’benzen’, ’C_6H_5Cl’, ... ’C_6H_5CL_2’, ’C_6H_5CL_3’);
Wprowadzenie do Matlaba, folia 52
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
subplot(’Position’,... [0.1 0.01 0.8 0.35]); axis off st{1} = sprintf([’Maksimum ’, ... ’C_6H_5Cl wynosi %6.4f’, ... ’ i odpowiada’], m); st{2} = sprintf(... [’koncentracji C_6H_6’, ... ’ rownej %6.4f.’], q(i)); st{3} = sprintf([’Ilosc ’, ... ’chloru na mol C_6H_6’, ... ’ potrzebna do jego’, ... ’ otrzymania’]); st{4} = sprintf(... ’wynosi %6.4f.’, p(i)); text(0, 0.55, st, ’fontsize’, 14)
Wprowadzenie do Matlaba, folia 53
Instytut Sterowania i Systemów Informatycznych, Politechnika Zielonogórska
function yp = tricl(theta,y) k = [ 240 30 1]; p0 = 2; v = 1; p = p0 -y(2)-2*y(3)... -3*(1-y(1)-y(2)-y(3)); q = y(1); r = y(2); s = y(3); yp = [ -k(1)*p*q; (k(1)*q - k(2)*r)*p; (k(2)*r -k(3)*s)*p; ]/v;
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 benzen C6H5Cl C6H5CL2 C6H5CL3
Maksimum C6H5Cl wynosi 0.7430 i odpowiada koncentracji C6H6 rownej 0.0907. Ilosc chloru na mol C6H6 potrzebna do jego otrzymania wynosi 1.0762.
Wprowadzenie do Matlaba, folia 54
