## The Newton Polynomial Interpolation

The Lagrange interpolation relies on the interpolation points , all of which need to be available to calculate each of the basis polynomials . If additional points are to be used when they become available, all basis polynomials need to be recalculated.

In comparison, in the Newton interpolation, when more data points are to be used, additional basis polynomials and the corresponding coefficients can be calculated, while all existing basis polynomials and their coefficients remain unchanged. Due to the additional terms, the degree of interpolation polynomial is higher and the approximation error may be reduced (e.g., when interpolating higher order polynomials).

Specifically, the basis polynomials of the Newton interpolation are calculated as below:

 (25)

and the Newton interpolating polynomial is constructed:
 (26)

When when the next data point is available, not used in any of the basis polynomials, but it is used for calculating the last coefficient , as shown below. For this nth degree polynomial to pass all points , it needs to satisfy the following equations:

 (27)

which can also be expressed in matrix form:

 (28)

The coefficients can be obtained by solving these equations in the triangular equation system progressively from top down:
 (29)

In general, we have

 (30)

which is the expanded form of the kth divided differences of the first points. Now the Newton polynomial interpolation can be written as

 (31)

Due to the uniqueness of the polynomial interpolation, this Newton interpolation polynomial is the same as that of the Lagrange and the power function interpolations: . They are the same nth degree polynomial but expressed in terms of different basis polynomials weighted by different coefficients.

We can now consider some important facts all related to the Newton polynomial interpolation.

• When an additional node point is avialable and to be used, all previous basis polynomials and their corresponding coefficients remain unchanged, we only need to obtain a new basis polynomial of degree :

 (32)

together with its coefficient, the (n+1)th divided difference . The new interpolation polynomial of degree can then be obtained by including an extra term in the summation above:

 (33)

As passes through the new point , we have

 (34)

But is just an arbitrary point, we can replace it by , and get

 (35)

Comparing this with the error function:

 (36)

we see that the error term can also be written as

 (37)

and we also get:

 (38)

• Specially, if all of the node points approach to a single position, , at the limit they become the same point repeated times, then we have

 (39)

and the Newton interpolation based on points becomes

 (40)

which is the first terms of the Taylor series expansion of the function with the truncation error:

 (41)

where is some point between and . We see that the Taylor series is actually a special case of the Newton interpolation at the limit when all node points approach to the same position .

• If all points are equally spaced, i.e.,

 (42)

then Newton's divided difference interpolation can take a simpler form. For any point , we let so that it can be represented as , and , now the Newton polynomial can be written as
 (43)

where

 (44)

Example: Approximate function by a polynomial of degree , based on the following points:

Based on , we can find all other divided differences recursively in tabular form as shown below. In general, can be found based on its left neighbor and top-left neighbor :

The coefficients are the four divided differences along the diagonal: , , , and . Alternatively, they can also be represented in the expanded form:

 (45)

Now the Newton interpolating polynomial can be obtained as:

The polynomial interpolations generated by the power series method, the Lagrange and Newton interpolations are exactly the same, , confirming the uniqueness of the polynomial interpolation, as plotted in the top panel below, together with the original function . We see that they indeed pass through all node points at , , and . Also, the weighted basis polynomials of each of the three methods are plotted in the subsequent panels, including power series in the second panel, the Lagrange basis polynomials in the third panel, and the Newton basis polynomials in the bottom panel. In each case, the weighted sum of these basis polynomials is the interpolating polynomial that approximates the given function.

The Matlab code that implements the Newton polynomial method is listed below. The coefficients can be generated in either the expanded form or the tabular form by recursion.

function [v N]=NI(u,x,y)  % Newton's Interpolation
% vectors x and y contain n+1 points and the corresponding function values
% vector u contains all discrete samples of the continuous argument of f(x)
n=length(x);          % number of interpolating points
k=length(u);          % number of discrete sample points
v=zeros(1,k);         % Newton interpolation
N=ones(n,k);          % all n Newton's polynomials (each of m elements)
N(1,:)=y(1);          % first Newton's polynomial
v=v+N(1,:);
for i=2:n             % generate remaining Newton's polynomials
for j=1:i-1
N(i,:)=N(i,:).*(u-x(j));
end
c=DividedDifference(x,y,i)  % get the ith coefficient c_i
v=v+c*N(i,:);     % weighted sum of all Newton's polynomials
end
end

function dd=DividedDifference(x,y,i) % generate f[x_0,...,x_i] in expanded form
dd=0;
for k=1:i             % loop for summation
de=1;
for l=1:i         % loop for product
if k~=l
de=de*(x(k)-x(l));
end
end
dd=dd+y(k)/de;    % ith coefficient c_i
end
end

function dd=DividedDifferenceMatrix(x,y) % generate divided difference matrix
n=length(x);                         % the coefficients are along diagonal
dd=zeros(n);                         % matrix of divided differences
dd(:,1)=y;
for i=1:n
fprintf('%6.3f\t',dd(i,1))
for j=2:i
dd(i,j)=(dd(i,j-1)-dd(i-1,j-1))/(x(i)-x(i-j+1));
fprintf('%6.3f\t',dd(i,j));
end
fprintf('\n');
end
end