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

🔗

Prev Top Next