Matlab Aproksymacja Interpolacja Inne metody obliczeniowe

1. MOiPP Wykład 6 • Matlab • Aproksymacja • Interpolacja • Inne metody obliczeniowe

2. Aproksymacja Metoda aproksymacji polega na znalezieniu funkcji f(x), której wykres przechodzi w pobliżu zbioru zadanych punktów, zaś sama funkcja minimalizuje wartość pewnego kryterium dopasowania Wykorzystanie np. gdy mamy wiele dyskretnych punktów pomiarowych i chcemy znaleźć funkcję aproksymującą. Funkcja aproksymująca nie musi wcale przyjmować identycznych wartości jak funkcja aproksymowana. Kryteria aproksymacyjne są różne, np. kryterium średniokwadratowe jak w metodzie najmniejszych kwadratów.

3. Niech w naszym przypadku funkcją aproksymującą f(x) będzie wielomian n-tego stopnia w(x). Do znalezienia współczynników wielomianu używamy w Matlabie funkcji polyfit, której składnia ma postać: p=polyfit(x,y,n) gdzie: x - jest zbiorem N posortowanych rosnąco wartości współrzędnych zmiennej niezależnej, y - jest zbiorem N odpowiadających wartości zmiennej zależnej, n - jest zadanym stopniem wielomianu aproksymującego, p - jest wektorem poszukiwanych wartości współczynników wielomianu aproksymującego.

4. x1=0, y1=0 x2=5, y2=1 x3=10, y3=0 % Obliczenie współczynników paraboli y(x)=p(1)*x^2 p(2)*x+p(3) x=[0,5,10]; y=[0,1,0]; • plot(x,y,'o') %wykres punktów p=polyfit(x,y,2) % sprawdzenie paraboli na wykresie x1=0:0.1:10; y1=p(1)*x1.^2+p(2)*x1+p(3); hold on plot(x1,y1) p = -0.0400 0.4000 0.0000

5. x1=0, y1=0 x2=2, y2=1 x3=3, y3=2 x4=10, y4=0 Nierównomierność kroku dla x – 4 punkty % Obliczenie współczynników paraboli y(x)=p(1)*x^3+p(2)*x1^2+p(3)*x1+p(4) clc; x=[0, 2, 3, 10]; y=[0, 1, 2, 0]; plot(x,y,'o') %wykres punktów – "kółeczka" p=polyfit(x,y,3) %teraz wyznaczamy wielomian 3-go stopnia % sprawdzenie paraboli na wykresie x1=0:0.1:10; y1=p(1)*x1.^3+p(2)*x1.^2+p(3)*x1+p(4); %taki zapis niewygodny hold on %zamrożenie wykresu plot(x1,y1) % wykres funkcji f(x) p = -0.0327 0.3304 -0.0298 0.0000

6. p1(0 0), p2(1 1) , p3(2 1) , p4(3 1) , p5(4 1) p6(5 0) • clc, clear • N=input('Podaj rząd wielomianu:') %interakcja: rząd wielomianu • x=[0 1 2 3 4 5] • y=[0 1 1 1 1 0] • p=polyfit(x,y,N)%współczynniki wielomianu • x1=0:0.1:5; • y1=0; • %Inny sposób: pętla sumująca elementy wielomianu • for m=1:N • y1=y1+p(m)*x1.^(N-m+1); • end; • %wykres punktów i funkcji • plot(x,y,'o',x1,y1) • axis([0 5 -1 2]) N=2 N=5

7. Wnioski: Wielomian aproksymujący powinien mieć rząd mniejszy niż liczba punktów – inaczej ostrzeżenie o niejednoznacznym rozwiązaniu reduce the degree of the polynomial

8. Wygodna metoda wstawiania wielu punktów do wielomianu o znanych współczynnikach -polyval • Wartości otrzymanego wielomianu, w dowolnych zadanych punktach wektorem xx oblicza funkcja polyval ypi = polyval(p, x) gdzie: x - tablica punktów zmiennej niezależnej p - wektor współczynników wielomianu ypi - wektor wartości wielomianu w zadanych punktach x

9. clc N=15 x=1:N y=round(10*rand(1,N)) %losowe wartości y plot(x, y, 'o') for k=4:7 % pętla stopni wielomianu od 4 do 7 hold on p=polyfit(x,y,k) xx=0:0.1:N; yy=polyval(p,xx); plot(xx,yy,'Color',[rand randrand],'LineWidth',k-3) end Przykładowe rozwiązanie

10. Aproksymacja - błędy clc,clear x = [0: 0.1: 2.5]; y=x.^6; p1 = polyfit(x,y,3) p2 = polyfit(x,y,6) f1 = polyval(p1,x); f2 = polyval(p2,x); subplot(2,1,1) plot(x,y,'o',x,f1,x, f2) subplot(2,1,2) plot(x,y-f1,'-',x, y-f2,'r-') błąd n=3 n=6 p1 = 53.0000 -126.3856 80.8523 -9.3899 p2 = 1.0000 -0.0000 0.0000 -0.00000.0000 -0.00000.0000

