8: Interpolation

Consider the following problem: Given the values of a known function \(y=f(x)\) at a sequence of ordered points \(x_{0}, x_{1}, \ldots, x_{n}\), find \(f(x)\) for arbitrary \(x .\) When \(x_{0} \leq x \leq x_{n}\), the problem is called interpolation. When \(x<x_{0}\) or \(x>x_{n}\), the problem is called extrapolation.

With \(y_{i}=f\left(x_{i}\right)\), the problem of interpolation is basically one of drawing a smooth curve through the known points \(\left(x_{0}, y_{0}\right),\left(x_{1}, y_{1}\right), \ldots,\left(x_{n}, y_{n}\right) .\) This is not the same problem as drawing a smooth curve that approximates a set of data points that have experimental error. This latter problem is called least-squares approximation, which is considered in the next chapter.

It is possible to interpolate the \(n+1\) known points by a unique polynomial of degree \(n\). When \(n=1\), the polynomial is a linear function; when \(n=2\), the polynomial is a quadratic function. Although low order polynomials are sometimes used when the number of points are few, higher-order polynomial interpolates tend to be unstable, and are not of much practical use.

Here, we will consider the more useful piece-wise polynomial interpolation. The two most used are piecewise linear interpolation, and cubic spline interpolation. The first makes use of linear polynomials, and the second cubic polynomials.

Piecewise linear interpolation

Here, we use linear polynomials. This is the default interpolation typically used when plotting data.

Suppose the interpolating function is \(y=g(x)\), and as previously, there are \(n+1\) points to interpolate. We construct the function \(g(x)\) out of \(n\) local linear polynomials. We write

\[g(x)=g_{i}(x), \quad \text { for } x_{i} \leq x \leq x_{i+1}, \nonumber \]

where

\[g_{i}(x)=a_{i}\left(x-x_{i}\right)+b_{i} \nonumber \]

and \(i=0,1, \ldots, n-1\)

We now require \(y=g_{i}(x)\) to pass through the endpoints \(\left(x_{i}, y_{i}\right)\) and \(\left(x_{i+1}, y_{i+1}\right) .\) We have

\[\begin{aligned} y_{i} &=b_{i} \\ y_{i+1} &=a_{i}\left(x_{i+1}-x_{i}\right)+b_{i} \end{aligned} \nonumber \]

Therefore, the coefficients of \(g_{i}(x)\) are determined to be

\[a_{i}=\frac{y_{i+1}-y_{i}}{x_{i+1}-x_{i}}, \quad b_{i}=y_{i} \nonumber \]

Although piecewise linear interpolation is widely used, particularly in plotting routines, it suffers from a discontinuity in the derivative at each point. This results in a function which may not look smooth if the points are too widely spaced. We next consider a more challenging algorithm that uses cubic polynomials.

Cubic spline interpolation

The \(n+1\) points to be interpolated are again

\[\left(x_{0}, y_{0}\right),\left(x_{1}, y_{1}\right), \ldots\left(x_{n}, y_{n}\right) . \nonumber \]

Here, we use \(n\) piecewise cubic polynomials for interpolation,

\[g_{i}(x)=a_{i}\left(x-x_{i}\right)^{3}+b_{i}\left(x-x_{i}\right)^{2}+c_{i}\left(x-x_{i}\right)+d_{i}, \quad i=0,1, \ldots, n-1, \nonumber \]

with the global interpolation function written as

\[g(x)=g_{i}(x), \quad \text { for } x_{i} \leq x \leq x_{i+1} . \nonumber \]

To achieve a smooth interpolation we impose that \(g(x)\) and its first and second derivatives are continuous. The requirement that \(g(x)\) is continuous (and goes through all \(n+1\) points) results in the two constraints

\[\begin{align} g_{i}\left(x_{i}\right) &=y_{i}, \quad i=0 \text { to } n-1 \\ g_{i}\left(x_{i+1}\right) &=y_{i+1}, \quad i=0 \text { to } n-1 \end{align} \nonumber \]

The requirement that \(g^{\prime}(x)\) is continuous results in

