Newton polynomial

From Wikipedia, the free encyclopedia

Jump to: navigation, search

In the mathematical field of numerical analysis, a Newton polynomial, named after its inventor Isaac Newton, is the interpolation polynomial for a given set of data points in the Newton form. The Newton polynomial is sometimes called Newton's divided differences interpolation polynomial because the coefficients of the polynomial are calculated using divided differences.

As there is only one interpolation polynomial for a given set of data points it is a bit misleading to call the polynomial Newton interpolation polynomial. The more precise name is interpolation polynomial in the Newton form.

Contents

[edit] Definition

Given a set of k + 1 data points

(x_0, y_0),\ldots,(x_k, y_k)

where no two xj are the same, the interpolation polynomial in the Newton form is a linear combination of Newton basis polynomials

N(x) := \sum_{j=0}^{k} a_{j} n_{j}(x)

with the Newton basis polynomials defined as

n_j(x) := \prod_{i=0}^{j-1} (x - x_i)

and the coefficients defined as

a_j := [y_0,\ldots,y_j]

where

[y_0,\ldots,y_j]

is the notation for divided differences.

Thus the Newton polynomial can be written as

N(x) := [y_0] + [y_0,y_1](x-x_0) + \cdots + [y_0,\ldots,y_k](x-x_0)(x-x_1)\cdots(x-x_{k-1}).

The Newton Polynomial above can be expressed in a simplified form when x0,x1,…,xk are arranged consecutively with equal space. Introducing the notation h=xi + 1-xi for each i=0,1,…,k-1 and x=x0+sh, the difference x-xi can be written as (s-i)h. So the Newton Polynomial above becomes:

N(x):= [y_0] + [y_0,y_1]sh + \cdots + [y_0,\ldots,y_k] s (s-1) \cdots (s-k+1){h}^{k} = \sum_{i=0}^{k}s(s-1) \cdots (s-i+1){h}^{i}[y_0,\ldots,y_i]= \sum_{i=0}^{k}{s \choose i}i!{h}^{i}[y_0,\ldots,y_i]
N(x):= \sum_{i=0}^{k}{s \choose i}i!{h}^{i}[y_0,\ldots,y_i]

is called the Newton Forward Divided Difference Formula.

If the nodes are reordered as xk,xk − 1,…,x0, the Newton Polynomial becomes:

N(x):=[y_k]+[{y}_{k}, {y}_{k-1}](x-{x}_{k})+\cdots+[{y}_{k},\ldots,{y}_{0}](x-{x}_{k})(x-{x}_{k-1})\cdots(x-{x}_{1})

If xk,xk − 1,…,x0 are equally spaced with x=xk+sh and xi=xk-(k-i)h for i=0,1,…,k, then,

N(x):= [{y}_{k}]+ [{y}_{k}, {y}_{k-1}]sh+\cdots+[{y}_{k},\ldots,{y}_{0}]s(s-1)\cdots(s-k+1){h}^{k}={(-1)}^{i}\sum_{i=0}^{k}{-s \choose i}i!{h}^{i}[{y}_{k},\ldots,{y}_{k-i}]
N(x):={(-1)}^{i}\sum_{i=0}^{k}{-s \choose i}i!{h}^{i}[{y}_{k},\ldots,{y}_{k-i}]

is called the Newton Backward Divided Difference Formula.

If nodes x0,x1,…,xk are equally spaced, we would choose to use the Newton Forward Divided Difference Formula to approximate x when x is close to the leading nodes. We would choose to use the Newton Backward Divided Difference Formula to approximate x when x is close to the end of the nodes. In this way, we can make the earliest use of the data points close to x to get a better approximation.

[edit] General case

For the special case of xi = i, there is a closely related set of polynomials, also called the Newton polynomials, that are simply the binomial coefficients for general argument. That is, one also has the Newton polynomials pn(z) given by

p_n(z)={z \choose n}= \frac{z(z-1)\cdots(z-n+1)}{n!}

In this form, the Newton polynomials generate the Newton series. These are in turn a special case of the general difference polynomials which allow the representation of analytic functions through generalized difference equations.

[edit] Main idea

Solving an interpolation problem leads to a problem in linear algebra where we have to solve a system of linear equations. Using a standard monomial basis for our interpolation polynomial we get the very complicated Vandermonde matrix. By choosing another basis, the Newton basis, we get a system of linear equations with a much simpler lower triangular matrix which can be solved faster.

For k + 1 data points we construct the Newton basis as

n_j(x) := \prod_{i=0}^{j-1} (x - x_i) \qquad j=0,\ldots,k.

