Processing math: 100%
Skip to main content

Section 21.3 Newton interpolating polynomial

Newton polynomial is a third way of constructing an interpolating polynomial of degree <n through n given points. It is still the same polynomial, since there is only one such polynomial. But the method is different because, yet again, we use a different basis. The Newton basis is somewhere between monomial basis and Lagrange basis, in the following sense. Let c be the vector of coefficients of interpolating polynomial in some basis; these coefficients are obtained from a linear system Ac=y where the square matrix A depends on the basis.

  • For the monomial basis, the matrix A is dense and often ill-conditioned (it has entries xi1j and is known as the Vandermonde matrix.)
  • For the Lagrange basis, the matrix A is the identity matrix, so the coefficients ck are simply the given values yk.
  • For the Newton basis, the matrix A is triangular. This means there is some work to be done to obtain c from y but it is straightforward.

The Newton polynomial basis consists of

ϕ1(x)=1ϕ2(x)=xx1ϕ3(x)=(xx1)(xx2)ϕ4(x)=(xx1)(xx2)(xx3)

and so on. An advantage of Newton polynomial is that it is easier to evaluate than Lagrange polynomial. Also, we will see that adding a new data point is very easy to handle because the computation processes one step at a time. A key fact is that ϕj(xk)=0 when k<j.

We are looking for c1,,cn such that the polynomial p=c1ϕ1++cnϕn has p(xk)=yk.

  • At x=x1 only ϕ1 is nonzero, so we must have c1=y1
  • At x=x2 only ϕ1,ϕ2 are nonzero, so c1+c2(xx1)=y2. We can solve this for c2.
  • At x=x3 only ϕ1,ϕ2,ϕ3 are nonzero, so c1+c2(xx1)+c3(xx1)(xx2)=y3. We can solve this for c3. And so on.

The process of solving for ck can be described algorithmically, in terms of getting the vector c from the vector y. Indeed, we are trying to match two sides of

y(x)=c1+c2(xx1)+c3(xx1)(xx2)+c4(xx1)(xx2)(xx3)+

As noted above, c1=y1. Next, subtract c1 from both sides of (21.3.1)and divide them by xx1:

y(x)c1xx1=c2+c3(xx2)+c4(xx2)(xx3)+

The coefficient c2 is the value of the left hand side at x=x2. In Matlab terms, it's the second entry of the array of y-values after we applied the subtract-and-divide process to it.

This repeats: we start with c=y and repeat the subtract-and-divide step to turn its term into the coefficients of the Newton polynomial.

c = y;
for k=2:n
    c(k:n) = (c(k:n)-c(k-1))./(x(k:n)-x(k-1))
end

At the end of the loop, the vector c contains the coefficients of Newton polynomial.

Evaluation of Newton's polynomial is done in a similar way, with repeated multiply-and-add process, using one coefficient at a time. Indeed, since

c3(xx2)(xx1)+c2(xx1)+c1=(c3(xx2)+c2)(xx1)+c1

(and similarly for higher degrees), the evaluation of p(t) goes as is: Thus: begin with constant function cn, then in the loop multiply by (xxk) and add ck, where k goes from n1 to 1.

p = c(n)*ones(size(t));
for k=n-1:-1:1
    p = p.*(t-x(k)) + c(k);
end

While being simpler than barycentric evaluation (21.2.5) of the Lagrange polynomial, the evaluation of Newton polynomial does not fit into an anonymous function because it uses a loop.

Plot the Newton polynomial for the points \((k, \arctan k)\text{,}\) \(k=-5, \dots, 5\text{.}\) Then repeat with 5 replaced by 8.

Answer
x = -5:5;
n = numel(x);
y = atan(x); 

c = y;
for k=2:n
    c(k:n) = (c(k:n)-c(k-1))./(x(k:n)-x(k-1));
end

t = linspace(min(x)-0.01, max(x)+0.01, 1000);
p = c(n)*ones(size(t));
for k=n-1:-1:1
    p = p.*(t-x(k)) + c(k);
end
plot(t, p, x, y, 'r*')