Newton interpolating polynomial

Section 21.3 Newton interpolating polynomial

Newton polynomial is a third way of constructing an interpolating polynomial of degree

$\lt 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

$\mathbf c$ be the vector of coefficients of interpolating polynomial in some basis; these coefficients are obtained from a linear system

$A\mathbf c=\mathbf 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 $x_j^{i-1}$ and is known as the Vandermonde matrix.)
For the Lagrange basis, the matrix $A$ is the identity matrix, so the coefficients $c_k$ are simply the given values $y_k\text{.}$
For the Newton basis, the matrix $A$ is triangular. This means there is some work to be done to obtain $\mathbf c$ from $\mathbf y$ but it is straightforward.

🔗

The Newton polynomial basis consists of

$\begin{align*} \phi_1(x) \amp = 1\\ \phi_2(x) \amp = x-x_1\\ \phi_3(x) \amp = (x-x_1)(x-x_2)\\ \phi_4(x) \amp = (x-x_1)(x-x_2)(x-x_3) \end{align*}$

🔗

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

$\phi_j(x_k) = 0$ when

$k \lt j\text{.}$

🔗

We are looking for

$c_1, \dots, c_n$ such that the polynomial

$p = c_1\phi_1 + \cdots + c_n\phi_n$ has

$p(x_k)=y_k\text{.}$

At $x=x_1$ only $\phi_1$ is nonzero, so we must have $c_1=y_1$
At $x=x_2$ only $\phi_1,\phi_2$ are nonzero, so $c_1 + c_2 (x-x_1) = y_2\text{.}$ We can solve this for $c_2\text{.}$
At $x=x_3$ only $\phi_1,\phi_2,\phi_3$ are nonzero, so $c_1 + c_2 (x-x_1)+c_3 (x-x_1)(x-x_2) = y_3\text{.}$ We can solve this for $c_3\text{.}$ And so on.

🔗

The process of solving for

$c_k$ can be described algorithmically, in terms of getting the vector

$c$ from the vector

$y\text{.}$ Indeed, we are trying to match two sides of

$\begin{equation} y(x) = c_1 + c_2(x-x_1) + c_3 (x-x_1)(x-x_2) + c_4 (x-x_1)(x-x_2)(x-x_3) + \cdots \label{eq-newton-coeffs}\tag{21.3.1} \end{equation}$

🔗

As noted above,

$c_1 = y_1\text{.}$ Next, subtract

$c_1$ from both sides of (21.3.1)and divide them by

$x-x_1\text{:}$

$\begin{equation} \frac{y(x)-c_1}{x-x_1} = c_2 + c_3 (x-x_2) + c_4 (x-x_2)(x-x_3) +\cdots\label{eq-newton-coeffs-2}\tag{21.3.2} \end{equation}$

🔗

The coefficient

$c_2$ is the value of the left hand side at

$x=x_2\text{.}$ 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

$\begin{equation*} c_3(x-x_2)(x-x_1) + c_2(x-x_1) + c_1 = (c_3 (x-x_{2}) + c_2) (x-x_1) + c_1 \end{equation*}$

🔗

(and similarly for higher degrees), the evaluation of p(t) goes as is: Thus: begin with constant function

$c_n\text{,}$ then in the loop multiply by

$(x-x_k)$ and add

$c_k\text{,}$ where

$k$ goes from

$n-1$ to

$1\text{.}$

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.

🔗

Example 21.3.1. Plot the Newton polynomial through 11 points.

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