KATEDRA URZĄDZEŃ I SYSTEMÓW AUTOMATYKI
LABORATORIUM
METOD NUMERYCZYCH
W TECHNICE
INTERPOLACJA
INSTRUKCJA LABORATORYJNA DO ĆWICZENIA NR 1.
WERSJA ROBOCZA
OPRACOWAŁ:
mgr inż. Tomasz KWAŚNIEWSKI
POLITECHNIKA ŚWIĘTOKRZYSKA
KIELCE 2011
Załóżmy że dane są wartości funkcji f ( x) na zbiorze punktów x , x , ..., x zwanych 0
1
n
wę złami interpolacji. Zadaniem interpolacji jest wyznaczenie przybliżonych wartości funkcji f ( x) zwanej funkcją interpolowaną w punktach nie będących węzłami interpolacji.
Przybliżoną wartość funkcji f ( x) obliczamy za pomocą funkcji F ( x) zwaną funkcją
interpolują cą, która w węzłach ma te same wartości co funkcja interpolowana. Funkcja interpolująca jest funkcją pewnej klasy. Najczęściej będzie to wielomian algebraiczny, wielomian trygonometryczny, funkcja wymierna lub funkcja sklejana. Interpolację stosuje się najczęściej gdy nie znamy analitycznej postaci funkcji f ( x) (jest ona tylko stablicowana) lub gdy jej postać analityczna jest zbyt skomplikowana.
2. Interpolacja wielomianowa, wzory Lagrange’a
Dane są węzły interpolacyjne x , x , ..., x , parami różne tzn. x ≠ x ⇔ i ≠ j . Dane 0
1
n
i
j
są wartości funkcji interpolowanej w węzłach f , f , ..., f gdzie f = f ( x ) i = , 1 ,
2 ..., n .
0
1
n
i
i
Zadanie interpolacji polega na znalezieniu wielomianu L ( x) stopnia co najwyżej n, by n
wartości tego wielomianu i funkcji interpolowanej w węzłach były sobie równe czyli L ( x ) = f
i = ,
1 ,
2 ..., .
n
[1]
n
i
i
Twierdzenie 1. Zadanie interpolacji wielomianowej posiada jednoznaczne rozwiązanie, czyli istnieje tylko jeden wielomian spełniający warunek [1].
n
( x− xj)
Konstruujemy funkcje pomocnicze: l x
i ( ) = ∏
j =0 ( x
x
i −
j )
j ≠ i
1 dla i = k
funkcje te są wielomianami stopnia n takimi, że l
i ( xk ) =
0 dla i ≠ k
Stąd wynika, że
n
n
n
( x− xj)
L( x) = ∑ f ( x l⋅ x = ∑ f x
⇒ L x
∏
= f x
[2]
i ) i ( )
( i )
=0
=0
=0
−
i
i
j
( xi xj) ( i) ( i)
j ≠ i
−
K
K
1 ⋅
− 2 ⋅ ⋅ −
− 0 ⋅ − 2 ⋅ ⋅ −
L( x)
( x x ) ( x x )
( x x
x
x
x
x
x
x
N )
(
) (
)
(
N )
= (
⋅ 0 +
⋅ 1 +
x
K
K
0 − x 1 ) ⋅ ( x 0 − x 2 ) ⋅
⋅ (
y
y
x 0 − x
1 −
0 ⋅
1 −
2 ⋅
⋅ 1 −
N )
( x x ) ( x x )
( x xN )
(
[3]
x − x
K
0 ) ⋅ ( x − x 1 ) ⋅
⋅ ( x − xN−1)
+ K + (
⋅ y
x
− x
K
0 ⋅ x
− x 1 ⋅ ⋅ x − x
N
) ( N
)
( N
−1 )
N
N
Wzór [2] i [3] noszą nazwę wzoru interpolacyjnego Lagrange’a.
2
3. Przykład zadania interpolacji wielomianowej
Stosując metody Lagrange’a zbudować wielomian interpolacyjny 4-go stopnia dla następującej tablicy .
x
0 2 3 5 6
f(x) 1 3 2 7 17
Szukamy wielomianu Lagrange’a stopnia o jeden mniejszego niż liczba węzłów interpolacji.
W tym przypadku będzie to wielomian czwartego stopnia w postaci:
L ( x)
2
3
4
= a + a x + a x + a x + a x
1
2
3
4
5
Następnie tworzymy wielomian Lagrange’a korzystając z wzoru [3].
x − 2 ⋅ x − 3 ⋅ x − 5 ⋅ x − 6
x − 0 ⋅ x − 3 ⋅ x − 5 ⋅ x − 6
L
x =
⋅1+
⋅3 +
4 (
) (
) (
) (
) (
)
(
) (
) (
) (
)
(0 2) (0 3) (0 5) (0 6) (2 0) (2 3) (2 5) (2 6)
− ⋅ − ⋅ − ⋅ −
− ⋅ − ⋅ − ⋅ −
( x − 0)⋅( x − 2)⋅( x −5)⋅( x − 6)
( x − 0)⋅( x − 2)⋅( x −3)⋅( x − 6)
+ ( ) ( ) ( ) ( ) ⋅2+
(
) (
) (
) (
) ⋅7 +
3 0
3 2
3 5
3 6
5 0
5 2
5 3
5 6
− ⋅ − ⋅ − ⋅ −
− ⋅ − ⋅ − ⋅ −
( x − 0)⋅( x − 2)⋅( x −3)⋅( x −5)
+ (
⋅17
6 − 0)⋅(6 − 2)⋅(6 − 3)⋅(6 − 5)
Po
wykonaniu
prostych
przekształceń
matematycznych
otrzymujemy
wielomian
interpolacyjny Lagrange’a czwartego stopnia w postaci:
L ( x)
47
481
19
1
2
3
4
= 1+
⋅ x −
⋅ x +
⋅ x −
⋅ x
4
10
180
45
180
4. Skrypt wielomianu interpolacyjnego Lagrange’a.
function [q, L] = lagrange(x, y)
%Input: x = [x0 x1 ... xN]
% y = [y0 y1 ... yN]
%Output: q = współczynniki wielomianu Lagrange'a stopnia N
% L = współczynniki Lagrange'a
N = length(x) - 1;
q = 0;
for m = 1:N + 1
P = 1;
for k = 1:N + 1
if k ~= m,
P = conv(P, [1 -x(k)])/(x(m) - x(k));
end
end
L(m, :) = P; % wspolczynniki wielomianowe Lagrange'a
q = q + y(m)*P; % wielomian Lagrange'a
end
3
5. Wzór interpolacyjny Newtona
Iloraz różnicowy
Funkcja f( x) jest określona za pomocą tablicy: x0, x1, …, xn są węzłami interpolacji, a f( x 0), f( x 1), …, f( xn) – odpowiadającymi tym węzłom wartościami funkcji f( x). Oczywiście xi ≠ xj dla i ≠ j ( i, j = 0, 1, …, n) i ponadto różnice ∆ xi = xi+ 1 – xi ( i = 0, 1, 2, …) nie są na ogół
stałe.
Wyrażenia:
f [
f x − f x
x , x =
0
1]
( 1) ( 0)
x − x
1
0
f [
f x
− f x
x , x
=
1
2 ]
( 2) ( 1)
x − x
2
1
K
K
K
K
K
K
K
K
K
f [
f x
f x
n −
n −
xn
xn =
1,
−
] ( ) ( 1)
x
x
n −
n 1
−
nazywamy ilorazami róż nicowymi pierwszego rzę du. Analogicznie definiujemy ilorazy róż nicowe drugiego rzę du:
f [
f x , x
− f x , x
x , x , x
=
0
1
2 ]
( 1 2) ( 0 1)
x − x
2
0
f [
f x , x
− f x , x
x , x , x
=
1
2
3 ]
( 2 3) ( 1 2)
x − x
3
1
K
K
K
K
K
K
K
K
K
f [
f x
, x
f x
x
n
n −
,
x
, x
,
−
n −
n −
x
n −
n
n =
2
1
−
] ( 1 ) ( 2 1)
x
x
n −
n −2
Ogólnie iloraz róż nicowy rzę du n tworzymy z ilorazu różnicowego n- 1 za pomocą wzoru rekurencyjnego:
+
K
K
1,
+2 ,
, + −
,
,
,
f [
f x
x
x
f x x
x
1
1
,
+ K
[4]
1,
, + =
+
+ −
i
x
i
x
i
x n ]
( i i
i n )
( i i
i n
)
x + − x
i n
i
dla n = 1, 2, … oraz i = 0, 1, 2, …
Z ilorazów różnicowych tworzy się tablice – oto zasada tworzenia:
xi f(xi)
Ilorazy różnicowe
rzędu 1
rzędu 2
rzędu 3
rzędu 4
x0 f(x0)
f(x0, x1)
x1 f(x1)
f(x0, x1, x2)
f(x1, x2)
f(x0, x1, x2, x3)
x2 f(x2)
f(x1, x2, x3)
f(x0, x1, x2, x3, x4)
f(x2, x3)
f(x1, x2, x3, x4)
x3 f(x3)
f(x2, x3, x4)
f(x3, x4)
x4 f(x4)
4
Ilorazem różnicowym rzędu k funkcji f opartym na parami różnych węzłach xl, xl+1, …, xl+k, w których określona jest funkcja f nazywamy wyrażenie:
l + k
f ( xi )
[ x , x
,K, x
; f ]
[5]
l
l +1
l + k
= ∑ l+ k
i = l
∏( x x )
i −
j
j = l j ≠ i
przy czym przez iloraz różnicowy rzędu zerowego nazywamy wartość tej funkcji w danym węźle: [ xl; f] = f( xl).
Własności ilorazu różnicowego
1. Iloraz różnicowy jest funkcją symetryczną, czyli
[ x , x
,K, x
; f ]
+
+
= [ x , x ,K, x ; ]
l
l 1
f
l k
i
i
i
l
l 1
+
l + k
gdzie x , x ,K x
jest dowolną permutacją liczb l, l+ 1, …, l+k
l
l +
,
1
l + k
2. Iloraz różnicowy jest funkcjonałem liniowym ze względu na funkcję f, tzn. jeśli f( x) = g( x) + h( x), to
[ x , x
,K, x
; f ] = [ x , x
,K, x
; g] [ x , x
,K, x
; h]
l
l 1
+
+
l + k
l
l 1
+
l + k
l
l 1
+
l + k
3. Iloraz różnicowy dla jednostajnej siatki ρ n (równoodległe węzły: xk+ 1 = x 0 + ( k + 1) h dla k = 0, 1, …, n–1, gdzie h jest stałą nazywaną długością kroku):
[
∆ f x
ρ ; f
0
=
n
] n ( )
n
n! h
gdzie n
∆ oznacza różnicę progresywną:
n
∆ f ( x = f x
0 )
( 0)
1
∆ f ( x = f
∆ x = f x + h − f x
0 )
( 0) ( 0 ) ( 0)
2
∆ f ( x = ∆ f
∆ x = ∆ f x + h − f x = f x + 2 h − f x + h − f x + h − f x 0 )
( ( 0)
( ( 0 ) ( 0)
( 0
) ( 0 ) ( ( 0 ) ( 0)
k
−
∆ f ( x ) = ∆( k 1
∆ f x
0
( 0))
Wzór interpolacyjny Newtona
Jeżeli teraz podstawimy wzór (5) do wzoru interpolacyjnego Lagrange’a przyjmując, że l = 0
to otrzymamy nowy wzór nazywany wzorem interpolacyjnym Newtona: N
=
+
,
ω
+
,
,
ω
+K+
,
,K,
ω
[6]
n ( x)
f ( x 0 ) f ( x x
0
1 ) 0 ( x)
f ( x x x
0
1
2 ) 1( x)
f ( x x
x
0
1
n ) n 1
− ( x)
gdzie ω x = x − x x
[7]
0
− x K x
1
− x
n ( )
(
)(
) (
n )
Niech ρ n oznacza dowolną siatkę bez węzłów wielokrotnych. Wielomianem interpolacyjnym Newtona Nn dla funkcji f na siatce ρ n nazywamy wielomian: N
K
n ( x) = f ( x
+ x x f x − x + x x x f x − x x − x + +
0 )
[ 0, 1; ](
0 )
[ 0, 1, 2; ](
0 )(
1 )
[8]
+ [ x x K x
K
n
f x − x
x − x
x −
0 , 1,
,
; ](
0 )(
1 )
( xn 1−)
5
6. Skrypt wielomianu interpolacyjnego Newtona.
function [n, DD] = newton(x, y)
%Input: x = [x0 x1 ... xN]
% y = [y0 y1 ... yN]
%Output: n = współczynniki wielomianu Newtona
N = length(x) - 1;
DD = zeros(N + 1, N + 1); %w pierwsza kolumnę macierzy kwadratowej N+1 x
N+1
DD(1:N + 1,1) = y';
for k = 2:N + 1
for m = 1:N + 2 - k % tablica różnic dzielonych
DD(m, k) = (DD(m + 1, k - 1) - DD(m, k - 1))/(x(m + k - 1) - x(m));
end
end
% określenie współczynników wielomianu interpolacyjnego Newtona
a = DD(1, :); % równanie (4) instrukcja laboratoryjna
n = a(N + 1);
for k = N:-1:1
n = [n a(k)] - [0 n*x(k)];
end
5. Przykład zadania interpolacji metodą Newtona
Korzystając ze wzoru interpolacyjnego Newtona, znaleźć wielomian interpolacyjny dla następującej tablicy .
x
0 2 3 5 6
f(x) 1 3 2 7 17
xi f(xi)
Ilorazy różnicowe
rzędu 1 rzędu 2 rzędu 3 rzędu 4
0 1
1
2 3
-2/3
-1
11/30
3 2
35/30
-1/180
2,5
1/3
5 7
2,5
10
6 17
2
11
1
N x = 1 + 1 x − 0 −
x − 0 x − 2 +
x − 0 x − 2 x − 3 −
x − 0 x − 2 x − 3 x − 5
4 ( )
(
) (
)(
)
(
)(
)(
)
(
)(
)(
)(
)
3
30
180
47
481 2
19 3
1
4
N ( x) = 1 +
x −
x +
x −
x
4
10
30
45
180
6
6. Błąd interpolacji wielomianowej i optymalny dobór węzłów Jak zaznaczono we wstępie, za pomocą wielomianu interpolacyjnego można wyznaczyć przybliżoną wartość funkcji interpolowanej w punktach nie będących węzłami.
Jeśli funkcja interpolowana f ( x) jest klasy Cn+1 w przedziale domkniętym 〈 a, b〉 oraz wszystkie węzły x , x ,..., x też należą do tego przedziału to błąd bezwzględny interpolacji 0
1
n
wielomianem Lagrange’a można oszacować wzorem
M
R ( x) = f ( x) − L ( x)
n 1
+
≤
ω ( x)
[9]
n
n
+
( n + )
1 !
n 1
gdzie
M
= max ( n+ )1
f
( x) a ω
( x) jest wielomianem czynnikowym określonym wzorem [7].
n 1
+
n 1
+
x 〈
∈ a, b〉
Zauważmy, że błąd interpolacji zależy nie tylko od (n+1)-szej pochodnej funkcji interpolowanej ale również od wielomianu czynnikowego ω
( x) , który z kolei zależy od
n 1
+
doboru węzłów interpolacji
x , x ,..., x . Problem optymalnego doboru węzłów 0
1
n
interpolacyjnych rozwiązał Czebyszew. Wykazał, że wartości węzłów w przedziale 〈 a, b〉 , optymalnie dobranych są określone
1
1
i
2 + 1
x =
( a + b) + ( b − a) cos(
π ) i = ,1
,
0 .. n
. .
[10]
i
2
2
2 n + 2
Wtedy najmniejsze oszacowanie błędu interpolacji wynosi
n 1
M
( b a) +
n 1
R ( x)
f x
L x
[11]
n
= ( ) − ( )
+
−
n
≤
2 n 1
( n + )
1 !
2 +
Optymalnie dobrane węzły wcale nie są równo odległe lecz zagęszczają się przy końcach przedziału.
7. Iteracyjna metoda Aitkena
Wielomian w postaci wzoru Lagrange’a jest niewygodny zarówno do wyznaczania jego wartości w dowolnym punkcie (stosuje się wzór Aitkena) jak i jego całkowania bądź
różniczkowania. Częściej wielomian interpolacyjny określa się w postaci wzoru Newtona przy czym obydwa wzory są sobie równoważne ponieważ, zgodnie z twierdzeniem 1, istnieje tylko jeden wielomian interpolacyjny dla węzłów x , x ,..., x .
0
1
n
Istnieje metoda obliczania wartości wielomianu Lagrange'a w zadanym punkcie, bez obliczania samego wielomianu interpolacyjnego. Służy do tego iteracyjna metoda Aitkena.
Oznaczmy przez W wielomian który w węzłach
x , x i ≠ j przyjmuje wartości
i
j (
)
i, j
f ( x ), f x :
j
( i )
y
x − x
i
i
y
x − x
j
j
W
=
i, j
x − x
j
i
7
W K
x − x
1,2,
, k 1
− , k
k
W K
x − x
−
W
K
x =
k m (
)
1,2,
, k 1, m
m
1,2,
, ,
x − x
m
k
Aby obliczyć wartość wielomianu interpolacyjnego opartego na n węzłach w dowolnym punkcie a różnym od węzłów, należy obliczyć wartość W
. Wszystkie wyniki
1,2,
,
K n
niezbędnych obliczeń wygodnie jest umieścić w macierzy trójkątnej wraz z węzłami oraz ich wartościami (schemat taki nazywamy schematem Aitkena). Rozwiązanie takie jest dogodne zarówno podczas rachunków ręcznych, jak i maszynowych, gdyż podczas obliczania każdej wartości zawsze korzystamy z wartości położonych na lewo w tym samym rzędzie i powyższych.
x
y
1
1
x
y
W
2
2
1,2
x
y
W
W
3
3
1,3
1,2,3
K K
K
K
x
y
W
W
W
n
n
1, n
1,2, n
1,2, ,
K n
Wyznaczenie f(4) metodą Aitkena :
1
0 − 4
1
0 − 4
1
0 − 4
1
0 − 4
3
2 − 4
2
3 − 4
7
7
5 − 4
29
17
6 − 4
70
W
=
= 5 , W =
= , W =
=
, W
=
=
1,2
2 − 0
1,3
3 − 0
3
1,4
5 − 0
5
1,5
6 − 0
6
5
2 − 4
5
2 − 4
5
2 − 4
7
29
70
3 − 4
5 − 4
6 − 4
3
1
5
83
6
25
W
=
= − , W
=
=
, W
=
=
1,2,3
3 − 2
3
1,2,4
5 − 2
15
1,2,5
6 − 2
3
1
−
1
3 − 4
−
3 − 4
3
3
83
25
5 − 4
6 − 4
15
13
3
23
W
=
=
, W
=
=
1,2,3,4
5 − 3
5
1,2,3,5
6 − 3
9
13
5 − 4
5
23
6 − 4
9
119
W
=
=
.
1,2,3,4,5
6 − 5
45
8
1
2
3
5
7
1
3
2
−
3
3
29
83
13
5
7
5
15
5
70
25
23
119
6 17
6
3
9
45
119
Stąd W
=
, jest wartością funkcji Lagrange’a w punkcie x = 4 .
1,2,3,4,5
45
8. MATLAB - rodzaje interpolacji
W programie MATLAB mamy do dyspozycji kilka metod interpolacji:
- wielomian pierwszego stopnia,
- wielomian trzeciego stopnia,
- metoda najbliższych sąsiadów,
- funkcje sklejane.
Interpolacje stosuje się do zagęszczania tabel. Zabieg ten pozwala na dokładniejsze przybliżenie danych zawartych w tabelach. Zadanie interpolacji możemy podzielić na: Interpolacja nieparametryczna
- algorytm najbliższy sąsiad - nearest
Interpolacja parametryczna
- wielomianowa – linear, cubic, pchip
- funkcjami sklejanymi - spline
Program MATLAB, w podstawowej bibliotece posiada funkcję o nazwie interp1.
Funkcja umożliwia rozwiązanie zadania interpolacji funkcji jednej zmiennej y = f(x) w punktach xi nie będących węzłami interpolacyjnymi.
yi = interp1(x, y, xi, 'metoda')
gdzie:
x, y - wektory współrzędnych węzłów interpolacji,
xi - wektor punktów na osi X dla których będą obliczane interpolowane wartości yi.
9
Zadanie 1.
Znaleźć wielomian interpolacyjny Lagrange’a dla funkcji określonej tablicą wartości oraz wyznaczyć f(7) wykorzystując metodę AITKENA.
x
-2 1 2
4 9
f(x) 3
1 -3 8 15
wyznaczyć f(7) wykorzystując metodę AITKENA
x
-6 -2 2
4
8
f(x) 2
5
10 12 13
wyznaczyć f(-4) wykorzystując metodę AITKENA
x
-5 -3 2 7 12
f(x) -2 -1 0 8 10
wyznaczyć f(5) wykorzystując metodę AITKENA
Zadanie 2.
Dokonać interpolacji funkcji:
a) - f ( x) = co (
s π ⋅ x)
b) - f ( x) = x 3 cos(π ⋅ x)
c) - f ( x) =
(π ⋅ )2
sin
x
w przedziale x ∈ − ,
2 6 z krokiem 0,5 funkcją wbudowaną w Matlab’a interp1 wszystkimi metodami. Narysować wykres danej funkcji i funkcji interpolującej na jednym układzie współrzędnych oraz wykres błędu interpolacji.
Zadanie 3.
Ocenić z jaką dokładnością można obliczyć wartość ln 100,5 przy użyciu wzoru interpolacyjnego Lagrange’a, jeżeli dane są wartości: ln 100, ln 101, ln 102 i ln 103.
Zadanie 4.
π
Z jaką dokładnością można obliczyć wartość sin
stosując wzór interpolacyjny Lagrange’a,
36
π
π
π
jeżeli dane są wartości sin ,
0 sin
,sin
,sin
.
6
4
3
Zadanie 5.
Wykorzystując ilorazy różnicowe podaj wzór ogólny wielomianu interpolacyjnego Newtona dla funkcji określonej tablicą w zadaniu 1.
10
Obliczyć wartość 119 za pomocą wzoru interpolacyjnego Lagrange,a dla funkcji y =
x i
węzłów interpolacji: x = 10 ,
0 x = 12 ,
1 x = 144. Oszacować błąd.
0
1
2
Program ćwiczenia
1.
Napisać program, który oblicza wartości wielomianu Lagrange’a dla dowolnych punktów x, węzłów równoodległych lub dobranych optymalnie [13] i zadanej przez prowadzącego funkcji interpolowanej.
2.
Napisać program, który oblicza wartości wielomianu Newtona dla dowolnych punktów x, węzłów równoodległych lub dobranych optymalnie [13] i zadanej przez prowadzącego funkcji interpolowanej.
3.
Napisać program, który oblicza wartości wielomianu wg schematu Aitkena dla dowolnych punktów x, węzłów równoodległych lub dobranych optymalnie [13] i zadanej przez prowadzącego funkcji interpolowanej.
4.
Sprawozdanie powinno zawierać:
a) kod programu i wyniki przeprowadzonych obliczeń,
b) opis przeprowadzonych badań,
c) analizę uzyskanych wyników,
d) wnioski.
7. Pytania kontrolne
1. Co to jest interpolacja i po co się ją stosuje?
2. Sformułuj zadanie interpolacji wielomianowej.
3. Podaj ogólną postać wzoru Lagrange’a. Znaleźć wielomian interpolacyjny Lagrange’a dla funkcji określonej przez prowadzącego (zastosuj wzór Lagrange’a).
4. Podaj wzór Newtona. Znaleźć wielomian interpolacyjny Newtona dla funkcji określonej przez prowadzącego (zastosuj wzór Newtona).
5. Co to jest iloraz różnicowy stopnia n-tego. Jak wygląda zależność rekurencyjna ilorazu różnicowego?
6. Podaj wzór na błąd interpolacji wielomianowej. Jak wygląda optymalny dobór węzłów interpolacji?
7. Scharakteryzuj algorytm Aitkena, po co się go stosuje?
11