Wykład 6 Dr Aneta Polewko-Klim

Wykład 6 Dr Aneta Polewko-Klim https://play.google.com

Obliczenia numeryczne. Wykonując obliczenia numeryczne zakładamy: operujemy na liczbach, dopuszczając możliwość występowania niedokładności w ich reprezentacji, podczas obliczeń mogą występować pewne niedokładności.

Miejsca zerowe i minima funkcji.

Opcje Wszystkie te funkcje można wywołać z opcjonalnym parametrem opcje, np.: x1 = fzero(‘cos(x) - 0.2’, 3, opcje) Opcje określają parametry wywołania tych funkcji. Ustawienia tych parametrów można zmienić poprzez funkcję optimset. opcje = optimset(parametr, wartość, ... ) Argumenty funkcji optimset to pary parametr – wartość, Parametr jest łańcuchem znakowym.

Funkcja optimset - lista parametrów i ich możliwych wartości (kilka):

Wywołanie z nazwą funkcji (wyświetla parametry dla konkretnej funkcji)

Przykłady parametrów funkcji optimset:

Miejsca zerowe Przykład: >> fzero('cos(2*x)-.2',8) ans = 8.7401 >> cos(2*ans)-.2 ans = 2.8866e-015

>> op=optimset('Diagnostics','on',... 'Display','Iter'); Looking for a zero in the interval [7.0949, 8.9051] 12 7.91238 -1.19319 interpolation 13 8.70207 -0.0749546 interpolation 14 8.74193 0.00367162 interpolation 15 8.74007 2.49516e-005 interpolation 16 8.74006 -1.30267e-009 interpolation 17 8.74006 2.88658e-015 interpolation 18 8.74006 -4.05231e-015 interpolation Zero found in the interval: [7.0949, 8.9051]. ans = 8.7401 >> fzero('cos(2*x)-.2',8,op) Func-count x f(x) Procedure 1 8 -1.15766 initial 2 7.77373 -1.18715 search 3 8.22627 -0.935369 search 4 7.68 -1.14007 search 5 8.32 -0.7962 search 6 7.54745 -1.01789 search 7 8.45255 -0.565028 search 8 7.36 -0.750391 search 9 8.64 -0.19876 search 10 7.0949 -0.252615 search 11 8.9051 0.30677 search

Minimum – funkcja 1 zmiennej Przykład >> funkcja=inline('sin(2*x)+cos(x)/2'); >>fplot(funkcja,[-2*pi, 2*pi]) >> x1=fminbnd(funkcja,-2,4) x1 =2.4375 >> funkcja(x1) ans = -1.3679

Minimum – funkcja 2 zmiennej Do znalezienia minimum musimy Podać postać tej funkcji lub funkcję inline oraz punkt startowy. Jak widać na rysunku funkcja ta ma dwa minima. W zależności od tego jaki punkt startowy podamy, znajdziemy albo pierwsze albo drugie. Przykład >> fun=inline('-exp(-(x(1)-2).^2-x(2).^2)-exp(-(x(1)+2).^2-x(2).^2)'); >> x0=fminsearch(fun,[2 2]) x0 = 2.0000 -0.0000 >> fx0=fun(x0) fx0 = -1.0000

Pierwiastki wielomianów. Ogólna postać wielomianu: W(x) = a1xn + a2xn-1 + ... + anx + an+1 współczynniki wielomianu a1, a2, ... są uszeregowane według malejących potęg zmiennej x.

Przykład: W(x) = x2 + 3x – 4 >> roots([1 3 -4]) ans = -4 1 >> poly(ans) ans = 1 3 -4 W(x) = 3x5 – 2x4 + 4x3 – 6x + 1 >> roots([3 -2 4 0 -6 1]) ans = 0.2084 + 1.4463i 0.2084 - 1.4463i -0.9198 1.0000 0.1697

