Chapter 3 — Interpolation and Polynomial Approximation

Given data — values of some function at a handful of points — we want a simple formula that reproduces them. Polynomials are the natural first choice: they are cheap to evaluate, trivial to differentiate and integrate, and by the Weierstrass approximation theorem they can come arbitrarily close to any continuous function. This chapter builds the interpolating polynomial three different ways, then abandons a single high-degree polynomial in favour of piecewise cubics, which behave far better in practice.

3.1 Interpolation and the Lagrange Polynomial

Building the polynomial from a basis

The Lagrange form is worth understanding as a change of basis rather than a formula to memorise. We look for \(P\) in the form \[ P(x) = f(x_0)L_{n,0}(x) + f(x_1)L_{n,1}(x) + \cdots + f(x_n)L_{n,n}(x), \] and ask what the functions \(L_{n,k}\) must do. If each one satisfies \[ L_{n,k}(x_j) = \begin{cases} 1 & \text{if } k = j,\\ 0 & \text{if } k \ne j, \end{cases} \] then evaluating \(P\) at \(x_j\) kills every term but the \(j\)-th and leaves exactly \(f(x_j)\). So the interpolation conditions are satisfied by construction, with no linear system to solve.

Such a function is easy to write down. To vanish at every node except \(x_k\), take a product with a factor \((x - x_i)\) for each \(i \ne k\); to equal \(1\) at \(x_k\), divide by the value of that product at \(x_k\): \[ L_{n,k}(x) = \prod_{i \ne k} \frac{x - x_i}{x_k - x_i}. \]

It is worth seeing the first few cases.

For \(n = 0\) there is one node and nothing to interpolate but a constant, so the basis is the single function \(L_{0,0}(x) = 1\).

For \(n = 1\) the two basis functions are straight lines, each rising to \(1\) at its own node and crossing zero at the other.

For \(n = 2\) the basis functions are quadratics, each vanishing at the two nodes that are not its own.

The pattern continues: for \(n\) nodes the basis functions are degree-\(n\) polynomials with \(n\) prescribed roots.

The error term

The error formula should look familiar: it is Taylor’s theorem with the single factor \((x-x_0)^{n+1}\) replaced by the product \(\prod_{i=0}^n (x - x_i)\), one factor for each node. In particular the error vanishes at every node, as it must, and grows between and beyond them.

3.3 Divided Differences

The Lagrange form has an awkward property: adding one more node means recomputing every basis function from scratch. The Newton form fixes this by writing the polynomial in nested products, \[ \begin{aligned} P_n(x) = a_0 &+ a_1(x-x_0) + a_2(x-x_0)(x-x_1) \\ &+ \cdots + a_n(x-x_0)(x-x_1)\cdots(x-x_{n-1}), \end{aligned} \] so that each new node appends a term and leaves the earlier ones untouched.

The coefficients now come out one at a time. Evaluating at \(x_0\), every term but the first vanishes: \[ P_n(x_0) = a_0 = f(x_0) \equiv f[x_0]. \] Evaluating at \(x_1\), only the first two survive: \[ P_n(x_1) = a_0 + a_1(x_1 - x_0) = f(x_1) \qquad\Longrightarrow\qquad a_1 = \frac{f(x_1) - f(x_0)}{x_1 - x_0} \equiv f[x_0, x_1]. \]

That pattern is the definition of the divided differences: \[ \begin{aligned} f[x_i] &= f(x_i), \\[2pt] f[x_i, x_{i+1}] &= \frac{f[x_{i+1}] - f[x_i]}{x_{i+1} - x_i}, \\[2pt] f[x_i, x_{i+1}, x_{i+2}] &= \frac{f[x_{i+1}, x_{i+2}] - f[x_i, x_{i+1}]}{x_{i+2} - x_i}, \end{aligned} \] and in general \[ f[x_i, x_{i+1}, \ldots, x_{i+k}] = \frac{f[x_{i+1}, \ldots, x_{i+k}] - f[x_i, \ldots, x_{i+k-1}]}{x_{i+k} - x_i}. \]

With these, the interpolating polynomial is Newton’s divided-difference formula: \[ \begin{aligned} P_n(x) = f[x_0] &+ f[x_0,x_1](x-x_0) + f[x_0,x_1,x_2](x-x_0)(x-x_1) \\ &+ \cdots + f[x_0,\ldots,x_n](x-x_0)\cdots(x-x_{n-1}). \end{aligned} \]

function F = divideddifference(x, f)
% Compute an interpolation table using divided differences.

n = length(x) - 1;
F = zeros(n + 1, n + 1);
F(:, 1) = f(:);
for i = 1:n
    for j = 1:i
        F(i + 1, j + 1) = (F(i + 1, j) - F(i, j)) / ...
            (x(i + 1) - x(i - j + 1));
    end
