Aproksymacja i interpolacja
PWSZ, ElblÄ…g 2002
Aproksymacja
Aproksymacja jest to przybli\
anie, zastępowanie jednych wielkości
drugimi. Aproksymacja mo\
e dotyczyć dowolnych wielkości
matematycznych  liczb, funkcji, krzywych, obszarów, wektorów, macierzy.
Zastępowanie danej wielkości inną obarczone jest pewnym błędem.
Oszacowanie wielkości błędu pozwala na ocenę, czy dane przybli\
enie
jest zadawalajÄ…ce, czy te\
nie.
y = F(x)
Niech poszukiwana jest krzywa dla zadanej liczby
punktów:
(xi, yi )
jest opisana równaniem:
F(x) = c1 f1(x) + c2 f2(x) +...+ cn fn (x)
Aproksymacja stosowana jest wówczas, gdy ilość zadanych punktów m
jest mniejsza od ilości nieznanych współczynników n krzywej F(x).
Zwykle nie mo\
na przeprowadzić krzywej przez wszystkie punkty.
Poszukiwana jest wówczas najbli\
sza krzywa w sensie minimum
kwadratu błędu.
ElbÄ…g, PWSZ 2002r. 2
Regresja liniowa
Najbardziej podstawowÄ… i najprostszÄ… metodÄ… aproksymacji
średniokwadratowej jest aproksymacja funkcją liniową czyli regresja
liniowa. Wówczas:
f1(x) = x, f2(x) =1
Pozostałe funkcje:
f (x) = 1
j
Dla kolejnych punktów otrzymujemy:
Å„Å‚
x1 1 y1
îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚
ôÅ‚
c1x1 + c2 = y1
ïÅ‚ ïÅ‚ śł
x2 1śł y2
ôÅ‚
ïÅ‚ śłîÅ‚c1 Å‚Å‚ = ïÅ‚ śł
c1x2 + c2 = y2 lub
òÅ‚
ïÅ‚ śłïÅ‚ śł ïÅ‚ śł
Å› Å›
ðÅ‚c2 ûÅ‚
ôÅ‚
Å›
ïÅ‚ śł ïÅ‚ śł
ôÅ‚c xm + c2 = ym
m m
ðÅ‚x 1ûÅ‚ ðÅ‚y ûÅ‚
ół n
ElbÄ…g, PWSZ 2002r. 3
Regresja liniowa
co mo\
na zapisać w postaci macierzowej Ac=y. Równanie to dla m>n nie
ma dokładnego rozwiązania, stąd:
r = y - Ac
gdzie: r  wektor pionowych odległości pomiędzy poszukiwaną krzywą a
zadanymi punktami. Szukane jest takie rozwiązanie, dla którego:
m
2
rT r
lub macierzowo osiÄ…ga minimum.
"r
i
i=1
StÄ…d:
T
rT r = (y - Ac) (y - Ac)
= yT y - yT Ac - cT AT y + cT AT Ac
= yT y - 2yT Ac + cT AT Ac
ElbÄ…g, PWSZ 2002r. 4
Regresja liniowa
Iloczyn ten osiągnie minimum jeśli:
d
(rT r)= 0 - AT y + AT Ac = 0
dc
stÄ…d:
(AT A)c = AT y
Powy\
sze równanie nazywane jest równaniem aproksymacji.
ElbÄ…g, PWSZ 2002r. 5
Regresja liniowa
Przykład 1:
Wyznaczenie prostej aproksymującej punkty o współrzędnych (1,1), (2,2),
(4,2), (5,3):
3
x=[1 2 4 5]; y=[1 2 2 3];
2.8
x=x'; y=y';
2.6
A=[x ones(size(x))]; 2.4
2.2
c=(A'*A)\(A'*y)
2
ya=c(1)*x+c(2)
1.8
plot(x,y,'o',x,ya,'-')
1.6
grid on
1.4
xlabel('x')
1.2
ylabel('y=F(x)')
1
1 1.5 2 2.5 3 3.5 4 4.5 5
x
ElbÄ…g, PWSZ 2002r. 6
y=F(x)
Regresja liniowa
W programie MATLAB dzielenie lewostronne jest równoznaczne operacji:
y=A\y=(A *A)\(A *y)
Je\
eli macierz A ma wymiar m x n, gdzie m>n program MATLAB
rozwiązuje zagadnienie regresji liniowej (znajduje współczynniki linii
prostej).
Przykład 2:
Przedstawione poni\
ej polecenia realizujÄ… aproksymacjÄ™ funkcjÄ… liniowÄ….
Wygenerowane dane umieszczono w wektorach kolumnowych x i y. Przy
generowaniu danych u\
yto funkcjÄ™ rand generujÄ…ca liczby pseudolosowe
o rozkładzie normalnym o zadanej wartości średniej i wariancji. Dane
generowane są w pętli, w ka\
dym jej kroku zwiększa się średnią
generowanych liczb.
ElbÄ…g, PWSZ 2002r. 7
Regresja liniowa
%% Przygotowanie danych do obliczen
10
clear
9
x=[0.1:0.1:10]'; 8
7
[m,n]=size(x);
6
for i=1:m
5
4
y(i,1)=i*0.1*rand(1,1);
3
end
2
1
plot(x,y,'.');
0
0 1 2 3 4 5 6 7 8 9 10
%% poszukiwanie funkcji y=ax
%% najlepiej w sensie sumy kwadratów
%% przyblizajace te dane
a=x\y;
plot(x,y,x,a*x);
ElbÄ…g, PWSZ 2002r. 8
Regresja liniowa
10
9
8
7
6
5
4
3
2
1
0
0 1 2 3 4 5 6 7 8 9 10
ElbÄ…g, PWSZ 2002r. 9
Aproksymacja
Równanie aproksymacji:
(AT A)c = AT y
jest prawdziwe dla dowolnej funkcji aproksymacji:
F(x) = c1 f1(x) + c2 f2(x) +...+ cn fn (x)
f (x)
cj
gdzie: - nieznane współczynniki, - funkcje bazowe, zaś
j
macierze A, c i y:
f1(x1) f2(x1) š fn(x1) c1 y1
îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚
ïÅ‚ ïÅ‚ śł
f1(x2) f2(x2) š fn(x2)śł, c = ïÅ‚c śł y2
2
ïÅ‚ śł ïÅ‚ śł, y = ïÅ‚ śł
A =
ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł
Å› Å› Å›
ïÅ‚ ïÅ‚ śł ïÅ‚ śł
f1(xm) f2(xm) š fn(xm)śł
n m
ðÅ‚ ûÅ‚ ðÅ‚c ûÅ‚ ðÅ‚y ûÅ‚
ElbÄ…g, PWSZ 2002r. 10
Aproksymacja
Przykład 3:
Dla zadanego zbioru punktów poszukiwana jest funkcja aproksymująca:
c1
y = + c2x
x
x=[0.955 1.380 1.854 2.093 2.674 3.006 3.255]';
y=[5.722 4.812 4.727 4.850 5.011 5.253 5.253]';
A=[1./x x]; %przygotowanie macierzy A
c=(A'*A)\(A'*y)
%funkcja linspace generuje wektor
%wierszowy o elementach równoodleglych
%w zadanym zakresie
xa=linspace(min(x),max(x),100);
xa=xa';
Aa=[1./xa xa];
ya=Aa*c;
ElbÄ…g, PWSZ 2002r. 11
Aproksymacja
plot(x,y,'o',xa,ya,'-');
xlabel('x'); ylabel('y=F(x)');
legend('dane','aproksymacja')
6
dane
aproksymacja
5.5
5
0.5 1 1.5 2 2.5 3 3.5
x
ElbÄ…g, PWSZ 2002r. 12
y=F(x)
Aproksymacja wielomianem
Aproksymacja funkcjÄ… liniowÄ… mo\
e okazać się nie wystarczająca
wówczas, gdy między danymi występuje bardziej zło\
ona zale\
ność.
Stosuje się wówczas zazwyczaj aproksymację wielomianem:
W(x)= ar xr + ar -1xr -1 + ... + a1x + a0
którą w programie MATLAB mo\
na zrealizować przy pomocy funkcji
polyfit:
a=polyfit(x,y,r)
dla danych wektorów x i y wyznaczającej współczynniki wielomianu
stopnia r przybli\
ającego najlepiej w sensie średniokwadratowym
zale\
ność między serią danych xa y.
ElbÄ…g, PWSZ 2002r. 13
Aproksymacja wielomianem
Przykład 4:
Przedstawione poni\
ej polecenia generują dane losowe, a następnie
przybli\
ajÄ… zale\
ność pomiędzy nimi wielomianami stopni od 1 do 10:
clear
close all
x=[0.1:0.1:10]';
[m,n]=size(x);
for i=1:m
y(i,1)=(sin(0.05*i)+2*cos(0.08*i))*rand(1,1);
end
d=1
for N=1:d:10
figure(N)
a=polyfit(x,y,N);
ya(:,N/d)=polyval(a,x);
plot(x,y,'.',x,ya);
end
ElbÄ…g, PWSZ 2002r. 14
Aproksymacja wielomianem
2.5
2
1.5
1
0.5
0
-0.5
-1
0 1 2 3 4 5 6 7 8 9 10
ElbÄ…g, PWSZ 2002r. 15
Interpolacja
Interpolacja polega na poszukiwaniu funkcji pomiędzy znanymi punktami
(podobnie jak aproksymacja). W odró\
nieniu jednak od aproksymacji
funkcja ta przechodzi przez te punkty. Je\
eli poszukiwana jest funkcja
poza zakresem zadanych punktów mamy do czynienia z ekstrapolacją.
Wybrane metody interpolacji:
1) interpolacja wielomianem:
Pn-1(x)= c1xn-1 + c2xn-2 +š + cn-1x + cn
2) interpolacja wielomianem Lagrange a:
n
x - xk
Pn-1(x)= y1L1(x)+ y2L2(x)+š + ynLn(x), gdzie Lj(x)=
"
xj - xk
k =1,k `" j
3) interpolacja wielomianem Newtona:
Pn-1(x)= c1 + c2(x - x1)+ c3(x - x1)(c - x2)+š + cn(x - x1)(x - x2)› (x - xn)
ElbÄ…g, PWSZ 2002r. 16
Interpolacja
ElbÄ…g, PWSZ 2002r. 17
Interpolacja
Interpolacja wielomianami sklejanymi
W interpolacji wielomianami sklejanymi zamiast stosowania jednego
wielomianu dla wszystkich danych punktów stosowane jest wiele
wielomianów niskiego poziomu dla danego przedziału danych:
Pi(x)
xi d" x d" xi+1 Punkty łączenia wielomianów nazywamy węzłami. We
węzłach sprawdzane są warunki ciągłości (np. ciągłość pochodnych).
ElbÄ…g, PWSZ 2002r. 18
Interpolacja
Interpolacja kawałkami liniowa
1
x=linspace(0,2*pi,100)
x=linspace(0,2*pi,6)
0.8
0.6
0.4
0.2
Przykład 5:
0
x1=linspace(0,2*pi,100);
-0.2
x2=linspace(0,2*pi,6);
-0.4
-0.6
plot(x1,sin(x1),x2,sin(x2))
-0.8
grid on
-1
0 1 2 3 4 5 6 7
xlabel('x')
x
ylabel('plot(x,sin(x)')
legend('x=linspace(0,2*pi,100)',...
'x=linspace(0,2*pi,6)')
ElbÄ…g, PWSZ 2002r. 19
plot(x,sin(x)
Interpolacja
1
x=linspace(0,2*pi,100)
x=linspace(0,2*pi,6)
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
0 1 2 3 4 5 6 7
x
ElbÄ…g, PWSZ 2002r. 20
plot(x,sin(x)
Interpolacja
Interpolacja kawałkami sześcienna  Hermite a
2 3
Pi(x)= ai + bi(x - xi )+ ci(x - xi ) + di(x - xi )
gdzie: są poszukiwanymi współczynnikami, dla których są
ai, bi, ci, di
spełnione we węzłach następujące warunki:
1) ciągłości:
Pi-1(xi)= Pi(xi)
2) znane są wartości pierwszej pochodnej i jej ciągłość:
2 2
Pi-1(xi)= Pi (xi)
Interpolacja kawałkami sześcienna  splajny
Metoda ta w odró\
nieniu od interpolacji Hermite a nie wymaga znajomości
pochodnych we wszystkich punktach węzłowych, muszą być jednak
spełnione następujące warunki:
ElbÄ…g, PWSZ 2002r. 21
Interpolacja
1) ciągłość drugiej pochodnej:
2 2 2 2
Pi-1(xi)= Pi (xi)
2) pierwsze pochodne (nachylenia krzywej) muszą być znane na końcach
przedziału
2 2
P1 (x1), Pn (xn)
- ustalone nachylenie: 2 2
P1 (x1)= st1, Pn (xn)= st2,
2 2 2 2
P1 (x1)= Pn (xn)= 0,
- naturalne nachylenie:
2 2 2 2 2 2 2 2 2 2 2 2
P1 (x2)= P2 (x2)i P1 (x2)= P2 (x2)
- nachylenie nieznane:
ElbÄ…g, PWSZ 2002r. 22
Interpolacja
Program MATLAB realizuje interpolację za pomocą następujących metod:
1) interpolacja kawałkami liniowa i sześcienna,
2) interpolacja za pomocÄ… funkcji sklejanych.
Funkcja interp1:
yi=interp1(x, y, x1,  metoda )
Funkcja interp1umo\
liwia wykonanie interpolacji funkcji jednej zmiennej
w punktach określonych wektorem xi. Węzły interpolacji określone są
parametrami x i y. Parametr  metoda umo\
liwia wybór metody
interpolacji:
1)  linear  interpolacja funkcją łamaną (kawałkami liniowa),
2)  spline  interpolacja funkcjami sklejanymi trzeciego stopnia,
3)  cubic  interpolacja wielomianami trzeciego stopnia (kawałkami
sześcienna),
ElbÄ…g, PWSZ 2002r. 23
Interpolacja
Elementy wektora xmuszą tworzyć ciąg rosnący, dodatkowo w przypadku
interpolacji wielomianami trzeciego stopnia przyrosty wartości elementów
wektora xmuszą być sobie równe.
Przykład 6:
Interpolacja ró\
nymi metodami:
x=0:20; y=sin(x)+sin(2*x);
xi=0:0.2:20;
yi=interp1(x,y,xi,'linear');
plot(x,y,'o',xi,yi,xi,sin(xi)+sin(2*xi));
xlabel('x');
ylabel('y');
title('Interpolacja kawalkami liniowa')
ElbÄ…g, PWSZ 2002r. 24
Interpolacja
figure(2)
yi=interp1(x,y,xi,'cubic');
plot(x,y,'o',xi,yi,xi,sin(xi)+sin(2*xi));
xlabel('x');
ylabel('y');
title('Interpolacja kawalkami szescienna')
figure(3)
yi=interp1(x,y,xi,'spline');
plot(x,y,'o',xi,yi,xi,sin(xi)+sin(2*xi));
xlabel('x');
ylabel('y');
title('Interpolacja funkcjami sklejanymi')
ElbÄ…g, PWSZ 2002r. 25
Interpolacja
Interpolacja kawalkami liniowa
2
1.5
1
0.5
0
-0.5
-1
-1.5
-2
0 2 4 6 8 10 12 14 16 18 20
x
ElbÄ…g, PWSZ 2002r. 26
y
Interpolacja
Interpolacja kawalkami szescienna
2
1.5
1
0.5
0
-0.5
-1
-1.5
-2
0 2 4 6 8 10 12 14 16 18 20
x
ElbÄ…g, PWSZ 2002r. 27
y
Interpolacja
Interpolacja funkcjami sklejanymi
2
1.5
1
0.5
0
-0.5
-1
-1.5
-2
0 2 4 6 8 10 12 14 16 18 20
x
ElbÄ…g, PWSZ 2002r. 28
y
Dziękuję za uwagę