Układy równań liniowych Zapis macierzowy, każdy układ równań możemy przedstawić: A*X = B lub X*A = B Matlab dysponuje dwoma operatorami dzielenia, które odpowiednio rozwiązują przypadki: A*X = B mamy X = A\B –musi być równa liczba wierszy A i B X*A = B mamy X = B/A – musi być równa liczba kolumn A i B W technice częściej spotyka się pierwszy z tych przypadków. Macierz A w ogólnym przypadku może mieć wwierszy i k kolumn. Istnieją więc trzy możliwości: w = k – macierz kwadratowa, jeśli jest nieosobliwa X można wyznaczyć dokładnie; w > k – układ jest nadokreślony; w < k – układ jest nieoznaczony (rozwiązanie zawiera najwyżej w niezerowych elementów;

Układ kwadratowy. Przykład >> A=[2 -3 -4; -1 2 1 ; 3 -2 -1] A = 2 -3 -4 -1 2 1 3 -2 -1 >> B=[-5 -1 5]' B = -5 -1 5 >> X=A\B X = 2.0000 -1.0000 3.0000 Mamy układ trzech równań z trzema niewiadomymi. Musimy zdefiniować macierze A i B, a następnie obliczyć wynik.

Układ nieoznaczony cz.1 Przykład >> A=[6 8 7 3;3 5 4 1] A = 6 8 7 3 3 5 4 1 >> B=[1 2]' B = 1 2 >> X=A\B X = 0 0.7143 0 -1.5714 Liczba równań jest mniejsza od liczby niewiadomych

Układ nieoznaczony cz.2 Ogólne rozwiązanie takiego układu otrzymamy dodając do otrzymanego rozwiązania iloczyn macierzy Z będącej jądrem macierzy A i dowolnego wektora q. Jądro: Z = null(A) A*Z = 0 Przykład cd. >> Z=null(A) Z = -0.4253 -0.6519 -0.3981 0.3942 0.8126 -0.1607 0.0163 0.6276 >> q=[2 17]' q =2 17 >> X1=X+Z*q X1 = -11.9324 6.6186 -1.1066 9.1306

Układ nadokreślony. Gdy chcemy otrzymać n współczynników funkcji opisującej przebieg jakiegoś zjawiska. Dokonujemy wtedy np. dużo więcej niż n pomiarów i na ich podstawie szukamy wyniku. Przykład >> A=[1 -1 2;1 -2 -1;3 -1 5;-2 2 3] A =1 -1 2 1 -2 -1 3 -1 5 -2 2 3 >> B=[1 2 3 -4]' B = 1 2 3 -4 >> X=A\B X = 1.4286 -0.1429 -0.2857

Interpolacja. Zagadnienie dotyczy przechodzenia funkcji przez zadane punkty, tzw. węzły interpolacji. Standardowe procedury Matlaba realizują interpolację za pomocą wielomianów pierwszego i trzeciego stopnia, oraz funkcjami sklejanymi trzeciego stopnia. Funkcje interpolacyjne:

Metody interpolacji • ‘nearest’— interpolacja najbliższy sąsiad • ‘linear’ – interpolacja liniowa; • ‘cubic’ – interpolacja wielomianem trzeciego stopnia; • ‘spline’ – interpolacja funkcjami sklejanymi trzeciego stopnia;

Przykład: x=0:.4:2; y=exp(-.5*x).*cos(5*x); xi=0:.05:2; y1=interp1(x,y,xi,'nearest'); y2=inter1(x,y,xi,'linear'); y3=interp1(x,y,xi,'cubic'); y4=interp1(x,y,xi,'spline'); plot(x,y,'*',xi,exp(-.5*xi).*cos(5*xi),... xi,y1,':',xi,y2,':',xi,y3,':',xi,y4,':') legend('wezly','funkcja','nearest',… 'linear','cubic','spline')

Aproksymacja. • Zagadnienie aproksymacji polega na znalezieniu pewnej funkcji, której wykres z pewnym przybliżeniem przechodzi w pobliżu pewnych punktów. • Sama funkcja minimalizuje wartość pewnego kryterium dopasowania, np. aproksymacja średniokwadratowa, czyli minimalizujemy błąd średniokwadratowy. • Standardową metodą aproksymacji w Matlabie jest aproksymacja wielomianem: polyfit(x,y,n) • polecenie to znajduje wektor współczynników wielomianu aproksymującego dane zawarte w wektorach x i y. Wartość wielomianu możemy obliczyć funkcją polyval.

Aproksymacja. Przykład 1: x=1:8; y=[1.45 1.95 2.43 3.02 3.53 4 4.4 5.1]; a=polyfit(x,y,1) xx=1:.01:8; yy=polyval(a,xx); plot(x,y,'o,xx,yy) Otrzymaliśmy funkcję y=0.5121x+0.9304

Przykład 2: x=0:.4:2*pi; y1=sin(x); blad=rand(1,length(x))/10; y=y1+blad; a1=polyfit(x,y,2) a2=polyfit(x,y,3) xx=0:.01:2*pi; yy1=polyval(a1,xx); yy2=polyval(a2,xx); plot(x,y,'o',xx,sin(xx),':',xx,yy1,xx,yy2) Dla wielomianu drugiego stopnia mamy: y = – 0.0323x2 – 0.0982x + 0.7252 Dla wielomianu trzeciego stopnia mamy: y = 0.0921x3 – 0.8611 x2 + 1.8274x – 0.0793

Pochodna funkcji • Do znajdowania pochodnej możemy posłużyć się funkcją: • diff • Zwraca ona wektor różnic sąsiednich elementów wektora będącego jej argumentem. • Czyli pochodna funkcji danej w postaci tablicy (np. pomiary), obliczymy jako: • diff(y)./diff(x). • Należy przy tym pamiętać, że liczba elementów wektora zwracanego przez funkcję diff, jest o jeden mniejsza, niż liczba elementów argumentu.

Pochodna funkcji Przykład x=72:.5:90; y=[7237 7224 7203 7183 7164 7149 7129 7107 7060 7011 6912 6746 6517 ... 6225 5875 5449 4950 4376 3789 3202 2614 2082 1639 1290 1000 ... 758 581 470 394 328 284 241 191 146 103 52 0]; yprim=diff(y)./diff(x); subplot(2,1,1),plot(x,y,'+'),grid on ylabel('Dane') subplot(2,1,2),plot(x, [yprim 0],'o'),grid on ylabel('Pochodna')

Pochodna wielomianu • Jeśli mamy dany wielomian, do obliczenia wartości pochodnej możemy użyć funkcji: • polyder • Zwraca ona współczynniki wielomianu pochodnego. Wywołanie: • a_prim = polyder(a) • [L_prim, M_prim] = polyder(L,M)

Pochodna wielomianu a=[-1 3 5 -2]; x=0:.05:3; y=polyval(a,x); a1=polyder(a); y1=polyval(a1,x); subplot(2,1,1); plot(x,y) grid on ylabel('Dane') subplot(2,1,2) plot(x,y1) grid on ylabel('Pochodna')

Pochodna wielomianu - funkcja wymierna Dodatkowo należy musimy zdefiniować licznik i mianownik: Przykład: L=[1 -3]; M=[1 5 -2 1]; x=-3:.05:3; y=polyval(L,x)./polyval(M,x); [L1,M1]=polyder(L,M); y1=polyval(L1,x)./polyval(M1,x); subplot(2,1,1),plot(x,y),grid on ylabel('Dane') subplot(2,1,2),plot(x,y1),grid on ylabel('Pochodna')

Całkowanie – całki oznaczone. Parametry: f – łańcuch znaków określający nazwę funkcji; a, b – przedział całkowania; tol – tolerancja błędu (domyślnie 10-3); tr – jeśli jest różny od 0, rysowany jest wykres obrazujący postęp w obliczeniach lub wypisywane kolejne iteracje.

Całkowanie Przykład: >> quad('sin(x)',1,10) ans = 1.3794 >> quad8('sin(x)',1,10) ans = 1.3794 Gdy dodamy kolejne parametry: >> quad('sin(x)',1,10,1e-3,1) 5 1.0000000000 9.00000000e+000 2.1980581837 7 1.0000000000 4.50000000e+000 -0.1666251940 9 1.0000000000 2.25000000e+000 1.5343093586 11 1.0000000000 1.12500000e+000 1.0665674753 13 2.1250000000 1.12500000e+000 0.4678628302 15 3.2500000000 2.25000000e+000 -1.7026633720 17 3.2500000000 1.12500000e+000 -0.6631045444 19 4.3750000000 1.12500000e+000 -1.0396930457 21 5.5000000000 4.50000000e+000 1.5317252037 23 5.5000000000 2.25000000e+000 0.6048270787 25 7.7500000000 2.25000000e+000 0.9427905376 ans = 1.3793

Całkowanie Przykład: W przypadku funkcji quad8 otrzymamy wykres: >> quad8('sin(x)',1,10,1e-3,1) ans = 1.3794

Równania różniczkowe Równania różniczkowe możemy podzielić na 3 grupy: - równania różniczkowe zwyczajne gdzie szukamy rozwiązania równania różniczkowego dla zadanego warunku początkowego; - równania różniczkowe zwyczajne gdzie szukamy rozwiązania równania dla zadanych warunków granicznych; - równania różniczkowe cząstkowe

Równania różniczkowe I rzędu W rozwiązaniach wielu problemów fizycznych, ekonomicznych wymagana jest znajomość funkcji y=y(t) przy znajomości funkcji y'=f(t, y) oraz warunków początkowych y(a)=y0, gdzie: a oraz y0 są liczbami rzeczywistymi f funkcją w postaci jawnej. W takim przypadku mamy do czynienia z równaniem różniczkowym zwyczajnym pierwszego rzędu

Przykład: Szukamy rozwiązania numerycznych y = y(t) dla wartości t = 0, .25, .5, .75, 1 dla równania różniczkowego y' = -2ty2, przy warunku początkowym y(0) = 1. dy = @(t,y)(-2*t.*y(1).^2); tspan = [0 .25 .5 .75 1]; y0 = 1; [t1 y1] = ode23('eq1', tspan, y0); [t2 y2] = ode45('eq1', tspan, y0); [t1 y1 y2] % porównanie wyników

Wykład 6 Dr Aneta Polewko-Klim