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 xi−1j 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.
ϕ1(x)=1ϕ2(x)=x−x1ϕ3(x)=(x−x1)(x−x2)ϕ4(x)=(x−x1)(x−x2)(x−x3)
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(x−x1)=y2. We can solve this for c2.
- At x=x3 only ϕ1,ϕ2,ϕ3 are nonzero, so c1+c2(x−x1)+c3(x−x1)(x−x2)=y3. We can solve this for c3. And so on.
y(x)=c1+c2(x−x1)+c3(x−x1)(x−x2)+c4(x−x1)(x−x2)(x−x3)+⋯
As noted above, c1=y1. Next, subtract c1 from both sides of (21.3.1)and divide them by x−x1:
y(x)−c1x−x1=c2+c3(x−x2)+c4(x−x2)(x−x3)+⋯
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)) endAt 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(x−x2)(x−x1)+c2(x−x1)+c1=(c3(x−x2)+c2)(x−x1)+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 (x−xk) and add ck, where k goes from n−1 to 1.
p = c(n)*ones(size(t)); for k=n-1:-1:1 p = p.*(t-x(k)) + c(k); endWhile 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*')