Using these polynomials as a basis for Πk we have to solve

 
\begin{bmatrix}
      1 &         &        &        & 0  \\
      1 & x_1-x_0 &        &        &    \\
      1 & x_2-x_0 & (x_2-x_0)(x_2-x_1) &        &    \\
 \vdots & \vdots  &        & \ddots &    \\
      1 & x_k-x_0 & \ldots & \ldots & \prod_{j=0}^{k-1}(x_k - x_j)
\end{bmatrix}
\begin{bmatrix}
     a_0 \\
     \vdots \\
     a_{k} 
\end{bmatrix}
=
\begin{bmatrix}
     y_0 \\
     \vdots \\
     y_{k}
\end{bmatrix}

to solve the polynomial interpolation problem.

This system of equations can be solved recursively by solving

 \sum_{i=0}^{j} a_{i} n_{i}(x_j) = y_j \qquad j = 0,\dots,k.

[edit] Taylor polynomial

The limit of the Newton polynomial if all nodes coincide is a Taylor polynomial, because the divided differences become derivatives.


  \lim_{(x_0,\dots,x_n)\to(z,\dots,z)}
     f[x_0] + f[x_0,x_1]\cdot(\xi-x_0) + \dots + f[x_0,\dots,x_n]\cdot(\xi-x_0)\cdot\dots\cdot(\xi-x_n)

   =  f(z) + f'(z)\cdot(\xi-z) + \dots + \frac{f^{(n)}(z)}{n!}\cdot(\xi-z)^n

[edit] Application

As can be seen from the definition of the divided differences new data points can be added to the data set to create a new interpolation polynomial without recalculating the old coefficients. And when a data point changes we usually do not have to recalculate all coefficients. Furthermore if the xi are distributed equidistantly the calculation of the divided differences becomes significantly easier. Therefore the Newton form of the interpolation polynomial is usually preferred over the Lagrange form for practical purposes.

[edit] Example

The divided differences can be written in the form of a table. For example, for a function f is to be interpolated on points x_0, \ldots, x_n. Write

\begin{matrix}
   x_0 & f(x_0) &                                 & \\
       &        & {f(x_1)-f(x_0)\over x_1 - x_0}  & \\
   x_1 & f(x_1) &                                 & {{f(x_2)-f(x_1)\over x_2 - x_1}-{f(x_1)-f(x_0)\over x_1 - x_0} \over x_2 - x_0} \\
       &        & {f(x_2)-f(x_1)\over x_2 - x_1}  & \\
   x_2 & f(x_2) &                                 & \vdots \\
       &        & \vdots                          & \\
\vdots &        &                                 & \vdots \\
       &        & \vdots                          & \\
   x_n & f(x_n) &                                 & \\
\end{matrix}

Then the interpolating polynomial is formed as above using the topmost entries in each column as coefficients.

For example, suppose we are to construct the interpolating polynomial to f(x) = tanx using divided differences, at the points

x0 = − 1.5 x1 = − 0.75 x2 = 0 x3 = 0.75 x4 = 1.5
f(x0) = − 14.1014 f(x1) = − 0.931596 f(x2) = 0 f(x3) = 0.931596 f(x4) = 14.1014

Using six digits of accuracy, we construct the table

\begin{matrix}
-1.5  & -14.1014  &         &          &          &\\
      &           & 17.5597 &          &          &\\
-0.75 & -0.931596 &         & -10.8784 &          &\\
      &           & 1.24213 &          & 4.83484  &  \\
0     & 0         &         & 0        &          & 0\\
      &           & 1.24213 &          & 4.83484  &\\
0.75  & 0.931596  &         & 10.8784  &          &\\
      &           & 17.5597 &          &          &\\
1.5   & 14.1014   &         &          &          &\\
\end{matrix}

Thus, the interpolating polynomial is

\ -14.1014+17.5597(x+1.5)-10.8784(x+1.5)(x+0.75)
\ +4.83484(x+1.5)(x+0.75)(x)+0(x+1.5)(x+0.75)(x)(x-0.75)
\ =-0.00005-1.4775x-0.00001x^2+4.83484x^3

Given more digits of accuracy in the table, the first and third coefficients will be found to be zero.

[edit] Octave Implementation

The following code is a simple Octave implementation of Newton interpolation. The input consists of the data vectors X and Y along with the requested interpolation points, x. The output, y, are the interpolated values at the points x.

function y = newtonInterpolation(X,Y,x)
  % interpolate values y at points x using data vectors X and Y
  n = length(X);
  % build coefficient tableau
  D = diag(Y);
  for m = 1:n
    for i = 1:n-m
      j = i+m;
      D(i,j) = (D(i+1,j)-D(i,j-1))/(X(j)-X(i));
    end
  end
  % return final value
  y = Y(1)*ones(size(x));
  for k = 2:n
    s = D(1,k);
    for i = 1:k-1
      s = s.*(x-X(i)*ones(size(x)));
    end
    y = y+s;
  end
end

[edit] See also

[edit] External links

Personal tools