Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne wykład 3
Interpolacja
funkcja przybliżana f ( x ),
siatka węzłów xi , i = 0,...,n, fi = f ( xi )
= =
= =
= =
Dla dowolnych, różnych n+1 punktów węzłowych istnieje
dokładnie jeden wielomian interpolacyjny stopnia, co
najwyżej n taki, że
P( xi ) = fi dla i=0,1,...,n
=
=
=
Wzór interpolacyjny Lagrange a
n n
x - xk
-
-
-
P( x ) = fi
=
=
=
" "
" "
" "
" "
-
i =0 k =0 - xk
= = xi -
= = -
= =
k `"
`"i
`"
`"
W3-1
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne wykład 3
cn , cn-1,Lc1, c0 to możemy obliczyć jego
Jeśli wielomian P(x) ma współczynniki
P( x0 ), P( x1 ),LP( xm )w punktach x0 , x1,Lxm :
wartości
n n-1
P( x0 ) îÅ‚ Å‚Å‚ cn
îÅ‚ Å‚Å‚ x0 x0 L x0 1 îÅ‚ Å‚Å‚
ïÅ‚
ïÅ‚ śł
n n-1
P( x1 )
x1 x1 L x1 1śłïÅ‚ śł
ïÅ‚ śł
ïÅ‚ śł ïÅ‚cn-1śł
ïÅ‚
ïÅ‚ M śł = ïÅ‚ M śł
M M O M Mśł
ïÅ‚ śł
ïÅ‚P( xm-1 )śł ïÅ‚ śł
n n
c1
xm-1 L xm-1 1śłïÅ‚ śł
m-1 -1
ïÅ‚x
ïÅ‚ śł
n n
ïÅ‚
ïÅ‚ śł
P( xm ) c0
xm xm-1 L xm 1śłïÅ‚ śł
ðÅ‚ ûÅ‚ ðÅ‚ ûÅ‚
ðÅ‚ ûÅ‚
Schemat Hornera:
n=3
( ) (( ) )
P( x ) = c3x3 + c2x2 + c1x + c0 c3x + c2 x2 + c1x + c0 c3x + c2 x + c1 x + c0
= =
więc:
c2 c1 c0
c3= a3 a3x a2x a1x
a2=c2+a3x a1=c1+a2x P(x)=c0+a1x
W3-2
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne wykład 3
Współczynniki wielomianu interpolacyjnego
P( x ) = cn xn + cn-1xn-1 +L+ c1x + c0
można obliczyć z:
n n-1
îÅ‚ Å‚Å‚ cn f0
x0 x0 L x0 1 îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚
ïÅ‚
n n-1
f1
x1 x1 L x1 1śłïÅ‚ śł ïÅ‚ śł
ïÅ‚ śł
ïÅ‚cn-1śł ïÅ‚ śł
ïÅ‚
ïÅ‚ M śł ïÅ‚ M śł
=
M M O M Mśł
ïÅ‚ śł
ïÅ‚ śł ïÅ‚
n n-1
macierz Vandermonde a,
c1 fn-1śł
xn-1 L xn-1 1śłïÅ‚ śł ïÅ‚ śł
n-1
ïÅ‚x
n n-1
ïÅ‚
c0 fn
xn xn L xn 1śłïÅ‚ śł ïÅ‚ śł
ðÅ‚ ûÅ‚ ðÅ‚ ûÅ‚
ðÅ‚ ûÅ‚
jest nieosobliwa jeśli węzły xi są różne, ale z
le uwarunkowana (trudno jÄ…
odwrócić)
W3-3
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne wykład 3
Interpolacja przez rodzinę trójkątną
Õ0( x ) = 1
Õ1( x ) = ( x - x0 )
Õ2( x ) = ( x - x0 )( x - x1 )
P( x ) = cnÕn( x )+ cn-1Õn-1( x )+L+ c1Õ1( x )+ c0
,
L
Õn( x ) = ( x - x0 )( x - x1 )L( x - xn-1 )
f0 = P( x0 ) = c0 Ò! c0 = f0
f1 - c0
f1 = P( x1 ) = c1( x1 - x0 ) + c0 Ò! c1 =
x1 - x0
f2 = P( x2 ) = c2( x2 - x0 )( x2 - x1 )+ c1( x2 - x0 )+ c0 Ò! c2 = L
& & & & ..
W3-4
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne wykład 3
Rekurencyjne tworzenie wielomianów interpolacyjnych
Pi (x)
Niech będzie wielomianem stopnia nie
0,i1,...,i
k
większego od k, spełniającym równania węzłów i0 ,i1 ,...,ik:
Pi (xi ) = fi j = 0,...,k
j j
0,i1,...,i
k
Wtedy zachodzi wzór rekurencyjny
Pi(x) = fi i = 0,...,n
(x - xi )Pi (x) - (x - xi )Pi (x)
0 1,i2,...,i 0k 0,i1,...,i
k k-1
Pi (x) =
xi - xi
0,i1,...,ik
k 0
W3-5
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne wykład 3
Metoda Aitken a
x0 P0 ( x ) = f0
=
=
=
x1 P1( x ) = f1 P0,1( x )
=
=
=
x2 P2( x ) = f2 P0,2( x ) P0,1,2( x )
=
=
=
x3 P3 ( x ) = f3 P0,3( x ) P0,1,3( x ) P0,1,2,3( x )
=
=
=
x4 P4 ( x ) = f4 P0,4( x ) P0,1,4( x ) P0,1,2 ,4 ( x ) P0 ,1,2 ,3,4 ( x )
=
=
=
W3-6
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne wykład 3
i xi xi - x yi = Pi (x) P0,i (x) P0,1,i (x) P0,1,2,i (x) K P0,1,K,m (x)
0 x0 x0 - x y0 = P0(x)
x0 - x P0(x)
x1 - x P1(x)
1 x1 x1 - x y1 = P1(x) P0,1(x) =
x0 - x1
x1 - x P0,1(x)
x0 - x P0(x)
x2 - x P0,2(x)
x2 - x P2(x)
2 x2 x2 - x y2 = P2(x) P0,2(x) = P0,1,2(x) =
x0 - x2 x1 - x2
x1 - x P0,1(x) x2 - x P0,1,2(x)
x0 - x P0(x)
x3 - x P0,3(x) x3 - x P0,1,3(x)
x3 - x P3(x)
3 x3 x3 - x y3 = P3(x) P0,3(x) = P0,1,3(x) = P0,1,2,3(x) =
x0 - x3 x1 - x3 x2 - x3
M M M M M M M O
x1 - x P0,1(x) x2 - x P0,1,2(x) xm-1 - x P0,1,K,m-1(x)
x0 - x P0(x)
xm - x P0,m(x) xm - x P0,1,m (x) xm - x P0,1,K,m-2,m (x)
xm - x Pm (x)
m xm xm - x ym = Pm(x) P0,m (x) = P0,1,m (x) = P0,1,2,m(x) = K P0,1,K,m (x) =
x0 - xm x1 - xm x2 - xm xm-1 - xm
W3-7
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne wykład 3
Reszta wzoru interpolacyjnego:
Jeżeli funkcja f (Å" ) ma ciÄ…gÅ‚e pochodne do rzÄ™du n+1 a
Å"
Å"
Å"
P(Å" ) jest wielomianem interpolacyjnym stopnia n, to
Å"
Å"
Å"
n
1
( n+1 )
+
+
+
f ( x ) - P( x ) = f (¾ ) ( x - xi )
- = ¾ -
- = ¾ -
- = ¾ -
"
"
"
"
( n + 1 )! i =0
+ =
+ =
+ =
gdzie ¾ jest pewnym punktem z najmniejszego przedziaÅ‚u
¾
¾
¾
domkniętego zawierającego x, x0 ,..., xn
W3-8
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne wykład 3
1
Przykład: y( x ) =
=
=
=
1 + ( 5x )2
+
+
+
węzły równoodległe w [-1,1] węzły Czebyszewa w [-1,1]
w=[];x=[];y=[];apr=[]; w=[];x=[];y=[];apr=[];
xx=-1:.01:1;yy=1./(1+(5*xx).^2); xx=-1:.01:1;yy=1./(1+(5*xx).^2);
for n=4:16 for n=4:16
h=2/n;
for i=1:n+1 for i=1:n+1
x(n,i)=-1+(i-1)*h; x(n,1:n+1)=-seqcheb(n+1,2);
end end
y(n,1:n+1)=1./(1+(5*x(n,1:n+1)).^2); y(n,1:n+1)=1./(1+(5*x(n,1:n+1)).^2);
w=polyfit(x(n,1:n+1),y(n,1:n+1),n); w=polyfit(x(n,1:n+1),y(n,1:n+1),n);
apr(n,:)=polyval(w,xx); apr(n,:)=polyval(w,xx);
end end
W3-9
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne wykład 3
1
0.8
0.6
0.4
0.2
0
-0.2
-1 -0.5 0 0.5 1
n=5,6,7
W3-10
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne wykład 3
2
1
0
-1
-2
-3
-4
-1 -0.5 0 0.5 1
n=8,9,10,11,12
W3-11
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne wykład 3
1
0.8
0.6
0.4
0.2
0
-0.2
-1 -0.5 0 0.5 1
n=5,6,7
W3-12
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne wykład 3
1.2
1
0.8
0.6
0.4
0.2
0
-1 -0.5 0 0.5 1
n=8,9,10,11,12
W3-13