Chapter 2 Polynomial Interpolation
The purpose of polynomial interpolation is to fit tada using polynomial model. If (x1, y1),…, (xn,yn) are given, then to find a polynomial P(x) of degree n-1 such that P(xi)=yi, 1<=I<=n is called polynomial interpolation, while to find a polynomial Q(x) of degree <n is called least squares fitting.
MATLAB uses the same built-in function "polyfit' for both polynomial interpolation and polynomial least squares fitting.
Function a=polyfit(x,y,n);
where (x,y) is the data need to interpolate and n is the degree of the interpolating polynomial. The output "a" is the coefficient vector of the polynomial.
In MATLAB, the polynomial is given as a vector, which represents the coefficients of the polynomial. To evaluate a polunomial at xvals, the following function should be called.
Function pval=polyval(a,xvals);
where "a" is the coefficient of the polynomial and "xvals " is the vector of the x values. The output "pval" is the polynomial value vector.
For example, the given data are (1,3), (2,4), (3,-2), (4,-1), (5, 3). To find its polynomial interpolation:
y=[3, 4, -2,-1,3];
p=polyfit(x,y,4)
-0.7500 9.8333 -43.7500 74.6667 -37.0000
To evaluate its value at a=[0, 1, 3, 10, 20], do the following.
0 1 3 10 20
1.0e+004 *
-0.0037 0.0003 -0.0002 -0.1332 -5.7377
In this chapter, we discuss the method to find the interpolating polynomail and the different representation of the same polynomial so that the evaluation will be more effective.
2.1 The Vandermonde Approach.
2.1.1 A Four-point Interpolation Problem
Example: To find a polynomial, which interpolates the data (-2,10), (-1,4), (1,6), and (2,3).
Assume that the polynomial is
P(x)=a1+a2*x+a3*x^2+a4*x^3,
then we have
a1+a2*(-2)+a3*(-2)^2+a4*(-2)^3=10
a1+a2*(-1)+a3*(-1)^2+a4*(-1)^3=4
a1+a2*1+a3*1^2+a4*1^3=6
a1+a2*x+a3*2^2+a4*2^3=3
Set C=[1,-2,4,-8; 1,-1,1,-1;1, 1, 1, 1;1,2,4,8] ;
1 -2 4 -8
1 -1 1 -1
1 1 1 1
1 2 4 8
10
4
6
3
Then a=C\y
4.5000
1.9167
0.5000
-0.9167
4.5000
1.9167
0.5000
-0.9167
We can verify this is the right polynomial.
10 4 6 3
2.1.2 The General Case
In general, if we want to interpolate the data (x1, y1),…, (xn,yn) using a polynomial P(x) of degree n-1.
Represent the polynomial P(x) by P=[an, a(n-1), …, a1].
Set a=[a1,a2,…,an] (=fliplr(P))
The coefficents matrix of the linear system is the Vendermonde matrix.
V=[x1, …,x1^(n-1); x2,…,x2^(n-1); xn,…,xn^(n-1)];
Then the vector "a" is the solution of the linear system.
a=V\y.
2.1.3 Setting Up and Solving the System
The whole program for creating p(x) is the following.
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;
Note there is a buildin function for creating Vandermode Matrix. But it uses thereverse order.
Example:
vander(x)
1 1 1 1
8 4 2 1
27 9 3 1
64 16 4 1
Using it to set the equation, the result vector "a" can be directly used for evaluate the polynomial by function "polyval".
2.1.4 Nested Multiplication
In this section, we will discuss how to aasign a better way to evaluate polynomial.
For example, we want to evaluate polynomial
p(x) = a(1) + a(2)x + ... a(n)x^(n-1)
at a vector "z".
The direct way to do so is set
powerz=[1;z;z.^2;…;z.^(n-1)] (need use a loop to create it.) Then get the values using
pval=a*powerz;
However this way often causes overflow. It also costs more computing time.
The nested multiplication is based on the following,which is called Horner's rule.
p(x) = (...( (a(n)*x+a(n-1))*x+a(n-2))*x+…+a(2))*x+a(1).
The following program is based on this rule.
function pVal = HornerV(a,z)
% pVal = HornerV(a,z)
% evaluates the Vandermonde interpolant on z where
% a is an n-vector and z is an m-vector.
%
% pVal is a vector the same size as z with the property that if
%
% p(x) = a(1) + .. +a(n)x^(n-1)
% then
% pVal(i) = p(z(i)) , i=1:m.
n = length(a);
m = length(z);
pVal = a(n)*ones(size(z));
for k=n-1:-1:1
pVal = z.*pVal + a(k);
end
The following script file demos the Vandermonde method.
% 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
2.2 The Newton Representation
Because of the unstability of Vandermonde matrix, the Vandermonde method has a certain limitation. A better way to represent interpolating polynomial is the Newton Method.
2.2.1 A Four-Point Example.
We start with a four-point example.
Let (x1,y1),(x2,y2),(x3,y3), and (x4,y4) are given. To find a polynomial p(x)=c1+c2(x-x1)+c3(x-x1)(x-x2)+c4(x-x1)(x-x2)(x-x3),
which interpolates the data.
We can get the following linear sysytem:
c1=y1
c1+c2(x2-x1)=y2
c1+c2(x3-x1)+c3(x3-x1)(x3-x2)=y3
c1+c2(x4-x1)+c3(x4-x1)(x4-x2)+c4(x4-x1)(x4-x2)(x4-x3)=y4.
We use the elimination method to solve it.
_________________________________________________
c1 =y1
c2 =(y2-y1)/(x2-x1)=y21
c2 +c3(x3-x2) =(y3-y1)/(x3-x1)=y31
c2 +c3(x4-x2)+c4(x4-x2)(x4-x3) =(y4-y1)/(x4-x1)=y41
_________________________________________________
c1 =y1
c2 =y21
c3 =(y31-y21)/(x3-x2)=y321
c3 +c4(4-x3) =(y41-y21)/(x4-x2)=y421
_________________________________________________
c1 =y1
c2 =y21
c3 =y321
c4 =(y421-y321)/(x4-x3)=y4321
__________________________________________________
In the process, the vectors we need are
y1 y2 y3 y4
y21 y31 y41
y321 y421
y4321
2.2.2 The General Case
In the general case, we can develop the idea in the previous sub section.
The vectors used to find "c" now is
y1 y2 y3 y4 ……….. yn
y21 y31 y41 ……….. yn1
y321 y421 ………. yn21
y4321 ……. yn321
…………………………………………….
yn…4321
The program for obtain Newton's representation is the following.
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;
When in a function each step is repeatable, we can use the recurrent function. For example, we can use the following recurrent function to represent Newton's polynomial.
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
In the function InterpNRecur(x,y), when you meet a calling for the function, the program will stop at the point and then return to the beginning of the program.
2.2.3 Nested Multiplication
The nested multiplication for evaluating Newton's polynomial can be performed as
p(x) = (…((c(n)(x-x(n-1))+c(n-1))(x-x(n-2))+...+c(2))(x-x(1))+c(1).
The program is the following.
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
2.3 Properties
We can interpolate data using Vandermonde approach or Newton's representation. Which one is better? Speed and accuracy are two main concerns for the cpmparison. In Subsection 2.3.1, we discuss the spped, in 2.3.2 we study the accuracy.
2.3.1 Efficiency
We first compare the number of operations used in the programs. The speeds of different arithmetic operations are different in a computer. We standard them by a uniform measure called "flop". Then we use flop counts the speed. The build-in funvtion "flops" can count the number. The following script compares the flops for Vandermonde and Newton approaches.
% Script File: InterpEff
% Compares the Vandermonde and Newton Approaches
disp('Flop Counts:')
disp(' ')
disp(' n InterpV InterpN InterpNRecur')
disp('--------------------------------------')
for n = [4 8 16]
x = linspace(0,1,n)'; y = sin(2*pi*x);
flops(0); a = InterpV(x,y); f1 = flops;
flops(0); c = InterpN(x,y); f2 = flops;
flops(0); c = InterpNRecur(x,y); f3 = flops;
disp(sprintf('%3.0f %7.0f %7.0f %7.0f',n,f1,f2,f3));
end
n InterpV InterpN InterpNRecur
--------------------------------------
4 179 28 18
8 847 106 84
16 4671 406 360
From the result, we can see that the most effective algorithm is InterNRecur. The number of flops in InterpV is O(n^3) while the number of flops in InterpN is O(n^2).
2.3.2 Accuracy.
InterpV and InterpN produce the same polynomial in different forms. Then the accuracy for them is the same. The error can be estimated by the following formula. Let p(x) be the polynomial interpolating the data (x1,y1), (x2,y2), …, (xn,yn), where f(xn)=yn for an existent function f(x). Let [a,b] be the interval containing x1, x2, …, xn, and x. Let Mn be the maximum of |f(n)(x)|. Then the error can be represented by
|f(x)-p(x)|<=Mn/n!*max|(x-x1)(x-x2)…(x-xn)|.
If the interpolation is on the equal spaced points
xi=a+(b-a)(i-1)/(n-1),
then the error bound is
|f(x)-p(x)|<=Mn/(4n)((b-a)/(n-1))^n.
Thus, if a function has ill-behaved higher derivatives, then the quality of the polynomial interpolants may actually decrease as the degree increases.
As an example we show the phemomenon for the function f(x)=1/(1+25x^2) across the interval [-1,1].
% Script File: RungeEg
% For n = 7:2:15, the equal-spacing interpolants of f(x) = 1/(1+25x^2) % % on [-1,1]
% are plotted.
% create the data for plot the function f(x) = 1/(1+25x^2).
x = linspace(-1,1,100)';
y = ones(100,1)./(1 + 25*x.^2);
% Make the polynomial interplants with degree 6,8,10,12, and 14.
% The Newton's representation is used for equal spaces data.
for n=7:2:15
figure
xEqual = linspace(-1,1,n)';
yEqual = ones(n,1)./(1+25*xEqual.^2);
cEqual=InterpN(xEqual,yEqual);
% Evaluate polynomial at x to create the data for plotting polynomial.
pValsEqual = HornerN(cEqual,xEqual,x);
% plot function, polynomial interplant, and the data in a figure.
plot(x,y,x,pValsEqual,xEqual,yEqual,'*')
title(sprintf('Equal Spacing (n = %2.0f)',n))
pause(3)
end
2.4 Special Topics
2.4.1 Divided Differences
For given function f(x), its divided differences at x=[x1,x2, ..,xn] are defined by the following.
f[xi]=f(xi), i=1…n. (The zero order)
f[x(i-1), x(i)]=(f[xi]-f[x(I-1)])/(xi-x(i-1)) (The first order)
f[x(i-1), xi, x(i+ 1)]=(f[x(i+ 1),xi]) - (f[x(i), x(i- 1)])/(x(i+ 1)-x(i-1)).
Let y=[f(x1);f(x2), …, f(xn)].
The following function can create the table of divided differences for a given function and x.
function T=fdivd(fname,x)
% T is a cell array containing the divided differences of a given function up to order n-1 at x, where n is the length of x.
n=length(x);
T{1}=fval(fname,x);
for j=2:n
T{j}=diff(T{j-1})./(x(j:n)-x(1:n-j+1));
end
The divided differences for data can be calculated by the following.
function T=divd(x,y)
% T is a cell array containing the divided differences up to order n-1 of data (x,y), where n is the length of % x and length(y)=length(x).
function T=divd(x,y)
n=length(x);
T{1}=y;
for j=2:n
T{j}=diff(T{j-1})./(x(j:n)-x(1:n-j+1));
end
Note that if p(x) is a polynomial of degree n. Then its n-th order divided difference is a constant, and therefore, its (n+1)-th divided difference is zero.
Thus the Newton's representation of interpolating polynomial for (x1,y1), (x2,y2) …, (xn,yn) is
p(x)=f[x1]+f[x1,x2](x-x1)+…+f[x1,x2,…xn](x-x1)(x-x2)…(x-x(n-1)).
The following program uses this formula to compute the coefficients of Newton representation.
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;
2.4.2 Inverse Interpolation
If a function f(x) has its inverse g(x) on an interval [a,b]. Then the y-intercept of g(x) is the x-intercept of f(x). Recall that the x-intercept of a function is the zero of this function. Thus, g(0) is the solution of the equation f(x)=0. This way is espacially useful for finding the zero of a function in an interval if the function has the inverse on that interval.
To get the interpolation of the inverse function, you only nedd to interchage the x and y vectors of data.
For example the following program give the interpolation of sqrt(x) using inverse interpolation.
%create the data from x^2 for inverse interpolation for sqrt(x).
x=linspace(0,1,6)';
y=x.*x;
% make the inverse interpolation.
a=InterpV(y,x);
% plot the interpolating polynomial.
yvals=linspace(y(1),y(6));
xvals=HornerV(a,yvals);
plot(yvals,xvals)