11. Jeszcze o wykresach… clc,clear for k=1:100 x(k)=k; y(k)=5*rand; end %wykres losowych punktów plot(x,y,'.','Color','k','MarkerSize',5) hold on %figura konturowa x=[10 10 90 90 10] y=[1 4 4 1 1] fill(x,y,'r','Linewidth',5) W funkcji fill domknięcie konturu jest automatyczne x=[10 10 90 90] y=[1 4 4 1]

12. Interpolacja Interpolacjato przybliżanie przy pomocy funkcji z pewnej klasy, np. wielomianami, funkcjami trygonometrycznymi, funkcjami sklejanymi itp. Znane są dokładne wartości funkcji w punktach węzłowych. Zadaniem interpolacji jest wyznaczenie funkcji interpolującej, co daje możliwość znalezienia przybliżonych wartości funkcji w punktach nie będących węzłami oraz oszacowanie błędów tych przybliżonych wartości.

13. W tym celu należy znaleźć funkcję W(x),nazywaną funkcją interpolującą, która w węzłach interpolacji przyjmuje takie same wartości co funkcja y=f(x). Interpolacja jest w pewnym sensie zadaniem odwrotnym do tablicowania funkcji. Przy tablicowaniu mając analityczną postać funkcji budujemy tablicę wartości, przy interpolacji natomiast na podstawie tablicy wartości funkcji określamy jej postać analityczną.

14. Interpolacją funkcji nazywa się wyznaczenie przybliżonych wartości funkcji f(x) dla dowolnego argumentu przy znanych jej wartościach w ustalonych punktach zwanych węzłami interpolacji.

15. Metody interpolacji 'nearest' - interpolacja "najbliższy sąsiad" 'linear' - liniowa 'spline' - przy pomocy funkcji sklejanych 'pchip' lub 'cubic' - sześcienna

16. Przykłady metod interpolacji clc clear x = 0:10; %węzeł y = sin(x); %węzeł xi = 0:0.25:10; %inne x figure(1) yi = interp1(x,y,xi,'linear'); plot(x,y,'o',xi,yi) figure(2) yi = interp1(x,y,xi,'cubic'); plot(x,y,'o',xi,yi) figure(3) yi = interp1(x,y,xi,'spline'); plot(x,y,'o',xi,yi)

17. Teraz dowolne węzły clc clear x = 0:10; y = [1 2.5 3 2.5 1 0 -1 -2.5 -3 -2.5 -1]; xi = 0:0.25:10; subplot(3,1,1) yi = interp1(x,y,xi,'linear'); plot(x,y,'o',xi,yi) subplot(3,1,2) yi = interp1(x,y,xi,'cubic'); plot(x,y,'o',xi,yi) subplot(3,1,3) yi = interp1(x,y,xi,'spline'); plot(x,y,'o',xi,yi)

18. EKSTRAPOLACJA – rozszerzenie poza przedział węzłów clear x = 0:10; y = sin(x); xi = -5:0.25:15; %szerszy przedział x subplot(3,1,1) yi = interp1(x,y,xi,'linear','extrap'); plot(x,y,'o',xi,yi) subplot(3,1,2) yi = interp1(x,y,xi,'spline','extrap'); plot(x,y,'o',xi,yi) subplot(3,1,3) yi = interp1(x,y,xi,'pchip','extrap'); plot(x,y,'o',xi,yi)

19. Interpolacja powierzchni [X,Y]=meshgrid(-3:0.4:3) Z=5*sin(X).*cos(Y) [XI,YI]=meshgrid(-3:0.1:3) ZI=interp2(X,Y,Z,XI,YI,'linear') mesh(X,Y,Z) hold mesh(XI,YI,ZI+25) holdoff axis([-3 3 -3 3 -5 30])

20. Iteracje ciągów factorial(n) = n! %Obliczenie liczby Eulera clc,format compact k=0:1000; e=sum(1./factorial(k)) %to samo z pętlą for e=0; for k=0:1000 e=e+1/factorial(k); end disp(e) %inny sposób k=0:1000; a=k+1; b=factorial(k); e=0.5*sum(a./b) e = 2.7183 e = 2.7183 e = 2.7183 suma innego ciągu

21. a nawet gdy nie znamy funkcji factorial, można sobie poradzić z liczeniem silni: N=6; S=1; for k=2:N S=S*k; end disp(S) 720

22. Iteracja (pętla) warunkowawhile (dopóki) whilewarunek_logiczny instrukcje end Badanie warunku (true/false) Jeśli warunek==true to instrukcje są wykonywane, jeśli warunek==false to koniec działania pętli warunkowej Ponowne sprawdzenie warunku … itd.

23. Algorytm iteracyjny z połowieniem przedziału dla znalezienia intervalbisectionmethod 12=1 za małe 22=4 za duże 1.52=2.25 za duże 1.252=1.5625 za małe 1.3752=1.8906 za małe itd. a b xi xi+1 xi+1 jeśli za duże jeśli za małe

