Newton interpolation

1The method¶

Consider $N+1$ distinct points $x_0,x_1,\ldots,x_N$ at which data about a function is given. The Newton form of interpolation makes use of the basis functions

so that the degree $N$ polynomial has the form

The two practical questions are:

1. How to determine the coefficients $\{ c_i \}$ ?

2. How to evaluate the resulting polynomial ?

The interpolation conditions are given by

We have a matrix equation with a lower triangular matrix which can be solved recursively starting from the first equation. There is no need to invert the matrix. The non-singularity of the matrix is clear since the determinant is the product of the diagonal terms and this is non-zero since the points $\{ x_i \}$ are distinct. The coefficients $\{ c_i \}$ can be computed in $O(N^2)$ operations but there is danger of overflow/ underflow errors if we do not do this carefully. The triangular structure means that if we add a new data point $(x_{N+1},f_{N+1})$, then the coefficients $c_0, c_1, \cdots c_N$ do not change; we have to compute the additional coefficient $c_{N+1}$.

2Divided differences¶

Define zeroth divided difference

and the first divided difference

The $k$’th divided difference is defined recursively by

We claim that

We can compute the divided differences by constructing the following table

The diagonal terms are the coefficients $\{c_i\}$ of Newton interpolation. To compute the the diagonal terms, it is not necessary to store all the other terms in memory.

for(i=0; i<=N; ++i) c[i] = f[i];
for(j=1; j<=N; ++j)
   for(i=N; i>=j; --i)
      c[i] = (c[i] - c[i-1])/(x[i] - x[i-j]);

This requires $N^2$ additions and $\frac{1}{2}N^2$ divisions.

3Evaluation of polynomials¶

Suppose given $x$ we have to compute

A direct term-by-term evaluation as in the above formula can suffer from overflow/ underflow error and the cost is also high. To see how to evaluate more efficiently, let us consider some simpler examples for quadratic

and cubic polynomials,

This suggests the following recursive algorithm called Horner’s method

which gives $p(x)$ at the end. This approach can be used for the Newton form also.

This requires $2N$ additions and $N$ multiplications so the cost is $O(N)$. The direct term-by-term computation

double x[N], c[N];
double p = 0.0;
for(int i=0; i<=N; ++i)
{
   double tmp = 1.0;
   for(int j=0; j<i; ++j) tmp *= (x - x[j]);
   p += c[i] * tmp;
}

would require $N$ additions and

multiplications for a total cost of $O(N^2)$ for each evaluation of $p(x)$. This can be improved by a simple modification

double x[N], c[N];
double tmp = 1.0, p = 0.0;
for(int i=0; i<=N; ++i)
{
   p += c[i] * tmp;
   tmp *= (x-x[i]);
}

4Proof of $c_k = f[x_0,x_1,\ldots,x_k]$¶

The polynomial interpolating the data at $\{ x_0,x_1,\ldots,x_k \}$ is

Note that $c_k$ depends on $x_0,x_1,\ldots,x_k$ so we indicate this dependence symbolically as

We will show that this symbol satisfies the recursion relation for divided differences. Note that $c_k$ is the coefficient of $x^k$ in the polynomial $p_k(x)$.

Let us also write the Lagrange polynomial interpolating the same data as above. Define

so that

Then the Lagrange form can be written as

The coefficient of $x^k$ in this form is

since the interpolating polynomial is unique.

Now, using above result, we can write a formula for $f[x_0,x_1,\ldots,x_{k-1}]$ but we have to knock off the factor $(x_j-x_k)$ from $\psi_k'(x_j)$, so that

Similarly

Hence

which is the recursion relation for divided differences.

5Interpolation error¶

The Newton form of interpolation to the data $(x_i,f_i)$, $i=0,1,\ldots,N$ is

Let $t \in \re$ be distinct from the nodes $x_0,x_1,\ldots,x_N$ but otherwise arbitrary. The polynomial interpolating $f(x)$ at $x_0,x_1,\ldots,x_N,t$ is

But $p_{N+1}(t) = f(t)$ so that

Since $t$ was arbitrary, we could call it $x$. Then the error is

Comparing this with the formula we have derived earlier, we conclude that

6Another proof of independence of order¶

We will show that $f[x_0,x_1,\ldots,x_N]$ is independent of the order of the arguments. Let $p_{N-1}(x)$ interpolate the data $\{ x_0, x_1, \ldots, x_{N-1}\}$. In Newton form

Consider a random permutation

and let $Q_{N-1}(x)$ interpolate this data

But these two polynomials are same $Q_{N-1} = P_{N-1}$ since they interpolate same data, and in particular the error at $x=x_N$ is same

Since the polynomial factors are same, we get

We can repeat this proof by leaving out $x_{N-1}, x_{N-2}, \ldots, x_0$ and then we obtain the complete proof.

7Forward difference operator¶

For $h > 0$, define the forward difference

We will use uniformly spaced data at

Then

or, in compact form

For $r \ge 0$, define

$\Delta^{r+1} f(x)$ is called the $r$-th order forward difference of $f$ at $x$.

8Newton form on uniform data¶

For uniformly spaced data, we can avoid the use of divided differences and write in terms of forward differences, which is more efficient since we avoid division, which is an expensive operation. Define

Newton interpolation formula has factors of the form $(x-x_0)\ldots(x-x_k)$. Since

hence

Using previous lemma, we write divided difference in terms of forward differences

Define the binomial coefficients

Then

This is the Newton forward difference form of the interpolating polynomial. A schematic of the computation is shown in Table 1.

9Errors in data and forward differences¶

We can use a forward difference table to detect noise in physical data, as long as the noise is larger relative to the experimental error.

10Some properties of divided differences¶

1. Recall that the divided difference can be written as

   $$f[x_0,\ldots,x_n] = \frac{f^{(n)}(\xi)}{n!}$$

   (63)

   Prove this result by induction. The Hermite-Gennochi formula relates the sames quantities in terms of an integral.

2. If $f^{(n)}$ is continuous in the interval $[\min_j x_j, \max_j x_j]$, then $f[x_0,\ldots,x_n]$ is a continuous function of $x_0, \ldots,x_n$.

3. From the Hermite-Gennochi formula, we see that the right hand side makes sense even if all the $x_j$ are not distinct. If we take all arguments to be equal to $x$, then

   $$\begin{aligned} f[\underbrace{x,x,\ldots,x}_{n+1 \textrm{ arguments}}] =& \int_{\tau_n} f^{(n)}(x) \ud t_1 \ldots \ud t_n \\ =& f^n(x) \textrm{ volume}(\tau_n) \\ =& \frac{f^{(n)}(x)}{n!} \end{aligned}$$

   (64)

4. Consider $f[x_0,\ldots,x_n,x]$ as a function of $x$ and differentiate this. By definition of the derivative

   $$\begin{aligned} \dd{}{x}f[x_0,\ldots,x_n,x] =& \lim_{h \to 0} \frac{f[x_0,\ldots,x_n,x+h] - f[x_0,\ldots,x_n,x]}{h} \\ =& \lim_{h \to 0} \frac{f[x_0,\ldots,x_n,x+h] - f[x,x_0,\ldots,x_n]}{h} \\ =& \lim_{h \to 0} f[x,x_0,\ldots,x_n,x+h] \\ =& f[x,x_0,\ldots,x_n,x] \\ =& f[x_0,\ldots,x_n,x,x] \end{aligned}$$

   (65)

   The quantity on the right has $n+3$ arguments and it is well defined if $f \in \cts^{n+2}$.