Introduction to Scientific
Computing
-- A Matrix Vector Approach Using Matlab
Written by Charles F.Van Loan
陈 文 斌
Wbchen@fudan.edu.cn
复旦大学
Chapter2
Polynomial Interpolation
? The vandermonde Approach
? The Newton Approach
? Properties
? Special Topics
Polynomial interpolation problem
:n fo r iy)(xps u c h th a t ( o r le s s )
r e e n -( x ) o f y n o m ia l pfi n d a p o l
,,.,,,y) a n d y ( d i s t in c t,.,,,xG iv e n x
iin-
n-
nn
1
1d e g
1
1
11
??
Vandermonder Approach
123211,..)( ?? ????? nnn xaxaxaaxp
i
n
iniiin yxaxaxaaxp ??????
?
?
12
3211,,,)(
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
nn
n
nnn
n
n
n
y
y
y
y
a
a
a
a
xxx
xxx
xxx
xxx
??
?
?????
?
?
?
3
2
1
3
2
1
12
1
3
2
33
1
2
2
22
1
1
2
11
1
1
1
1
Base 1,x,x2,…
V
n = length(x); V = ones(n,n)
for i=1:n
% Set up row i,
for j=1:n
V(i,j) = x(i)^(j-1);
end
end
n = length(x); V = ones(n,n);
for i=1:n
% Set up row i,
V(i,:) = (x(i)*ones(1,n)).^(j-
1);
end
n = length(x); V = ones(n,n)
for j=1:n
% Set up column i,
for i=1:n
V(i,j) = x(i)^(j-1);
end
end
n = length(x); V = ones(n,n)
for j=2:n
% Set up column i,
for i=1:n
V(i,j) = x(i)*V(i,j-1);
end
end
function a = InterpV(x,y)
% a = InverpV(x,y)
% This computes the Vandermonde polynomial interpolant where
% x is a column n-vector with distinct components and y is a
% column n-vector.
%
% a is a column n-vector with the property that if
%
% p(x) = a(1) + a(2)x +,.,a(n)x^(n-1)
% then
% p(x(i)) = y(i),i=1:n
n = length(x);
V = ones(n,n);
for j=2:n
% Set up column j.
V(:,j) = x.*V(:,j-1);
end
a = V\y; 2n3/3 flops
Nested Multipication
zxxaaxp nnn ???? ?? a t,,,)( 111
n=length(a); zpower=1;
pval=a(1);
for i=2:n
zpower=z*zpower;
pval=pval+a(i)*zpower;
end
Horner's ruler
1234
3
4
2
3211
))((
)(
axaxaxa
xaxaxaaxp n
????
?????
n=length(a);
pval=a(n);
for i=n-1:-1:1
pval=z*pval+a(i);
end
m=length(z);n=length(a);
pval=zeros(m,1);
For j=1:m
pval(j)=a(n);
for i=n-1:-1:1
pval(j)=z(j)*pval(j)+a(i);
end;end
Many points
0 2 4 6-2
-1
0
1
2
0 2 4 6-2
-1
0
1
2
0 2 4 6-2
-1
0
1
2
0 2 4 6-2
-1
0
1
2
Evalute at many points
Function pval=HornerV(a,Z)
m=length(z);n=length(a);
pval=a(n)*ones(m,1);
for k=n-1:-1:1
pval=z.*pval+a(k);
end
2m flops
All:2mn flops
Script see next
page
% Script File,ShowV
% Plots 4 random cubic interpolants of sin(x) on [0,2pi].
% Uses the Vandermonde method.
close all
x0 = linspace(0,2*pi,100)';
y0 = sin(x0);
for eg=1:4
x = 2*pi*sort(rand(4,1));
y = sin(x);
a = InterpV(x,y);
pVal = HornerV(a,x0);
subplot(2,2,eg)
plot(x0,y0,x0,pVal,'--',x,y,'*')
axis([0 2*pi -2 2])
end
Display cubic interpolants of
sin(x),The abscissas are chosen
randomly
Newton representation
Four point example
))()((
))((),(,1 b a s i s
321
211
xxxxxx
xxxxxx
???
???
))()((
))(()()(
3214
2131213
xxxxxxc
xxxxcxxccxp
????
??????
))()((
))(()()(
))(()(
)(,
3424144
2414314214
2313313213
1221211
xxxxxxc
xxxxcxxccxy
xxxxcxxccy
xxccycy
????
??????
??????
????
))()((
)))(()((
))((
))((
,
342414
2414314214
4
2313
13213
3
12
12
211
xxxxxx
xxxxcxxccy
c
xxxx
xxccy
c
xx
cy
cyc
???
??????
?
??
???
?
?
?
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??????
???
?
4
3
2
1
4
3
2
1
342414241414
231313
12
))()(())(()(1
0))(()(1
00)(1
0001
y
y
y
y
c
c
c
c
xxxxxxxxxxxx
xxxxxx
xx
Matrix-vector notation
General n Case
)()(...)()( 12232 ???????? nn xxxxcxxccxq ?
)()()( 11 xqxxyxp ???
ni
xx
yyx
i
i
i,2,,
1
1 ?
???
?
???
?
?
?
nccc ?,,32
ncccc ?,,,321
function c = InterpNRecur(x,y)
% c = InterpNRecur(x,y)
% The Newton polynomial interpolant.
% x is a column n-vector with distinct components and y is
% a column n-vector,c is a column n-vector with the property that if
%
% p(x) = c(1) + c(2)(x-x(1))+...+ c(n)(x-x(1))...(x-x(n-1))
% then
% p(x(i)) = y(i),i=1:n.
n = length(x);
c = zeros(n,1);
c(1) = y(1);
if n > 1
c(2:n) = InterpNRecur(x(2:n),(y(2:n)-y(1))./(x(2:n)-x(1)));
end
function c = InterpN(x,y)
% c = InterpN(x,y)
% The Newton polynomial interpolant.
% x is a column n-vector with distinct components and y is a column
% n-vector,c is a column n-vector with the property that if
%
% p(x) = c(1) + c(2)(x-x(1))+...+ c(n)(x-x(1))...(x-x(n-1))
% then
% p(x(i)) = y(i),i=1:n.
n = length(x);
for k = 1:n-1
y(k+1:n) = (y(k+1:n)-y(k)),/ (x(k+1:n) - x(k));
end
c = y;
Nonrecursive
implementation
2
1
1
1
1 2
3
2
)1(33)(3 nnnkkn n
k
n
k
????? ??
?
?
?
?
flops
function pVal = HornerN(c,x,z)
% pVal = HornerN(c,x,z)
% Evaluates the Newton interpolant on z where
% c and x are n-vectors and z is an m-vector.
%
% pVal is a vector the same size as z with the property that if
%
% p(x) = c(1) + c(2)(x-x(1))+,.,+ c(n)(x-x(1))...(x-x(n-1))
% then
% pval(i) = p(z(i)),i=1:m.
n = length(c);
pVal = c(n)*ones(size(z));
for k=n-1:-1:1
pVal = (z-x(k)).*pVal + c(k);
end
Nested Multiplication
Efficiency
? Time efficiency
– Vandermode approach,2n3/3 flops(cubic)
– Newton approach,3/2n2 flops(quadratic)
? Memory efficiency
– InterpV requires n-by-n array
– InterpN needs just a few n-vectors
– InterpNRecur requires a couple of n-vectors
with each level of recursion.
Accuracy
b.ηw h e r e a
)x(x)x(x
n!
( ηf
( x )pf ( x )
I,n y xt h e n f o r a
h e x,n t a i n i n g te r v a l I c o a b l e o n a nd i f f e r e n t i
n t i n u o u s l yn t i m e s c o,I f f i s,.,,,xs xi n tt i n c t p oa t t h e d i s
o n f ( x )t h e f u n c t ie r p o l a t e s ( x ) S u p p o s e P
n
( n )
n-
n
n-
??
????
?
?
11
1
1
)
2 T h e o r e m
int
int
If a function has ill-behaved higher derivatives,
then the quality of the polynomial interpolants
may actually decrease as the degree increases.
for n=7:35
xEqual = linspace(-1,1,n)';
yEqual = ones(n,1)./(1+25*xEqual.^2);
cEqual=InterpN(xEqual,yEqual);
pValsEqual = HornerN(cEqual,xEqual,x);
plot(x,y,x,pValsEqual,xEqual,yEqual,'*')
title(sprintf('Equal Spacing (n = %2.0f)',n))
end
the equal-spacing interpolants
of f(x) = 1/(1+25x^2) on [-1,1]'
The Runge phenomenon
),( 11 xfc ?
Divided differences
34
23
12
12
13
13
24
12
12
14
14
4
)()()()()()()()(
xx
xx
xx
xfxf
xx
xfxf
xx
xx
xfxf
xx
xfxf
c
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
23
12
12
13
13
3
)()()()(
xx
xx
xfxf
xx
xfxf
c
?
?
?
?
?
?
?
12
12
2
)()(
xx
xfxfc
?
??
],.,,,[ 1 kk xxfc ?
??
?
?
??
?
?
?? ??
?
??
?
1
11
11 )(],.,,,[)(
k
j
j
n
k
kn xxxxfxp
l y,r e s p e c t i v e},,.,,,{ a n d },.,,,{on
of tsi n t e r p o l a n t h ea r e )( a n d )( S u p p o s e
211 kk-
RL
xxxx
fxpxp
k
RLk
xx
xpxxxpxx
xp
?
???
?
)()()()(
)( 1
kixfxp ii,1),()( ??
)()](,.,,,[
)](,[][)(
111
1211
????
????
kk xxxxxxf
xxxxfxfxp
?
?
)()](,.,,,[
)](,[][)(
2111
1211
?? ???
????
kk
L
xxxxxxf
xxxxfxfxp
?
?
)()](,.,,,[
)](,[][)(
1212
2322
?? ???
????
kk
R
xxxxxxf
xxxxfxfxp
?
?
1
112
1
],...,[],...,[
],...,[
xx
xxfxxf
xxf
k
kk
k ?
?
? ?
]),,,,([ 54321 xxxxxf
]),,,([ 4321 xxxxf ]),,,([
5432 xxxxf
]),,([ 321 xxxf ]),,([ 432 xxxf
]),([ 21 xxf ]),([ 32 xxf
])([ 1xf ])([ 2xf
])([ 2xf ])([ 3xf
]),,([ 432 xxxf ]),,([
543 xxxf
]),([ 32 xxf ]),([ 43 xxf ]),([ 43 xxf ]),([ 54 xxf
])([ 3xf ])([ 4xf
])([ 4xf ])([ 5xf
]),([ 32 xxf ]),([ 43 xxf
])([ 2xf ])([ 3xf
])([ 2xf
])([ 3xf ])([ 4xf
])([ 3xf ])([ 4xf])([ 3xf
Divided differences
])([ 1xf
])([ 2xf
])([ 3xf
])([ 4xf
])([ 5xf
]),([ 21 xxf
]),([ 32 xxf
]),([ 43 xxf
]),([ 54 xxf
]),,([ 321 xxxf
]),,([ 432 xxxf
]),,([ 543 xxxf
]),,,([ 4321 xxxxf
]),,,([ 5432 xxxxf
]),,,,([ 54321 xxxxxf
Efficient computation of divided differences
function c = InterpN2(x,y)
% c = InterpN2(x,y)
% The Newton polynomial interpolant.
% x is a column n-vector with distinct components and y is
% a column n-vector,c is a column n-vector with the property that if
%
% p(x) = c(1) + c(2)(x-x(1))+...+ c(n)(x-x(1))...(x-x(n-1))
% then
% p(x(i)) = y(i),i=1:n.
n = length(x);
for k = 1:n-1
y(k+1:n) = (y(k+1:n)-y(k:n-1)),/ (x(k+1:n) - x(1:n-k));
end
c = y;
hixx i )1(1 ??? Equally spaced
)1(
],...,[],...,[],...,[ 112
1 ?
?? ?
kh
xxfxxfxxf kk
k
?
?
?
????
???
? 1],.,,,[],.,,,[
1)(],.,,,[
112
1
1
kxxfxxf
kxfxxf
kk
k
differences
)!1(
],.,,,[],.,,,[
1
1
1
?
??
? kh
xxfxxf
k
k
k
123452345345455
1234234344
123233
122
1
464332
332
2
fffffffffffffff
ffffffffff
ffffff
fff
f
??????????
??????
???
?
? ?6
?
n?
???
?
?
??
?
?
]2/[ n
n
differences Matlab diff
d=diff(y) d=y(2:n)-y(1:n-1)
d=diff(y,2) d=y(3:n)-2*y(2:n-1)+y(1:n-1)
x=linspace(0,1,6)';
y=x.*x;
a=InterpV(y,x);
yvals=linspace(y(1),y(6));
xvals=HornerV(a,yvals);
plot(yvals,xvals,'r');
hold on
ezplot('sqrt',[0,1]);
plot(y,x,'*')
2)( xxf ? xxg ?)(
xxfg ?))((
Inverse
interpolation
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
x
sqrt(x)
Interpolation in two dimensions
-0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
(a,c) ((1-m)a+? b,c) (b,c)
(a,d) (b,d)((1-m)a+? b,d)
(b,(1-l)c+? d)(a+(1-l)c+? d)
},:),{( dycbxayxR ?????
bdadxd
bcacxc
fff
fff
ll
ll
???
???
)1(
)1(
)~,~()1( yxfffcz xdxc ???? mm
))1(,)1((
))1(())1)((1(
dcbaf
ffffz bdadbcac
mmll
llmllm
?????
???????
function fA = SetUp(fname,a,b,n,c,d,m)
x = linspace(a,b,n);
y = linspace(c,d,m);
fA = zeros(n,m);
for i=1:n
for j=1:m
fA(i,j) = feval(fname,x(i),y(j));
end
end
function u=fun1(x,y)
u=sin(x+y);
a=setup('fun1',0,10,100,0,10,100)
surf(a)
Test code
get the value of
Function fname
function z = LinInterp2D(xc,yc,a,b,c,d,fA)
[n,m]=size(fA);
% xc=a+(i-1+dx)*hx 0<=dx<=1
hx=(b-a)/(n-1);
i=max([1 ceil((xc-a)/hx)]);
dx=(xc (a+(i-1)*hx))/hx;
hy=(d-c)/(m-1);
j=max([1 ceil((yc-c)/hy)]);
dy=(yc (c+(j-1)*hy))/hy;
z=(1-dy)*((1-dx)*fA(i,j)+dx*fA(i+1,j)) + dy*((1-…
dx)*fA(i,j+1)+dx*fA(i+1,j+1));
CEIL Round towards
plus infinity.
Matlab's Polynomial Tools
close all
x = linspace(0,4*pi);
y = exp(-.2*x).*sin(x);
x0 = [.5 4 8 12];
y0 = exp(-.2*x0).*sin(x0);
p = polyfit(x0,y0,length(x0)-1);
pvals = polyval(p,x);
plot(x,y,x,pvals,x0,y0,'o')
axis([0 4*pi -1 1])
set(gca,'XTick',[])
set(gca,'YTick',[])
title('Interpolating {\ite}^{-.2{\itx}}sin({\itx}) with a Cubic
Polynomial')
polyfit
polyval
Interpolating e-.2xsin(x) with a Cubic Polynomial
POLYFIT Fit polynomial to data.
POLYFIT(X,Y,N) finds the coefficients of a polynomial P(X) of
degree N that fits the data,P(X(I))~=Y(I),in a least-squares sense.
[P,S] = POLYFIT(X,Y,N) returns the polynomial coefficients P and
a structure S for use with POLYVAL to obtain error estimates on
predictions,If the errors in the data,Y,are independent normal
with constant variance,POLYVAL will produce error bounds which
contain at least 50% of the predictions.
The structure S contains the Cholesky factor of the Vandermonde
matrix (R),the degrees of freedom (df),and the norm of the
residuals (normr) as fields,
This centering and scaling transformation improves the numerical
properties of both the polynomial and the fitting algorithm,
Computing
-- A Matrix Vector Approach Using Matlab
Written by Charles F.Van Loan
陈 文 斌
Wbchen@fudan.edu.cn
复旦大学
Chapter2
Polynomial Interpolation
? The vandermonde Approach
? The Newton Approach
? Properties
? Special Topics
Polynomial interpolation problem
:n fo r iy)(xps u c h th a t ( o r le s s )
r e e n -( x ) o f y n o m ia l pfi n d a p o l
,,.,,,y) a n d y ( d i s t in c t,.,,,xG iv e n x
iin-
n-
nn
1
1d e g
1
1
11
??
Vandermonder Approach
123211,..)( ?? ????? nnn xaxaxaaxp
i
n
iniiin yxaxaxaaxp ??????
?
?
12
3211,,,)(
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
nn
n
nnn
n
n
n
y
y
y
y
a
a
a
a
xxx
xxx
xxx
xxx
??
?
?????
?
?
?
3
2
1
3
2
1
12
1
3
2
33
1
2
2
22
1
1
2
11
1
1
1
1
Base 1,x,x2,…
V
n = length(x); V = ones(n,n)
for i=1:n
% Set up row i,
for j=1:n
V(i,j) = x(i)^(j-1);
end
end
n = length(x); V = ones(n,n);
for i=1:n
% Set up row i,
V(i,:) = (x(i)*ones(1,n)).^(j-
1);
end
n = length(x); V = ones(n,n)
for j=1:n
% Set up column i,
for i=1:n
V(i,j) = x(i)^(j-1);
end
end
n = length(x); V = ones(n,n)
for j=2:n
% Set up column i,
for i=1:n
V(i,j) = x(i)*V(i,j-1);
end
end
function a = InterpV(x,y)
% a = InverpV(x,y)
% This computes the Vandermonde polynomial interpolant where
% x is a column n-vector with distinct components and y is a
% column n-vector.
%
% a is a column n-vector with the property that if
%
% p(x) = a(1) + a(2)x +,.,a(n)x^(n-1)
% then
% p(x(i)) = y(i),i=1:n
n = length(x);
V = ones(n,n);
for j=2:n
% Set up column j.
V(:,j) = x.*V(:,j-1);
end
a = V\y; 2n3/3 flops
Nested Multipication
zxxaaxp nnn ???? ?? a t,,,)( 111
n=length(a); zpower=1;
pval=a(1);
for i=2:n
zpower=z*zpower;
pval=pval+a(i)*zpower;
end
Horner's ruler
1234
3
4
2
3211
))((
)(
axaxaxa
xaxaxaaxp n
????
?????
n=length(a);
pval=a(n);
for i=n-1:-1:1
pval=z*pval+a(i);
end
m=length(z);n=length(a);
pval=zeros(m,1);
For j=1:m
pval(j)=a(n);
for i=n-1:-1:1
pval(j)=z(j)*pval(j)+a(i);
end;end
Many points
0 2 4 6-2
-1
0
1
2
0 2 4 6-2
-1
0
1
2
0 2 4 6-2
-1
0
1
2
0 2 4 6-2
-1
0
1
2
Evalute at many points
Function pval=HornerV(a,Z)
m=length(z);n=length(a);
pval=a(n)*ones(m,1);
for k=n-1:-1:1
pval=z.*pval+a(k);
end
2m flops
All:2mn flops
Script see next
page
% Script File,ShowV
% Plots 4 random cubic interpolants of sin(x) on [0,2pi].
% Uses the Vandermonde method.
close all
x0 = linspace(0,2*pi,100)';
y0 = sin(x0);
for eg=1:4
x = 2*pi*sort(rand(4,1));
y = sin(x);
a = InterpV(x,y);
pVal = HornerV(a,x0);
subplot(2,2,eg)
plot(x0,y0,x0,pVal,'--',x,y,'*')
axis([0 2*pi -2 2])
end
Display cubic interpolants of
sin(x),The abscissas are chosen
randomly
Newton representation
Four point example
))()((
))((),(,1 b a s i s
321
211
xxxxxx
xxxxxx
???
???
))()((
))(()()(
3214
2131213
xxxxxxc
xxxxcxxccxp
????
??????
))()((
))(()()(
))(()(
)(,
3424144
2414314214
2313313213
1221211
xxxxxxc
xxxxcxxccxy
xxxxcxxccy
xxccycy
????
??????
??????
????
))()((
)))(()((
))((
))((
,
342414
2414314214
4
2313
13213
3
12
12
211
xxxxxx
xxxxcxxccy
c
xxxx
xxccy
c
xx
cy
cyc
???
??????
?
??
???
?
?
?
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??????
???
?
4
3
2
1
4
3
2
1
342414241414
231313
12
))()(())(()(1
0))(()(1
00)(1
0001
y
y
y
y
c
c
c
c
xxxxxxxxxxxx
xxxxxx
xx
Matrix-vector notation
General n Case
)()(...)()( 12232 ???????? nn xxxxcxxccxq ?
)()()( 11 xqxxyxp ???
ni
xx
yyx
i
i
i,2,,
1
1 ?
???
?
???
?
?
?
nccc ?,,32
ncccc ?,,,321
function c = InterpNRecur(x,y)
% c = InterpNRecur(x,y)
% The Newton polynomial interpolant.
% x is a column n-vector with distinct components and y is
% a column n-vector,c is a column n-vector with the property that if
%
% p(x) = c(1) + c(2)(x-x(1))+...+ c(n)(x-x(1))...(x-x(n-1))
% then
% p(x(i)) = y(i),i=1:n.
n = length(x);
c = zeros(n,1);
c(1) = y(1);
if n > 1
c(2:n) = InterpNRecur(x(2:n),(y(2:n)-y(1))./(x(2:n)-x(1)));
end
function c = InterpN(x,y)
% c = InterpN(x,y)
% The Newton polynomial interpolant.
% x is a column n-vector with distinct components and y is a column
% n-vector,c is a column n-vector with the property that if
%
% p(x) = c(1) + c(2)(x-x(1))+...+ c(n)(x-x(1))...(x-x(n-1))
% then
% p(x(i)) = y(i),i=1:n.
n = length(x);
for k = 1:n-1
y(k+1:n) = (y(k+1:n)-y(k)),/ (x(k+1:n) - x(k));
end
c = y;
Nonrecursive
implementation
2
1
1
1
1 2
3
2
)1(33)(3 nnnkkn n
k
n
k
????? ??
?
?
?
?
flops
function pVal = HornerN(c,x,z)
% pVal = HornerN(c,x,z)
% Evaluates the Newton interpolant on z where
% c and x are n-vectors and z is an m-vector.
%
% pVal is a vector the same size as z with the property that if
%
% p(x) = c(1) + c(2)(x-x(1))+,.,+ c(n)(x-x(1))...(x-x(n-1))
% then
% pval(i) = p(z(i)),i=1:m.
n = length(c);
pVal = c(n)*ones(size(z));
for k=n-1:-1:1
pVal = (z-x(k)).*pVal + c(k);
end
Nested Multiplication
Efficiency
? Time efficiency
– Vandermode approach,2n3/3 flops(cubic)
– Newton approach,3/2n2 flops(quadratic)
? Memory efficiency
– InterpV requires n-by-n array
– InterpN needs just a few n-vectors
– InterpNRecur requires a couple of n-vectors
with each level of recursion.
Accuracy
b.ηw h e r e a
)x(x)x(x
n!
( ηf
( x )pf ( x )
I,n y xt h e n f o r a
h e x,n t a i n i n g te r v a l I c o a b l e o n a nd i f f e r e n t i
n t i n u o u s l yn t i m e s c o,I f f i s,.,,,xs xi n tt i n c t p oa t t h e d i s
o n f ( x )t h e f u n c t ie r p o l a t e s ( x ) S u p p o s e P
n
( n )
n-
n
n-
??
????
?
?
11
1
1
)
2 T h e o r e m
int
int
If a function has ill-behaved higher derivatives,
then the quality of the polynomial interpolants
may actually decrease as the degree increases.
for n=7:35
xEqual = linspace(-1,1,n)';
yEqual = ones(n,1)./(1+25*xEqual.^2);
cEqual=InterpN(xEqual,yEqual);
pValsEqual = HornerN(cEqual,xEqual,x);
plot(x,y,x,pValsEqual,xEqual,yEqual,'*')
title(sprintf('Equal Spacing (n = %2.0f)',n))
end
the equal-spacing interpolants
of f(x) = 1/(1+25x^2) on [-1,1]'
The Runge phenomenon
),( 11 xfc ?
Divided differences
34
23
12
12
13
13
24
12
12
14
14
4
)()()()()()()()(
xx
xx
xx
xfxf
xx
xfxf
xx
xx
xfxf
xx
xfxf
c
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
23
12
12
13
13
3
)()()()(
xx
xx
xfxf
xx
xfxf
c
?
?
?
?
?
?
?
12
12
2
)()(
xx
xfxfc
?
??
],.,,,[ 1 kk xxfc ?
??
?
?
??
?
?
?? ??
?
??
?
1
11
11 )(],.,,,[)(
k
j
j
n
k
kn xxxxfxp
l y,r e s p e c t i v e},,.,,,{ a n d },.,,,{on
of tsi n t e r p o l a n t h ea r e )( a n d )( S u p p o s e
211 kk-
RL
xxxx
fxpxp
k
RLk
xx
xpxxxpxx
xp
?
???
?
)()()()(
)( 1
kixfxp ii,1),()( ??
)()](,.,,,[
)](,[][)(
111
1211
????
????
kk xxxxxxf
xxxxfxfxp
?
?
)()](,.,,,[
)](,[][)(
2111
1211
?? ???
????
kk
L
xxxxxxf
xxxxfxfxp
?
?
)()](,.,,,[
)](,[][)(
1212
2322
?? ???
????
kk
R
xxxxxxf
xxxxfxfxp
?
?
1
112
1
],...,[],...,[
],...,[
xx
xxfxxf
xxf
k
kk
k ?
?
? ?
]),,,,([ 54321 xxxxxf
]),,,([ 4321 xxxxf ]),,,([
5432 xxxxf
]),,([ 321 xxxf ]),,([ 432 xxxf
]),([ 21 xxf ]),([ 32 xxf
])([ 1xf ])([ 2xf
])([ 2xf ])([ 3xf
]),,([ 432 xxxf ]),,([
543 xxxf
]),([ 32 xxf ]),([ 43 xxf ]),([ 43 xxf ]),([ 54 xxf
])([ 3xf ])([ 4xf
])([ 4xf ])([ 5xf
]),([ 32 xxf ]),([ 43 xxf
])([ 2xf ])([ 3xf
])([ 2xf
])([ 3xf ])([ 4xf
])([ 3xf ])([ 4xf])([ 3xf
Divided differences
])([ 1xf
])([ 2xf
])([ 3xf
])([ 4xf
])([ 5xf
]),([ 21 xxf
]),([ 32 xxf
]),([ 43 xxf
]),([ 54 xxf
]),,([ 321 xxxf
]),,([ 432 xxxf
]),,([ 543 xxxf
]),,,([ 4321 xxxxf
]),,,([ 5432 xxxxf
]),,,,([ 54321 xxxxxf
Efficient computation of divided differences
function c = InterpN2(x,y)
% c = InterpN2(x,y)
% The Newton polynomial interpolant.
% x is a column n-vector with distinct components and y is
% a column n-vector,c is a column n-vector with the property that if
%
% p(x) = c(1) + c(2)(x-x(1))+...+ c(n)(x-x(1))...(x-x(n-1))
% then
% p(x(i)) = y(i),i=1:n.
n = length(x);
for k = 1:n-1
y(k+1:n) = (y(k+1:n)-y(k:n-1)),/ (x(k+1:n) - x(1:n-k));
end
c = y;
hixx i )1(1 ??? Equally spaced
)1(
],...,[],...,[],...,[ 112
1 ?
?? ?
kh
xxfxxfxxf kk
k
?
?
?
????
???
? 1],.,,,[],.,,,[
1)(],.,,,[
112
1
1
kxxfxxf
kxfxxf
kk
k
differences
)!1(
],.,,,[],.,,,[
1
1
1
?
??
? kh
xxfxxf
k
k
k
123452345345455
1234234344
123233
122
1
464332
332
2
fffffffffffffff
ffffffffff
ffffff
fff
f
??????????
??????
???
?
? ?6
?
n?
???
?
?
??
?
?
]2/[ n
n
differences Matlab diff
d=diff(y) d=y(2:n)-y(1:n-1)
d=diff(y,2) d=y(3:n)-2*y(2:n-1)+y(1:n-1)
x=linspace(0,1,6)';
y=x.*x;
a=InterpV(y,x);
yvals=linspace(y(1),y(6));
xvals=HornerV(a,yvals);
plot(yvals,xvals,'r');
hold on
ezplot('sqrt',[0,1]);
plot(y,x,'*')
2)( xxf ? xxg ?)(
xxfg ?))((
Inverse
interpolation
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
x
sqrt(x)
Interpolation in two dimensions
-0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
(a,c) ((1-m)a+? b,c) (b,c)
(a,d) (b,d)((1-m)a+? b,d)
(b,(1-l)c+? d)(a+(1-l)c+? d)
},:),{( dycbxayxR ?????
bdadxd
bcacxc
fff
fff
ll
ll
???
???
)1(
)1(
)~,~()1( yxfffcz xdxc ???? mm
))1(,)1((
))1(())1)((1(
dcbaf
ffffz bdadbcac
mmll
llmllm
?????
???????
function fA = SetUp(fname,a,b,n,c,d,m)
x = linspace(a,b,n);
y = linspace(c,d,m);
fA = zeros(n,m);
for i=1:n
for j=1:m
fA(i,j) = feval(fname,x(i),y(j));
end
end
function u=fun1(x,y)
u=sin(x+y);
a=setup('fun1',0,10,100,0,10,100)
surf(a)
Test code
get the value of
Function fname
function z = LinInterp2D(xc,yc,a,b,c,d,fA)
[n,m]=size(fA);
% xc=a+(i-1+dx)*hx 0<=dx<=1
hx=(b-a)/(n-1);
i=max([1 ceil((xc-a)/hx)]);
dx=(xc (a+(i-1)*hx))/hx;
hy=(d-c)/(m-1);
j=max([1 ceil((yc-c)/hy)]);
dy=(yc (c+(j-1)*hy))/hy;
z=(1-dy)*((1-dx)*fA(i,j)+dx*fA(i+1,j)) + dy*((1-…
dx)*fA(i,j+1)+dx*fA(i+1,j+1));
CEIL Round towards
plus infinity.
Matlab's Polynomial Tools
close all
x = linspace(0,4*pi);
y = exp(-.2*x).*sin(x);
x0 = [.5 4 8 12];
y0 = exp(-.2*x0).*sin(x0);
p = polyfit(x0,y0,length(x0)-1);
pvals = polyval(p,x);
plot(x,y,x,pvals,x0,y0,'o')
axis([0 4*pi -1 1])
set(gca,'XTick',[])
set(gca,'YTick',[])
title('Interpolating {\ite}^{-.2{\itx}}sin({\itx}) with a Cubic
Polynomial')
polyfit
polyval
Interpolating e-.2xsin(x) with a Cubic Polynomial
POLYFIT Fit polynomial to data.
POLYFIT(X,Y,N) finds the coefficients of a polynomial P(X) of
degree N that fits the data,P(X(I))~=Y(I),in a least-squares sense.
[P,S] = POLYFIT(X,Y,N) returns the polynomial coefficients P and
a structure S for use with POLYVAL to obtain error estimates on
predictions,If the errors in the data,Y,are independent normal
with constant variance,POLYVAL will produce error bounds which
contain at least 50% of the predictions.
The structure S contains the Cholesky factor of the Vandermonde
matrix (R),the degrees of freedom (df),and the norm of the
residuals (normr) as fields,
This centering and scaling transformation improves the numerical
properties of both the polynomial and the fitting algorithm,
