MOiPP
Wykład 6
Matlab
" Aproksymacja
" Interpolacja
" Inne metody obliczeniowe
1
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
Wykorzystanie np. gdy mamy wiele dyskretnych
punktów pomiarowych i chcemy znalez
ć 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.
2
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:
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.
3
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
% sprawdzenie paraboli na wykresie
x1=0:0.1:10;
y1=p(1)*x1.^2+p(2)*x1+p(3);
hold on
1.4
plot(x1,y1) 1.2
1
0.8
0.6
0.4
0.2
p =
0 4
0 1 2 3 4 5 6 7 8 9 10
-0.0400 0.4000 0.0000
x1=0, y1=0
Nierównomierność kroku dla x  4 punkty
x2=2, y2=1
x3=3, y3=2
x4=10, y4=0
% 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
hold on %zamrożenie wykresu
plot(x1,y1) % wykres funkcji f(x)
5
4
3
p =
2
-0.0327 0.3304 -0.0298 0.0000
1
0
5
-1
0 1 2 3 4 5 6 7 8 9 10
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
%wykres punktów i funkcji
plot(x, y, 'o', x1, y1)
axis([0 5 -1 2])
2
2
N=2 N=5
1.5 1.5
1
1
0.5
0.5
0
0
-0.5
-0.5
6
-1
-1
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
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
7
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:
gdzie:
x - tablica punktów zmiennej niezależnej
p - wektor współczynników wielomianu
ypi - wektor wartości wielomianu w zadanych punktach x
8
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 rand rand],'LineWidth',k-3)
plot(xx,yy,'Color',[rand rand rand],'LineWidth',k-3)
end
25
20
15 Przykładowe
rozwiązanie
10
5
0
9
0 5 10 15
Aproksymacja - błędy
300
clc,clear
x = [0: 0.1: 2.5];
200
y=x.^6;
100
p1 = polyfit(x,y,3)
p2 = polyfit(x,y,6)
0
f1 = polyval(p1,x);
-100
0 0.5 1 1.5 2 2.5
f2 = polyval(p2,x);
subplot(2,1,1)
błąd
15
plot(x,y,'o',x,f1,x, f2)
10
10
subplot(2,1,2)
subplot(2,1,2)
n=3
n=3
5
plot(x,y-f1,'-',x, y-f2,'r-')
n=6
0
-5
-10
0 0.5 1 1.5 2 2.5
p1 =
53.0000 -126.3856 80.8523 -9.3899
p2 =
1.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000
10
Jeszcze o wykresach&
clc,clear
5
4.5
for k=1:100
4
x(k)=k;
3.5
3
y(k)=5*rand;
2.5
end
2
%wykres losowych punktów
1.5
1
plot(x,y,'.','Color','k','MarkerSize',5)
0.5
hold on
0
0 20 40 60 80 100
%figura konturowa
%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]
11
y=[1 4 4 1]
Interpolacja
Interpolacja to 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.
12
W tym celu należy znalez
ć 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.
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ą.
13
Interpolacją funkcji
nazywa się wyznaczenie przybliżonych wartości funkcji
f(x) dla dowolnego argumentu
przy znanych jej wartościach
przy znanych jej wartościach
w ustalonych punktach
zwanych węzłami interpolacji.
14
Metody interpolacji
'nearest' - interpolacja "najbliższy sąsiad"
'linear' - liniowa
'spline' - przy pomocy funkcji sklejanych
'pchip' lub 'cubic' - sześcienna
15
1
0.8
Przykłady metod interpolacji
0.6
0.4
0.2
clc
0
-0.2
clear
-0.4
-0.6
x = 0:10; %węzeł
-0.8
-1
y = sin(x); %węzeł 0 1 2 3 4 5 6 7 8 9 10
1
xi = 0:0.25:10; %inne x
0.8
0.6
figure(1)
0.4
0.2
yi = interp1(x,y,xi,'linear');
0
-0.2
plot(x,y,'o',xi,yi)
plot(x,y,'o',xi,yi)
-0.4
-0.6
figure(2)
-0.8
-1
yi = interp1(x,y,xi,'cubic'); 0 1 2 3 4 5 6 7 8 9 10
1
plot(x,y,'o',xi,yi)
0.8
0.6
figure(3)
0.4
0.2
yi = interp1(x,y,xi,'spline');
0
-0.2
plot(x,y,'o',xi,yi)
-0.4
-0.6
-0.8
16
-1
0 1 2 3 4 5 6 7 8 9 10
Teraz dowolne węzły
4
clc
2
0
clear
-2
x = 0:10;
-4
0 2 4 6 8 10
y = [1 2.5 3 2.5 1 0 -1 -2.5 -3 -2.5 -1];
xi = 0:0.25:10;
4
2
subplot(3,1,1)
0
yi = interp1(x,y,xi,'linear');
-2
plot(x,y,'o',xi,yi)
plot(x,y,'o',xi,yi)
-4
0 2 4 6 8 10
subplot(3,1,2)
4
yi = interp1(x,y,xi,'cubic');
2
plot(x,y,'o',xi,yi)
0
subplot(3,1,3)
-2
yi = interp1(x,y,xi,'spline');
-4
0 2 4 6 8 10
plot(x,y,'o',xi,yi)
17
EKSTRAPOLACJA
 rozszerzenie poza przedział węzłów