\[g_{i}^{\prime}\left(x_{i+1}\right)=g_{i+1}^{\prime}\left(x_{i+1}\right), \quad i=0 \text { to } n-2 \nonumber \]

And the requirement that \(g^{\prime \prime}(x)\) is continuous results in

\[g_{i}^{\prime \prime}\left(x_{i+1}\right)=g_{i+1}^{\prime \prime}\left(x_{i+1}\right), \quad i=0 \text { to } n-2 . \nonumber \]

There are \(n\) cubic polynomials \(g_{i}(x)\) and each cubic polynomial has four free coefficients; there are therefore a total of \(4 n\) unknown coefficients. The number of constraining equations from (8.7)-(8.10) is \(2 n+2(n-1)=4 n-2\). With \(4 n-2\) constraints and \(4 n\) unknowns, two more conditions are required for a unique solution. These are usually chosen to be extra conditions on the first \(g_{0}(x)\) and last \(g_{n-1}(x)\) polynomials. We will discuss these extra conditions later.

We now proceed to determine equations for the unknown coefficients of the cubic polynomials. The polynomials and their first two derivatives are given by

\[\begin{align} g_{i}(x) &=a_{i}\left(x-x_{i}\right)^{3}+b_{i}\left(x-x_{i}\right)^{2}+c_{i}\left(x-x_{i}\right)+d_{i} \\ g_{i}^{\prime}(x) &=3 a_{i}\left(x-x_{i}\right)^{2}+2 b_{i}\left(x-x_{i}\right)+c_{i} \\ g_{i}^{\prime \prime}(x) &=6 a_{i}\left(x-x_{i}\right)+2 b_{i} \end{align} \nonumber \]

We will consider the four conditions (8.7)-(8.10) in turn. From (8.7) and (8.11), we have

\[d_{i}=y_{i}, \quad i=0 \text { to } n-1, \nonumber \]

which directly solves for all of the \(d\)-coefficients.

To satisfy (8.8), we first define

\[h_{i}=x_{i+1}-x_{i}, \nonumber \]

and

\[f_{i}=y_{i+1}-y_{i} . \nonumber \]

Now, from (8.8) and (8.11), using (8.14), we obtain the \(n\) equations

\[a_{i} h_{i}^{3}+b_{i} h_{i}^{2}+c_{i} h_{i}=f_{i}, \quad i=0 \text { to } n-1 . \nonumber \]

From (8.9) and (8.12) we obtain the \(n-1\) equations

\[3 a_{i} h_{i}^{2}+2 b_{i} h_{i}+c_{i}=c_{i+1}, \quad i=0 \text { to } n-2 \nonumber \]

From (8.10) and (8.13) we obtain the \(n-1\) equations

\[3 a_{i} h_{i}+b_{i}=b_{i+1} \quad i=0 \text { to } n-2 \text {. } \nonumber \]

It will be useful to include a definition of the coefficient \(b_{n}\), which is now missing. (The index of the cubic polynomial coefficients only go up to \(n-1\).) We simply extend (8.19) up to \(i=n-1\) and so write

\[3 a_{n-1} h_{n-1}+b_{n-1}=b_{n}, \nonumber \]

which can be viewed as a definition of \(b_{n}\).

We now proceed to eliminate the sets of \(a\) - and c-coefficients (with the \(d\)-coefficients already eliminated in (8.14)) to find a system of linear equations for the \(b\)-coefficients. From (8.19) and \((8.20)\), we can solve for the \(n a\)-coefficients to find

\[a_{i}=\frac{1}{3 h_{i}}\left(b_{i+1}-b_{i}\right), \quad i=0 \text { to } n-1 . \nonumber \]

From (8.17), we can solve for the \(n\) c-coefficients as follows:

\[\begin{align} \nonumber c_{i} &=\frac{1}{h_{i}}\left(f_{i}-a_{i} h_{i}^{3}-b_{i} h_{i}^{2}\right) \\ \nonumber &=\frac{1}{h_{i}}\left(f_{i}-\frac{1}{3 h_{i}}\left(b_{i+1}-b_{i}\right) h_{i}^{3}-b_{i} h_{i}^{2}\right) \\ &=\frac{f_{i}}{h_{i}}-\frac{1}{3} h_{i}\left(b_{i+1}+2 b_{i}\right), \quad i=0 \text { to } n-1 \end{align} \nonumber \]

