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ę