2
clear 0
x = 0:10;
-2
y = sin(x);
-4
xi = -5:0.25:15; %szerszy przedział x
-6
-5 0 5 10 15
subplot(3,1,1)
15
yi = interp1(x,y,xi,'linear','extrap');
10
plot(x,y,'o',xi,yi)
plot(x,y,'o',xi,yi)
5
subplot(3,1,2)
0
yi = interp1(x,y,xi,'spline','extrap');
-5
plot(x,y,'o',xi,yi)
-10
-5 0 5 10 15
subplot(3,1,3)
40
yi = interp1(x,y,xi,'pchip','extrap');
30
plot(x,y,'o',xi,yi)
20
10
0
18
-10
-5 0 5 10 15
Interpolacja powierzchni
30
25
20
15
10
5
[X,Y]=meshgrid(-3:0.4:3)
0
-5
Z=5*sin(X).*cos(Y)
2
3
[XI,YI]=meshgrid(-3:0.1:3)
[XI,YI]=meshgrid(-3:0.1:3)
2
2
0 1
0
-1
ZI=interp2(X,Y,Z,XI,YI,'linear')
-2
-2
-3
mesh(X,Y,Z)
hold
mesh(XI,YI,ZI+25)
hold off
axis([-3 3 -3 3 -5 30])
19
factorial(n) = n!
Iteracje ciągów
%Obliczenie liczby Eulera
clc,format compact
k=0:1000;
e=sum(1./factorial(k))
%to samo z pętlą for
e=0;
e =
for k=0:1000
for k=0:1000
2.7183
e=e+1/factorial(k);
end
e =
disp(e)
2.7183
suma innego ciągu
%inny sposób
e =
k=0:1000;
2.7183
a=k+1;
b=factorial(k);
e=0.5*sum(a./b)
20
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
end
disp(S)
720
21
Iteracja (pętla) warunkowa while (dopóki)
while warunek_logiczny
instrukcje
end
1. Badanie warunku (true/false)
2. Jeśli warunek==true to instrukcje są wykonywane,
jeśli warunek==false to koniec działania pętli
warunkowej
3. Ponowne sprawdzenie warunku
4. & itd.
22
Algorytm iteracyjny z połowieniem
przedziału dla znalezienia
interval bisection method
2
1 =1 za małe
2
2 =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
1.3752=1.8906 za małe
itd.
a b
xi
xi+1 xi+1
jeśli za jeśli za
duże małe
23
Realizacja
format long
M = 2
a = 1
b = 2
k = 0;
d=1e-7 % dokładność
while b-a > d
x = (a + b)/2;
2.2204e-016
if x^2 > M
if x^2 > M
a =
a =
b = x;
else
1.4142
a = x;
b =
end
k = k + 1;
1.4142
end
disp(a)
disp(b)
24
Badanie zmiany znaku funkcji  też połowienie przedziału
warunek: sign(f(a))~=sign(f(b))
funkcja (signum) - znak
clc,clear
k = 0;
sign(x)= (-1 || +1)
a = 3;
b = 4;
if sign(sin(a))~=sign(sin(b))
while abs(b-a) > 1e-10
a =
x = (a + b)/2; %środek
~=
if sign(sin(x)) == sign(sin(b)) 3.1416
if sign(sin(x)) == sign(sin(b)) 3.1416
nierówne
nierówne
b = x; % koniec=środek
b =
else
a = x; % początek=środek
3.1416
end
k = k + 1;
Warunek
end
sign(f(a))~=sign(f(b))
a
można też zapisać:
b
f(a)*f(b)<0
else
disp('Ustal inne a i b')
end
25
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:
1. W przedziale [a,b] znajduje się dokładnie jeden pierwiastek
2. Funkcja ma różne znaki na krańcach przedziału, tj. f(a)*f(b)<0
3. Pierwsza i druga pochodna funkcji mają stały znak w tym przedziale.
minus bo f(x)<0
26
u
ą
f(x1)
-f(x1)/u = tg(ą) = f '(x1)
u = x2-x1= -f(x1)/f '(x1)
x2 = x1-f(x1)/f '(x1)
27
Wyznaczenie miejsca zerowego funkcji ln(x)
f =
log(x)
clc, clear, format compact, format long
df =
1/x
syms x
x =
f=log(x)
0.330258509299405
x =
df=diff(f)
0.696145164463784
a=0.1;
x =
b=4;
0.948286903858718
x=a;
x =
0.998639213816122
xp=b;
x =
while abs(x - xp) > eps
0.999999073710225
xp = x;
x =
x = x - subs(f)/subs(df)
0.999999999999571
x =
end
1
x =
1
>>
28
Paproć Barnsley'a
85%
7%
7%
1%
29
Realizacja w Matlabie
krok=500000;
x=zeros(1, krok);
y=zeros(1, krok);
for n=1:krok
r=rand(); %liczba losowa
if r<=0.01
x(n+1)=0;
y(n+1)=0.16*y(n);
elseif r<=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;
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)
30
31
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
32
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;
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
33
naszePI =
3.1411
34
Transformacje geometryczne
x1 x2 x3 & & .
Kontur zamknięty opisany
y1 y2 y3 & & .
współrzędnymi punktów w postaci:
1 1 1 & ..
Macierz przesunięcia Prz 1 0 px
0 1 py
0 0 1
Macierz skalowania Ska sx 0 0
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
35
Realizacja w Matlabie
clc,clear
x=[0 0 1 1 3 3 0; 0 5 5 1 1 0 0; 1 1 1 1 1 1 1]
x(1,:)
fill(x(1,:),x(2,:),'c')
%axis([-1 6 -1 6])
axis equal
syms px py sx sy kat
Prz=[1 0 px;0 1 py;0 0 1]
Ska=[sx 0 0;0 sy 0; 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')
36
37