We can now find an equation for the \(b\)-coefficients by substituting \((8.21)\) and \((8.22)\) into \((8.18)\) :

\[\begin{array}{r} 3\left(\frac{1}{3 h_{i}}\left(b_{i+1}-b_{i}\right)\right) h_{i}^{2}+2 b_{i} h_{i}+\left(\frac{f_{i}}{h_{i}}-\frac{1}{3} h_{i}\left(b_{i+1}+2 b_{i}\right)\right) \\ =\left(\frac{f_{i+1}}{h_{i+1}}-\frac{1}{3} h_{i+1}\left(b_{i+2}+2 b_{i+1}\right)\right), \nonumber \end{array} \nonumber \]

which simplifies to

\[\frac{1}{3} h_{i} b_{i}+\frac{2}{3}\left(h_{i}+h_{i+1}\right) b_{i+1}+\frac{1}{3} h_{i+1} b_{i+2}=\frac{f_{i+1}}{h_{i+1}}-\frac{f_{i}}{h_{i}} \nonumber \]

an equation that is valid for \(i=0\) to \(n-2\). Therefore, (8.23) represent \(n-1\) equations for the \(n+1\) unknown \(b\)-coefficients. Accordingly, we write the matrix equation for the \(b\)-coefficients, leaving the first and last row absent, as

\[\left(\begin{array}{cccccccc} \ldots & \ldots & \ldots & \ldots & \text { missing } & \ldots & \ldots \\ \frac{1}{3} h_{0} & \frac{2}{3}\left(h_{0}+h_{1}\right) & \frac{1}{3} h_{1} & \ldots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \ldots & \frac{1}{3} h_{n-2} & \frac{2}{3}\left(h_{n-2}+h_{n-1}\right) & \frac{1}{3} h_{n-1} \\ \ldots & \ldots & \ldots & \ldots & \text { missing } & \ldots & \ldots \end{array}\right)\left(\begin{array}{c} b_{0} \\ b_{1} \\ \vdots \\ b_{n-1} \\ b_{n} \end{array}\right)\left(\begin{array}{c} \text { missing } \\ \frac{f_{1}}{h_{1}}-\frac{f_{0}}{h_{0}} \\ \vdots \\ \frac{f_{n-1}}{h_{n-1}}-\frac{f_{n-2}}{h_{n-2}} \\ \text { missing } \end{array}\right) . \nonumber \]

Once the missing first and last equations are specified, the matrix equation for the \(b\)-coefficients can be solved by Gaussian elimination. And once the \(b\)-coefficients are determined, the \(a\) - and \(c\)-coefficients can also be determined from (8.21) and (8.22), with the \(d\)-coefficients already known from (8.14). The piecewise cubic polynomials, then, are known and \(g(x)\) can be used for interpolation to any value \(x\) satisfying \(x_{0} \leq x \leq x_{n}\).

The missing first and last equations can be specified in several ways, and here we show the two ways that are allowed by the MATLAB function spline.m. The first way should be used when the derivative \(g^{\prime}(x)\) is known at the endpoints \(x_{0}\) and \(x_{n} ;\) that is, suppose we know the values of \(\alpha\) and \(\beta\) such that

\[g_{0}^{\prime}\left(x_{0}\right)=\alpha, \quad g_{n-1}^{\prime}\left(x_{n}\right)=\beta . \nonumber \]

From the known value of \(\alpha\), and using (8.12) and (8.22), we have

\[\begin{aligned} \alpha &=c_{0} \\ &=\frac{f_{0}}{h_{0}}-\frac{1}{3} h_{0}\left(b_{1}+2 b_{0}\right) \end{aligned} \nonumber \]

Therefore, the missing first equation is determined to be

\[\frac{2}{3} h_{0} b_{0}+\frac{1}{3} h_{0} b_{1}=\frac{f_{0}}{h_{0}}-\alpha \nonumber \]

From the known value of \(\beta\), and using \((8.12),(8.21)\), and \((8.22)\), we have