24. Realizacja format long M = 2 a = 1 b = 2 k = 0; d=1e-7 % dokładność whileb-a > d x = (a + b)/2; if x^2 > M b = x; else a = x; end k = k + 1; end disp(a) disp(b) 2.2204e-016 a = 1.4142 b = 1.4142

25. Badanie zmiany znaku funkcji– też połowienie przedziału warunek: sign(f(a))~=sign(f(b)) funkcja (signum) - znak sign(x)= (-1 || +1) clc,clear k = 0; a = 3; b = 4; ifsign(sin(a))~=sign(sin(b)) whileabs(b-a) > 1e-10 x = (a + b)/2; %środek ifsign(sin(x)) == sign(sin(b)) b = x; % koniec=środek else a = x; % początek=środek end k = k + 1; end a b else disp('Ustal inne a i b') end a = 3.1416 b = 3.1416 ~= nierówne Warunek sign(f(a))~=sign(f(b)) można też zapisać: f(a)*f(b)<0

26. Metoda Newtona Zadaniem jest znalezienie pierwiastka równania zadanej funkcji ciągłej f(x) W metodzie Newtona przyjmuje się następujące założenia: W przedziale [a,b] znajduje się dokładnie jeden pierwiastek Funkcja ma różne znaki na krańcach przedziału, tj. f(a)*f(b)<0 Pierwsza i druga pochodna funkcji mają stały znak w tym przedziale. minus bo f(x)<0

27. u α f(x1) -f(x1)/u = tg(α) = f '(x1) u = x2-x1= -f(x1)/f '(x1) x2 = x1-f(x1)/f '(x1)

28. Wyznaczenie miejsca zerowego funkcji ln(x) f = log(x) df = 1/x x = 0.330258509299405 x = 0.696145164463784 x = 0.948286903858718 x = 0.998639213816122 x = 0.999999073710225 x = 0.999999999999571 x = 1 x = 1 >> clc, clear, format compact, format long syms x f=log(x) df=diff(f) a=0.1; b=4; x=a; xp=b; whileabs(x - xp) > eps xp = x; x = x - subs(f)/subs(df) end

29. Paproć Barnsley'a 85% 7% 7% 1%

30. Realizacja w Matlabie krok=500000; x=zeros(1, krok); y=zeros(1, krok); for n=1:krok r=rand(); %liczba losowa ifr<=0.01 x(n+1)=0; y(n+1)=0.16*y(n); elseifr<=0.08 x(n+1)=0.2*x(n)-0.26*y(n); y(n+1)=0.23*x(n)+0.22*y(n)+1.6; elseif r<=0.15 x(n+1)=-0.15*x(n)+0.28*y(n); y(n+1)=0.26*x(n)+0.24*y(n)+0.44; else x(n+1)=0.85*x(n)+0.04*y(n); y(n+1)=-0.04*x(n)+ 0.85*y(n)+1.6; end end plot(x,y,'.','Color','g','MarkerSize',1)

32. Metoda Monte Carlo "Strzelamy" losowo N razy w kwadratową tarczę o boku A. Liczymy liczbę trafień k w koło wpisane w kwadrat. Pole kwadratu: A*A Pole koła:  A2/4 stąd =4*k/N

33. Realizacja w Matlabie %Monte Carlo pi clc,clear N=1000000; x=zeros(1,N); y=zeros(1, N); k=0; for n=1:N r1=rand()-0.5; r2=rand()-0.5; if r1^2+r2^2<=0.25 x(n+1)=r1; y(n+1)=r2; k=k+1; end end plot(x,y,'.','Color','g','MarkerSize',1) naszePI=4*k/N

34. naszePI = 3.1411

35. Transformacje geometryczne x1 x2 x3 ……. y1 y2 y3 ……. 1 11 ….. Kontur zamknięty opisany współrzędnymi punktów w postaci: Macierz przesunięcia Prz 1 0 px 0 1 py 0 0 1 Macierz skalowania Ska sx 0 0 0 1 py 0 0 1 Macierz obrotu Obr cossin 0 -sincos 0 0 0 1 Przekształcenie= Psz*Ska*Obr

36. Realizacja w Matlabie clc,clear x=[0 0 1 1 3 3 0; 0 5 5 1 1 0 0; 1 111111] x(1,:) fill(x(1,:),x(2,:),'c') %axis([-1 6 -1 6]) axisequal symspxpysxsy kat Prz=[1 0 px;0 1 py;0 0 1] Ska=[sx 0 0;0 sy 0; 0 0 1] Obr=[cos(kat) sin(kat) 0; -sin(kat) cos(kat) 0; 0 0 1] Cal=Prz*Ska*Obr; x=Cal*x px=2; py=2; sx=0.5; sy=0.5; kat=pi/4; x=subs(x) hold on fill(x(1,:),x(2,:),'r')