end

Divided differences are derivatives

The notation is suggestive, and the suggestion is correct. The first divided difference is a difference quotient, so the mean value theorem gives \[ f[x_0,x_1] = \frac{f(x_1) - f(x_0)}{x_1 - x_0} = f'(\xi) \qquad\text{for some } \xi \in [x_0, x_1]. \]

This generalises to any order, and the general statement is the reason divided differences appear in every error formula in the chapter: if \(f \in C^n[a,b]\) and \(x_0, \ldots, x_n\) are distinct points of \([a,b]\), then there exists \(\xi \in [a,b]\) with \[ f[x_0, \ldots, x_n] = \frac{f^{(n)}(\xi)}{n!}. \]

Equally spaced nodes

When the nodes are equally spaced the formula simplifies considerably. Write \(h = x_{i+1} - x_i\) and introduce a new variable \(s\) through \[ x = x_0 + s h, \] so that \(s\) counts nodes rather than measuring distance.

Then \(x - x_i = (s - i)h\), and every product in the Newton form becomes a product of integers times a power of \(h\): \[ \begin{aligned} P_n(x) = P_n(x_0 + sh) &= f[x_0] + sh\,f[x_0,x_1] + s(s-1)h^2 f[x_0,x_1,x_2] \\ &\quad + \cdots + s(s-1)\cdots(s-n+1)h^n f[x_0,\ldots,x_n] \\ &= f[x_0] + \sum_{k=1}^n \binom{s}{k} k!\, h^k f[x_0,\ldots,x_k], \end{aligned} \] where \[ \binom{s}{k} = \frac{s(s-1)\cdots(s-k+1)}{k!}. \]

The forward difference notation removes the factorials and powers of \(h\) altogether. Writing \(\Delta f(x_i) = f(x_{i+1}) - f(x_i)\), \[ f[x_0,x_1] = \frac{f(x_1) - f(x_0)}{x_1 - x_0} = \frac{1}{h}\,\Delta f(x_0), \] \[ f[x_0,x_1,x_2] = \frac{1}{2h}\cdot\frac{\Delta f(x_1) - \Delta f(x_0)}{h} = \frac{1}{2h^2}\,\Delta^2 f(x_0), \] and in general \[ f[x_0,\ldots,x_k] = \frac{1}{k!\,h^k}\,\Delta^k f(x_0). \]

Substituting back, the \(k!\) and \(h^k\) cancel and the Newton forward-difference formula is simply \[ P_n(x) = f(x_0) + \sum_{k=1}^n \binom{s}{k}\Delta^k f(x_0). \]

Example 3.4

A forward-difference table.

Interpolate the equally spaced data \[ (x_i, f(x_i)) = (0,1),\ \left(\tfrac13, 2\right),\ \left(\tfrac23, 1\right),\ (1,1), \] so \(h = \tfrac13\) and \(s = x/h = 3x\).

Differencing repeatedly:

\(f\)	\(\Delta f\)	\(\Delta^2 f\)	\(\Delta^3 f\)
\(1\)
	\(1\)
\(2\)		\(-2\)
	\(-1\)		\(3\)
\(1\)		\(1\)
	\(0\)
\(1\)

The top of each column supplies a coefficient: \[ P_3(s) = 1 + \binom{s}{1}\cdot 1 + \binom{s}{2}\cdot(-2) + \binom{s}{3}\cdot 3 . \] Expanding the binomial coefficients, \[ P_3(s) = 1 + s + \frac{s(s-1)}{2}(-2) + \frac{s(s-1)(s-2)}{3!}\cdot 3 = 1 + s - s(s-1) + \frac{s(s-1)(s-2)}{2}. \] Finally, substituting \(s = 3x\), \[ P_3(x) = 1 + 3x - 3x(3x-1) + \frac{3x(3x-1)(3x-2)}{2}. \]

3.4 Hermite Interpolation

So far we have seen two ways to make a polynomial agree with \(f\). Taylor’s theorem matches many derivatives at a single point; Lagrange interpolation matches the value at many points. Each is the extreme case of a single idea.

Osculating polynomials are the general case: at each of several nodes, match the value and the first few derivatives. Taylor is the case of one node with many derivatives; Lagrange is the case of many nodes with no derivatives. The case that matters in practice is the one in between — match the value and the first derivative at every node — and that is Hermite interpolation.

Example 3.5

A Hermite cubic by brute force.

Find \(P\) with \[ P(0) = 0, \quad P'(0) = 0, \quad P(1) = 1, \quad P'(1) = 0 . \]

Four conditions determine four coefficients, so we look for a cubic, \[ P_3(x) = ax^3 + bx^2 + cx + d, \qquad P_3'(x) = 3ax^2 + 2bx + c . \]

Imposing the conditions in turn: \[ \begin{aligned} P_3(0) &= d = 0, \\ P_3'(0) &= c = 0, \\ P_3(1) &= a + b + c + d = 1 &&\Longrightarrow\quad a + b = 1, \\ P_3'(1) &= 3a + 2b + c = 0 &&\Longrightarrow\quad 3a + 2b = 0 . \end{aligned} \] The last two give \(-b = -3\), so \(b = 3\) and \(a = -2\). Hence \[ P_3(x) = -2x^3 + 3x^2 . \]

Hermite polynomials from a basis

Solving a linear system each time is unnecessary: as with Lagrange, there is a basis that makes the conditions automatic. Now we need two families — one to carry the values and one to carry the derivatives.

Example 3.6

The Hermite basis for \(n = 1\).

Take \(x_0 = 0\) and \(x_1 = 1\), so the Lagrange basis is \[ L_{1,0}(x) = 1 - x, \quad L_{1,0}' = -1, \qquad L_{1,1}(x) = x, \quad L_{1,1}' = 1 . \]

Substituting into the definitions from the theorem, \[ \begin{aligned} H_{1,0}(x) &= \bigl[1 - 2x(-1)\bigr](1-x)^2 = (1 + 2x)(1-x)^2, \\ H_{1,1}(x) &= \bigl[1 - 2(x-1)(1)\bigr]x^2 = (3 - 2x)x^2, \\ \hat H_{1,0}(x) &= x(1-x)^2, \\ \hat H_{1,1}(x) &= (x-1)x^2 . \end{aligned} \]

Hermite polynomials from divided differences

There is a much slicker route. Build a divided-difference table in which every node is repeated, listing \(z_{2i} = z_{2i+1} = x_i\). The repeated entries would give a zero denominator, so wherever that happens we use the limiting value from the previous section, \[ f[x_i, x_i] = f'(x_i), \] and the ordinary Newton formula then produces the Hermite polynomial with no new machinery at all.

Example 3.7

A Hermite cubic by divided differences.

Interpolate \(f(0) = 0\), \(f'(0) = 0\), \(f(1) = 0\), \(f'(1) = 1\).

List each node twice, \(z = 0, 0, 1, 1\):

\(z_i\)	\(f(z_i)\)	1st DD	2nd DD	3rd DD
\(0\)	\(0\)
		\(f[0,0] = f'(0) = 0\)
\(0\)	\(0\)		\(f[0,0,1] = \frac{0-0}{1-0} = 0\)
		\(f[0,1] = \frac{0-0}{1-0} = 0\)		\(f[0,0,1,1] = \frac{1-0}{1-0} = 1\)
\(1\)	\(0\)		\(f[0,1,1] = \frac{1-0}{1-0} = 1\)
		\(f[1,1] = f'(1) = 1\)
\(1\)	\(0\)

The two entries marked \(f[0,0]\) and \(f[1,1]\) are where the derivative data enters. Reading off the coefficients, \[ H_3(x) = f[z_0] + \sum_{k=1}^{3} f[z_0,\ldots,z_k](x - z_0)\cdots(x - z_{k-1}), \] which here is \[ H_3(x) = 0 + 0\,(x-0) + 0\,(x-0)(x-0) + 1\,(x-0)^2(x-1) = x^2(x-1). \]

Example 3.8

A Hermite quintic.

Now take three nodes, with \[ f(0) = 0,\ f'(0) = 0, \qquad f(1) = 1,\ f'(1) = -1, \qquad f(2) = 1,\ f'(2) = 0 . \]

Six conditions give a quintic. Listing \(z = 0,0,1,1,2,2\) and running the table produces the leading coefficients \[ 0,\quad 0,\quad 1,\quad -3,\quad \tfrac94,\quad -\tfrac74, \] so that \[ H_5(x) = 0 + 0\,x + 1\,x^2 - 3\,x^2(x-1) + \tfrac94\,x^2(x-1)^2 - \tfrac74\,x^2(x-1)^2(x-2), \] which simplifies to \[ H_5(x) = -\frac{x^2}{4}\left(7x^3 - 37x^2 + 65x - 39\right). \]

3.5 Cubic Spline Interpolation

A single interpolating polynomial through many nodes is usually a bad idea. For \(f(x) = \dfrac{1}{1+x^2}\) with the eleven equally spaced nodes \(x_0, \ldots, x_{10} = -5, -4, \ldots, 5\), the degree-10 interpolant is wildly oscillatory near the ends of the interval — it matches the data and misses the function.

The alternatives are all piecewise:

1. Piecewise linear — join the dots. Interpolates the data, but is not smooth: the slope jumps at every node.
2. Piecewise Hermite (cubic) — smooth, but it needs the derivative values \(f'(x_i)\), which the data usually does not supply.
3. Cubic spline — a piecewise cubic whose zeroth, first and second derivatives are all continuous, using only the values \(f(x_i)\).

The third gets smoothness for free by solving for the derivative information rather than demanding it as input.

The conditions defining a cubic spline \(S\) on nodes \(x_0 < \cdots < x_n\) are:

1. \(S(x)\) is a cubic polynomial \(S_j\) on each \([x_j, x_{j+1}]\);
2. \(S_j(x_j) = f(x_j)\) and \(S_j(x_{j+1}) = f(x_{j+1})\) — it interpolates;
3. \(S_{j+1}(x_{j+1}) = S_j(x_{j+1})\) — the pieces agree in value;
4. \(S_{j+1}'(x_{j+1}) = S_j'(x_{j+1})\) — and in slope;
5. \(S_{j+1}''(x_{j+1}) = S_j''(x_{j+1})\) — and in curvature;
6. one of two sets of boundary conditions holds: natural (or free), \(S''(x_0) = S''(x_n) = 0\), or clamped, \(S'(x_0) = f'(x_0)\) and \(S'(x_n) = f'(x_n)\).

Counting degrees of freedom

It is worth checking that this is neither under- nor over-determined. With \(n+1\) points there are \(n\) intervals, hence \(n\) cubics and \[ 4n \text{ coefficients.} \]

Counting equations:

Exactly \(4n\) conditions for \(4n\) unknowns. Condition (c) is worth a second look: it costs nothing, because (b) already forces both pieces meeting at \(x_{j+1}\) to take the value \(f(x_{j+1})\) there.

Example 3.9

Solving for the coefficients directly.

Find the natural cubic spline through \((0,0)\), \((1,1)\), \((2,1)\).

Two intervals means two cubics, written about the left end of each: \[ \begin{aligned} S_0(x) &= a_0 + b_0 x + c_0 x^2 + d_0 x^3, \\ S_1(x) &= a_1 + b_1(x-1) + c_1(x-1)^2 + d_1(x-1)^3, \end{aligned} \] so there are eight coefficients to find.

(b) Interpolation. \[ \begin{aligned} S_0(0) &= a_0 = 0, \\ S_0(1) &= a_0 + b_0 + c_0 + d_0 = 1, \\ S_1(1) &= a_1 = 1, \\ S_1(2) &= a_1 + b_1 + c_1 + d_1 = 1 . \end{aligned} \]

(d) Continuity of the first derivative. Differentiating, \[ S_0'(x) = b_0 + 2c_0 x + 3d_0 x^2, \qquad S_1'(x) = b_1 + 2c_1(x-1) + 3d_1(x-1)^2, \] and matching at \(x=1\) gives \[ b_0 + 2c_0 + 3d_0 = b_1 . \]

(e) Continuity of the second derivative. Similarly \[ S_0''(x) = 2c_0 + 6d_0 x, \qquad S_1''(x) = 2c_1 + 6d_1(x-1), \] so matching at \(x=1\) gives \[ 2c_0 + 6d_0 = 2c_1 . \]

(f) Natural boundary conditions. \[ S_0''(0) = 2c_0 = 0, \qquad S_1''(2) = 2c_1 + 6d_1 = 0 . \]

That is eight linear equations in eight unknowns, which we hand to MATLAB as a single system \(A\mathbf{c} = \mathbf{b}\) with \[ \mathbf{c} = (a_0, b_0, c_0, d_0, a_1, b_1, c_1, d_1)^T . \] The solution is \[ (a_0, b_0, c_0, d_0, a_1, b_1, c_1, d_1) = \left(0, \tfrac54, 0, -\tfrac14, 1, \tfrac12, -\tfrac34, \tfrac14\right), \] giving \[ \begin{cases} S_0(x) = \dfrac54 x - \dfrac14 x^3, & 0 \le x \le 1, \\[6pt] S_1(x) = 1 + \dfrac12 (x-1) - \dfrac34 (x-1)^2 + \dfrac14 (x-1)^3, & 1 \le x \le 2 . \end{cases} \]

Note how the natural boundary condition shows up in the answer: \(c_0 = 0\) is exactly \(S_0''(0) = 0\).

Solving an \(8 \times 8\) system by hand does not scale. The systematic version reduces the problem to a tridiagonal system in the second derivatives alone, which is what the algorithm on the slides does — and being tridiagonal, it costs \(O(n)\) rather than \(O(n^3)\) to solve.