\[\begin{aligned} \beta &=3 a_{n-1} h_{n-1}^{2}+2 b_{n-1} h_{n-1}+c_{n-1} \\ &=3\left(\frac{1}{3 h_{n-1}}\left(b_{n}-b_{n-1}\right)\right) h_{n-1}^{2}+2 b_{n-1} h_{n-1}+\left(\frac{f_{n-1}}{h_{n-1}}-\frac{1}{3} h_{n-1}\left(b_{n}+2 b_{n-1}\right)\right), \end{aligned} \nonumber \]

which simplifies to

\[\frac{1}{3} h_{n-1} b_{n-1}+\frac{2}{3} h_{n-1} b_{n}=\beta-\frac{f_{n-1}}{h_{n-1}} \nonumber \]

to be used as the missing last equation.

The second way of specifying the missing first and last equations is called the not-a-knot condition, which assumes that

\[g_{0}(x)=g_{1}(x), \quad g_{n-2}(x)=g_{n-1}(x) \nonumber \]

Considering the first of these equations, from (8.11) we have

\[a_{0}\left(x-x_{0}\right)^{3}+b_{0}\left(x-x_{0}\right)^{2}+c_{0}\left(x-x_{0}\right)+d_{0} \nonumber \]

\[=a_{1}\left(x-x_{1}\right)^{3}+b_{1}\left(x-x_{1}\right)^{2}+c_{1}\left(x-x_{1}\right)+d_{1} . \nonumber \]

Now two cubic polynomials can be proven to be identical if at some value of \(x\), the polynomials and their first three derivatives are identical. Our conditions of continuity at \(x=x_{1}\) already require that at this value of \(x\) these two polynomials and their first two derivatives are identical. The polynomials themselves will be identical, then, if their third derivatives are also identical at \(x=x_{1}\), or if

\[a_{0}=a_{1} . \nonumber \]

From (8.21), we have

\[\frac{1}{3 h_{0}}\left(b_{1}-b_{0}\right)=\frac{1}{3 h_{1}}\left(b_{2}-b_{1}\right), \nonumber \]

or after simplification

\[h_{1} b_{0}-\left(h_{0}+h_{1}\right) b_{1}+h_{0} b_{2}=0, \nonumber \]

which provides us our missing first equation. A similar argument at \(x=x_{n-1}\) also provides us with our last equation,

\[h_{n-1} b_{n-2}-\left(h_{n-2}+h_{n-1}\right) b_{n-1}+h_{n-2} b_{n}=0 . \nonumber \]

The MATLAB subroutines spline.m and ppval.m can be used for cubic spline interpolation (see also interp1.m). I will illustrate these routines in class and post sample code on the course web site.

Multidimensional interpolation

Suppose we are interpolating the value of a function of two variables,

\[z=f(x, y) . \nonumber \]

The known values are given by

\[z_{i j}=f\left(x_{i}, y_{j}\right), \nonumber \]

with \(i=0,1, \ldots, n\) and \(j=0,1, \ldots, n .\) Note that the \((x, y)\) points at which \(f(x, y)\) are known lie on a grid in the \(x-y\) plane.

Let \(z=g(x, y)\) be the interpolating function, satisfying \(z_{i j}=g\left(x_{i}, y_{j}\right) .\) A two-dimensional interpolation to find the value of \(g\) at the point \((x, y)\) may be done by first performing \(n+1\) one-dimensional interpolations in \(y\) to find the value of \(g\) at the \(n+1\) points \(\left(x_{0}, y\right),\left(x_{1}, y\right), \ldots\), \(\left(x_{n}, y\right)\), followed by a single one-dimensional interpolation in \(x\) to find the value of \(g\) at \((x, y)\).

In other words, two-dimensional interpolation on a grid of dimension \((n+1) \times(n+1)\) is done by first performing \(n+1\) one-dimensional interpolations to the value \(y\) followed by a single one-dimensional interpolation to the value \(x\). Two-dimensional interpolation can be generalized to higher dimensions. The MATLAB functions to perform two-and three-dimensional interpolation are interp2.m and interp3.m.