Aproksymacja i
interpolacja
PWSZ, Elbląg 2002
Aproksymacja i
interpolacja
PWSZ, Elbląg 2002
Elbąg, PWSZ 2002r.
2
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.
Niech poszukiwana jest krzywa dla zadanej liczby
punktów:
jest opisana równaniem:
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.
)
(x
F
y
)
,
(
i
i
y
x
)
(
...
)
(
)
(
)
(
2
2
1
1
x
f
c
x
f
c
x
f
c
x
F
n
n
Elbąg, PWSZ 2002r.
3
Regresja liniowa
Najbardziej podstawową i najprostszą metodą aproksymacji
średniokwadratowej jest aproksymacja funkcją liniową czyli
regresja liniowa. Wówczas:
Pozostałe funkcje:
Dla kolejnych punktów otrzymujemy:
1
)
(
,
)
(
2
1
x
f
x
x
f
1
)
(
x
f
j
m
m
m
m
n
y
y
y
c
c
x
x
x
y
c
x
c
y
c
x
c
y
c
x
c
2
1
2
1
2
1
2
2
2
2
1
1
2
1
1
1
1
1
lub
Elbąg, PWSZ 2002r.
4
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:
gdzie: r – wektor pionowych odległości pomiędzy poszukiwaną
krzywą a zadanymi punktami. Szukane jest takie rozwiązanie, dla
którego:
lub
macierzowo
osiąga
minimum.
Stąd:
Ac
y
r
m
i
i
r
1
2
r
r
T
Ac
A
c
Ac
y
y
y
Ac
A
c
y
A
c
Ac
y
y
y
Ac
y
Ac
y
r
r
T
T
T
T
T
T
T
T
T
T
T
T
2
Elbąg, PWSZ 2002r.
5
Regresja liniowa
Iloczyn ten osiągnie minimum jeśli:
stąd:
Powyższe równanie nazywane jest równaniem aproksymacji.
0
0
Ac
A
y
A
r
r
dc
d
T
T
T
y
A
c
A
A
T
T
Elbąg, PWSZ 2002r.
6
Regresja liniowa
Przykład 1:
Wyznaczenie prostej aproksymującej punkty o współrzędnych
(1,1), (2,2), (4,2), (5,3):
x=[1 2 4 5]; y=[1 2 2 3];
x=x'; y=y';
A=[x ones(size(x))];
c=(A'*A)\(A'*y)
ya=c(1)*x+c(2)
plot(x,y,'o',x,ya,'-')
grid on
xlabel('x')
ylabel('y=F(x)')
1
1.5
2
2.5
3
3.5
4
4.5
5
1
1.2
1.4
1.6
1.8
2
2.2
2.4
2.6
2.8
3
x
y=
F
(x
)
Elbąg, PWSZ 2002r.
7
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.
8
Regresja liniowa
%% Przygotowanie danych do obliczen
clear
x=[0.1:0.1:10]';
[m,n]=size(x);
for i=1:m
y(i,1)=i*0.1*rand(1,1);
end
plot(x,y,'.');
%% poszukiwanie funkcji y=ax
%% najlepiej w sensie sumy kwadratów
%% przyblizajace te dane
a=x\y;
plot(x,y,x,a*x);
0
1
2
3
4
5
6
7
8
9
10
0
1
2
3
4
5
6
7
8
9
10
Elbąg, PWSZ 2002r.
9
Regresja liniowa
0
1
2
3
4
5
6
7
8
9
10
0
1
2
3
4
5
6
7
8
9
10
Elbąg, PWSZ 2002r.
10
Aproksymacja
Równanie aproksymacji:
jest prawdziwe dla dowolnej funkcji aproksymacji:
gdzie: - nieznane współczynniki,
- funkcje bazowe,
zaś macierze A, c i y:
y
A
c
A
A
T
T
)
(
...
)
(
)
(
)
(
2
2
1
1
x
f
c
x
f
c
x
f
c
x
F
n
n
j
c
x
f
j
m
n
m
n
m
m
n
n
y
y
y
y
c
c
c
c
x
f
x
f
x
f
x
f
x
f
x
f
x
f
x
f
x
f
A
2
1
2
1
2
1
2
2
2
2
1
1
1
2
1
1
,
,
Elbąg, PWSZ 2002r.
11
Aproksymacja
Przykład 3:
Dla
zadanego
zbioru
punktów
poszukiwana
jest
funkcja
aproksymująca:
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;
x
c
x
c
y
2
1
Elbąg, PWSZ 2002r.
12
Aproksymacja
plot(x,y,'o',xa,ya,'-');
xlabel('x'); ylabel('y=F(x)');
legend('dane','aproksymacja')
0.5
1
1.5
2
2.5
3
3.5
5
5.5
6
x
y=
F
(x
)
dane
aproksymacja
Elbąg, PWSZ 2002r.
13
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:
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 x a y.
0
1
1
1
...
a
x
a
x
a
x
a
x
W
r
r
r
r
Elbąg, PWSZ 2002r.
14
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.
15
Aproksymacja wielomianem
0
1
2
3
4
5
6
7
8
9
10
-1
-0.5
0
0.5
1
1.5
2
2.5
Elbąg, PWSZ 2002r.
16
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:
2) interpolacja wielomianem Lagrange’a:
3) interpolacja wielomianem Newtona:
n
n
n
n
n
c
x
c
x
c
x
c
x
P
1
2
2
1
1
1
n
j
k
k
k
j
k
j
n
n
n
x
x
x
x
x
L
gdzie
x
L
y
x
L
y
x
L
y
x
P
,
1
2
2
1
1
1
,
n
n
n
x
x
x
x
x
x
c
x
c
x
x
c
x
x
c
c
x
P
2
1
2
1
3
1
2
1
1
Elbąg, PWSZ 2002r.
17
Interpolacja
Elbąg, PWSZ 2002r.
18
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:
Punkty łączenia wielomianów nazywamy węzłami.
We węzłach sprawdzane są warunki ciągłości (np. ciągłość
pochodnych).
x
P
i
1
i
i
x
x
x
Elbąg, PWSZ 2002r.
19
Interpolacja
Interpolacja kawałkami liniowa
Przykład 5:
x1=linspace(0,2*pi,100);
x2=linspace(0,2*pi,6);
plot(x1,sin(x1),x2,sin(x2))
grid on
xlabel('x')
ylabel('plot(x,sin(x)')
legend('x=linspace(0,2*pi,100)',...
'x=linspace(0,2*pi,6)')
0
1
2
3
4
5
6
7
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
x
pl
ot
(x
,s
in
(x
)
x=linspace(0,2*pi,100)
x=linspace(0,2*pi,6)
Elbąg, PWSZ 2002r.
20
Interpolacja
0
1
2
3
4
5
6
7
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
x
pl
ot
(x
,s
in
(x
)
x=linspace(0,2*pi,100)
x=linspace(0,2*pi,6)
Elbąg, PWSZ 2002r.
21
Interpolacja
Interpolacja kawałkami sześcienna – Hermite’a
gdzie:
są poszukiwanymi współczynnikami, dla
których są spełnione we węzłach następujące warunki:
1) ciągłości:
2) znane są wartości pierwszej pochodnej i jej ciągłość:
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:
3
2
i
i
i
i
i
i
i
i
x
x
d
x
x
c
x
x
b
a
x
P
i
i
i
i
d
c
b
a
,
,
,
i
i
i
i
x
P
x
P
1
i
i
i
i
x
P
x
P
1
Elbąg, PWSZ 2002r.
22
Interpolacja
1) ciągłość drugiej pochodnej:
2) pierwsze pochodne (nachylenia krzywej) muszą być znane na
końcach
przedziału
- ustalone nachylenie:
- naturalne nachylenie:
- nachylenie nieznane:
i
i
i
i
x
P
x
P
1
n
n
x
P
x
P
,
1
1
,
2
,
1
1
1
st
x
P
st
x
P
n
n
,
0
1
1
n
n
x
P
x
P
2
2
2
1
2
2
2
1
x
P
x
P
i
x
P
x
P
Elbąg, PWSZ 2002r.
23
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 interp1 umoż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.
24
Interpolacja
Elementy wektora x muszą tworzyć ciąg rosnący, dodatkowo w
przypadku interpolacji wielomianami trzeciego stopnia przyrosty
wartości elementów wektora x muszą 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.
25
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.
26
Interpolacja
0
2
4
6
8
10
12
14
16
18
20
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
x
y
Interpolacja kawalkami liniowa
Elbąg, PWSZ 2002r.
27
Interpolacja
0
2
4
6
8
10
12
14
16
18
20
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
x
y
Interpolacja kawalkami szescienna
Elbąg, PWSZ 2002r.
28
Interpolacja
0
2
4
6
8
10
12
14
16
18
20
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
x
y
Interpolacja funkcjami sklejanymi
Dziękuję za uwagę