28.11.02
INTERPOLACJA
Wstęp
Załóżmy, że mamy zbiór danych pomiarowych np. populację Krakowa w latach : 1970,1980,1990. Lata są węzłami, natomiast liczba mieszkańców to wartości w odpowiednich węzłach. Zagadnienie interpolacji polega na dopasowaniu funkcji do posiadanych danych. W naszym przykładzie pozwala nam oszacować populację Krakowa w latach np. 1975 czy nawet 1995. Jedną z najbardziej użytecznych i dobrze znanych klas funkcji używanych do interpolacji jest klasa wielomianów tzn. funkcji postaci : Pn( x) = a 0 + a 1 x + . . . + anxn, gdzie n jest nieujemną liczbą całkowitą i a 0 , . . . , an są rzeczywistymi współczynnikami. Jedną z głównych cech wielomianów jest to, że aproksymują one w sposób jednoznaczny dowolną, ciągłą funkcję. Traktuje o tym twierdzenie Weierstrass’a :
Jeżeli dana jest funkcja f ciągła na przedziale [ a, b] oraz > 0 to istnieje wielomian P taki, że
|f ( x) − P ( x) | <
dla każ dego x ∈ [ a, b]
Z pośród wielomianów, do problemów interpolacji najczęściej używa się wielomianów Lagrange, Hermite, Newtona oraz funkcji sklejanych.
Interpolacja Lagrange
Załóżmy, że mamy n węzłów {x 0 , x 1 , . . . , xn} oraz wartości {y 0 , y 1 , . . . , yn} w tych węzłach. Szukamy pewnej klasy funkcji F ( x) takiej, że F ( xi) = yi.
Niech Wn( xi) = P n
a
i=0
ixi = a 0 + a 1 x + . . . + anxn. Chcemy tak dobrać ai aby Wn( xi) = yi. Tym sposobem otrzymujemy następujący układ równań :
W 0( x 0) = a 0 + a 1 x 0 + a 2 x 2 + . . . + anxn
0
0
W 1( x 1) = a 0 + a 1 x 1 + a 2 x 2 + . . . + a 1
nxn
1
Aa = y < = >
..
.
Wn( xn) = a 0 + a 1 xn + a 2 x 2 + . . . + a n
nxn
n
, gdzie :
a
0
y 0
1
x 0
. . .
xn
0
a 1
y 1
1
x 1
. . .
xn
1
a =
.
,
y = . ,
A = .
.
.
.
.
.
.
.
.
.
.
.
.
an
yn
1
xn
. . .
xn
n
A jest wyznacznikiem macierzy Vandermode’a :
Y
detA =
( xi − xj)
i,j
i>j
Przy założeniu, że xi 6= xj układ ma dokładnie jedno rozwiązanie ( detA 6= 0). Poszukiwany wielomian Lagrange’a otrzymujemy zatem ze wzoru :
n
X
Ln( x)
=
yjϕj( x)
j=0
Ω n( x)
ϕj
=
( x − xj)Ω 0 ( x)
n
n
Y
Ω n( x)
=
( x − xi)
i=0
Szacowanie błędu :
M
|
n+1
f ( x) − Ln( x) | ¬
|Ω n( x) |
( n + 1)!
Mn+1 = sup|f ( n+1)( ξ) | ,gdzie ξ ∈ [ x 0 , xn]
Przykład 1
x
1
2
3
4
Wielomian Legendre’a dla :
ma postać : L
( − 4 + x)( − 1 + x)2
y
0
1
2
0
3( x) = − 1
2
2
1.5
1
0.5
1.5
2
2.5
3
3.5
4
Interpolacja Newtona
Metody polegające na bezpośrednim wyznaczaniu wielomianów interpolacyjnych z danych stablicowanych nazywane są metodami ilorazów różnicowych. Metody te były powszechnie używane zanim komputery cyfrowe stały się powszechnie dostępne.
Załóżmy, że mamy wielomian Lagrange’a stopnia co najwyżej n przechodzący przez stablicowane wartości funkcji f . Aby wyznaczyć iloraz różnicowy funkcji f przedstawmy wielomian Lagrange’a w postaci : Pn( x) = a 0 + a 1( x − x 0) + a 2( x − x 0)( x − x 1) + . . . + an( x − x 0)( x − x 1) . . . ( x − xn− 1) Współczynnik a 0 jest równy f ( x 0) ponieważ dalsze wyrazy Pn( x) się zerują. Analogicznie możemy wyznaczyć a 1 :
f ( x 1) − f ( x 0)
a 1 =
x 1 − x 0
Ogólnie, oznaczając iloraz różnicowy jako f [ xi], otrzymujemy : f [ xi+1] − f [ xi]
f [ xi, xi+1] =
xi+1 − xi
Gdy znamy ( k − 1)-wszy iloraz różnicowy to możemy wyznaczyć k-ty z : f [ xi+1 , xi+2 , . . . , xi+ k] − f [ xi, xi+1 , . . . , xi+ k− 1]
f [ xi, x − i + 1 , . . . , xi+ k] =
xi+ k − xi
Tak więc wielomian Pn( x) można wyrazić jako : n
X
Pn( x) = f [ x 0] +
f [ x 0 , x 1 , . . . , xk]( x − x 0) . . . ( x − xk− 1) k=1
Powyższy wzór nosi nazwę wzoru interpolacyjnego Newtona.
Interpolacja Hermite’a
Przedstawione do tej pory metody interpolacji wymagają jedynie znajomości węzłów i odpowiednich wartości w tych węzłach. Nie mogą one zatem w żaden sposób oddać kształtu nieznanej zależności. Metodą, która to umożliwia jest interpolacja Hermite’a.
Jeżeli f ∈ C 1[ a, b] i x 0 , . . . , xn ∈ [ a, b] są różne to wielomianem najniższego stopnia zgodnym z f i f 0 w x 0 , . . . , xn jest wielomian stopnia co najwyżej 2 n + 1 dany wzorem : n
n
X
X
H 2 n+1( x) =
f ( xj) Hn,jx +
f 0( xj) ˆ
Hn,jx,
j=0
j=0
gdzie
Hn,j( x) = [1 − 2( x − xj) L0 ( x ( x)
(1)
n,j
j )] L 2
n,j
i
ˆ
Hn,j( x) = ( x − xj) L 2 ( x) (2)
n,j
Ln,j jest j-tym współczynnikiem wielomianu Lagrange’a stopnia n. Ponadto, jeżeli f ∈ C 2 n+2[ a, b] to ( x − x 0)2 . . . ( x − xn)2
f ( x) − H 2 n+1( x) =
f 2 n+2( ξ) ,
(2 n + 2)!
gdzie a < ξ < b
Aby pokazać, że :
H 2 n+1( xk) = f ( xk) i
H0
( x
2 n+1
k ) = f 0( xk )
dla każdego k = 0 , 1 , . . . , n wystarczy pokazać, że Hn,j i ˆ
Hn,j, zdefiniowane w (1) i (2) spełniają warunki :
0 ,
dla j 6= k,
a)
Hn,j( xk) =
b)
d H
1 ,
dla j = k
dx
n,j ( xk ) = 0 , dla wszystkich k
0 ,
jeżeli j 6= k,
c)
ˆ
H
ˆ
n,j ( xk ) = 0 dla wszystkich k
d)
d H
dx
n,j ( xk ) =
1 ,
jeżeli j = k
Rozważając najpierw wielomian ˆ
Hn,j, z warunków c) i ( d) otrzymujemy, że ˆ
Hn,j musi mieć podwójny pierwiastek
w xk dla k 6= j i pojedynczy pierwiastek w xj. Wielomianem stopnia co najwyżej (2 n + 1) spełniającym te warunki i ponadto posiadającym pochodną o odpowiedniej wartości jest : ˆ
( x − x 0)2 . . . ( x − xj− 1)2( x − xj)( x − xj+1)2 . . . ( x − xn)2
Hn,j( x) =
= ( x − xj) L 2 ( x)
( x
n,j
j − x 0)2 . . . ( xj − xj− 1)2(1)( xj − xj+1)2 . . . ( xj − xn)2
Warunki a) i b) implikują, że xk, dla każdego k 6= j, musi być podwójnym pierwiastkiem Hn,j( x). Wielomianem stopnia co najwyżej (2 n + 1) spełniającym a) i b) jest wielomian : Hn,j( x) = ( x − x 0)2 . . . ( x − xj− 1)2( x − xj+1)2 . . . ( x − xn)2(ˆ
ax + ˆ
b) ,
gdzie ˆ
a i ˆ
b są stałymi. Podstawiając :
n
Y
Y
a = ˆ
a
( xi − xj)2
i
b = ˆ
b
i=0
i=0
i6= j
i6= jn ( xi−xj )2
otrzymujemy
Hn,j( x) = L 2 ( x)( ax + b) .
n,j
Warunek a) implikuje, że
1 = Hn,j( xj) = L 2 ( x
n,j
j )( axj + b) = axj + b
(3)
i korzystając z warunku b) oraz równania (3) dostajemy dHn,j( xj)
0 =
=
2 Ln,j( xj) L0 ( xj)( axj + b) + L 2 ( xj)( a) dx
n,j
n,j
=
2 L0
( x
n,j
j )( axj + b) + a
=
2 L0
( x
n,j
j )(1) + a
lub a = − 2 L0
( x
n,j
j ). Z równania (3),
b = 1 − axj = 1 + 2 L0 ( x
n,j
j ) · xj
mamy
ax + b = − 2 L0
( x
( x
( x
n,j
j ) x + 1 + 2 L0n,j
j ) xj = 1 − 2( x − xj ) L0n,j j )
i
Hn,j( x) = ( ax + b) L 2 ( x) = [1 − 2( x − x ( x
( x).
n,j
j ) L0n,j
j )] L 2
n,j
Otrzymaliśmy zatem równania (1) i (2).
Przykład 2
x
0
2
3
y
0
2
5
Wielomian Hermite’a dla :
ma postać :
y’
0
4
6
y”
-
3
-
H 6( x) = 0 . 879 x 2(9 . 849 − 6 . 194 x + x 2)(2 . 443 − 3 . 027 z + z 2) 5
4
3
2
1
0.5
1
1.5
2
2.5
3
Interpolacja funkcjami sklejanymi
Oscylacyjna natura wielomianów interpolacyjnych wysokiego stopnia może prowadzić do znacznych fluk-tuacji na duzych przedziałach zmienności argumentu. Sposobem na zaradzenie temu jest podzielenie naszego przedzału na podprzedziały i okreslenie wielomianu interpolacyjnego, dla kazdego z nich. W ten sposób otrzymujemy metodę zwaną interpolacją za pomocą funkcji sklejanych. Niech będzie dana funkcja f : [ a, b] → R.
Przedział [ a, b] dzielimy na n podprzedziałów : a = x 0 < x 1 < . . . < xn = b i przyjmujemy, że w tych punktach (węzłach) znane są wartości funkcji f ( xi) = yi. Funkcja sklejana S trzeciego stopnia dla funkcji f to funkcja S : R → R ∈ C 2( R) spełniająca następujące warunki : a) S jest wielomianem trzeciego stopnia na podprzedziale [ xj, xj+1] dla każdego j = 0 , 1 , . . . , n − 1, b) S( xj) = f ( xj) dla każdego j = 0 , 1 , . . . , n c) Sj+1( xj+1) = Sj( xj+1) dla każdego j = 0 , 1 , . . . , n − 2
d) S0
( x
( x
j+1
j+1) = S0j
j+1) dla każdego j = 0 , 1 , . . . , n − 2
e) S00
( x
( x
j+1
j+1) = S00
j
j+1) dla każdego j = 0 , 1 , . . . , n − 2
f) jeden z poniższych warunków brzegowych jest spełniony
(i) S00( x 0) = S00( xn) = 0
(swobodne warunki brzegowe)
(ii) S0( x 0) = f 0( x 0) i S0( xn) = f 0( xn) (sztywne warunki brzegowe)
Można pokazać, że dla x ∈ [ xi, xi+1], gdzie i = 0 , 1 , . . . , n − 1 sklejkę naturalną (swobodne warunki brzegowe) można przedstawić w postaci :
S( x) = yi + bi( x − xi) + ci( x − xi)2 + di( x − xi)3
gdzie Sj( xj) = f ( xj) = a 0 oraz na podstawie c) : aj+1 = Sj+1( xj+1)
=
Sj( xj+1)
=
aj + bj( xj+1 − xj) + cj( xj+1 − xj)2 + dj( xj+1 − xj)3
(4)
dla każdego j = 0 , 1 , . . . , n − 2. Przyjmijmy zapis hj = xj+1 − xj oraz an = f ( xn). Równanie (4) przyjmuje postać :
aj+1 = aj + bjhj + cjh 2 + d j = 0 , 1 , . . . , n − 1
(5)
j
j h 3
j
podobnie, gdy S0 ( x
j
j ) = bj dostajemy :
bj+1 = bj + 2 cjhj + 3 djh 2
j = 0 , 1 , . . . , n − 1
(6)
j
gdy cn = S00( xn) to :
2
cj+1 = cj + 3 djhj
j = 0 , 1 , . . . , n − 1
(7)
wyliczjąc z (7) dj i wstawiając do (5) i (6) dostajemy : h 2 j
aj+1 = aj + bjhj +
(2 cj + cj+1)
(8)
3
i
bj+1 = bj + hj( cj + cj+1) (9)
1
hj
bj =
( aj+1 − aj) −
(2 cj + cj+1)
(10)
hj
3
redukując indeks o 1 otrzymujemy :
1
hj− 1
bj− 1 =
( aj − aj− 1) −
(2 cj− 1 + cj)
(11)
hj− 1
3
wstawiając powyższe wartości do (9) (po zmniejszeniu indeksu) : 3
3
hj− 1 cj− 1 + 2( hj + hj) cj + hjcj+1 =
( aj+1 − aj) −
( aj − aj− 1)
(12)
hj
hj− 1
Jedyną niewiadomą w tym układzie równań jest {cj}n . Możemy go przepisać w postaci Ax = b, lub macier-j=0
zowo :
1
0
0
. . .
. . .
0
.
.
h
. .
..
0
2( h 0 + h 1)
h 1
.
.
0
h
. .
..
A =
1
2( h 1 + h 2)
h 2
.
.
.
.
.
.
. .
. .
. .
. .
.
0
.
.
. ..
.
hn− 2
2( hn− 2 + hn− 1)
hn− 1
0
. . .
. . .
0
0
1
0
3
c
( a
( a
0
h
2 − a 1) − 3
1 − a 0)
1
h 0
.
c 1
b =
.
i
x = .
.
.
3
.
( a
( a
h
n − an− 1) −
3
n− 1 − an− 2)
n− 1
hn− 2
cn
0
Przykład 3
Funkcje sklejane dla
krzywa 1
krzywa 2
krzywa 3
i x
7
i
yi
i xi yi
i xi yi
0, 1, 3.0
0, 17, 4.5
0, 27.7, 4.1
6
1, 2, 3.7
1, 20, 7.0
1, 28, 4.3
5
2, 5, 3.9
2, 23, 6.1
2, 29, 4.1
4
3, 6, 4.2
3, 24, 5.6
3, 30, 3.0
5
10
15
20
25
30
4, 7, 5.7
4, 25, 5.8
5, 8, 6.6
5, 27, 5.2
6, 10, 7.1
6, 27.7, 4.1
7, 13, 6.7
8, 17, 4.5