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) |

(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) |

(28) |

(29) |

In general, we have

(30) |

(31) |

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) (33) (34) (35) (36) (37) (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) (40) (41) - If all points
are
equally spaced, i.e.,
(42)

(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) |

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