5.3: 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{gathered} 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{gathered} \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 (5.1)-(5.4) 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{aligned} 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{aligned} \nonumber \]

We will consider the four conditions (5.1)-(5.4) in turn. From (5.1) and (5.5), 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 \((5.2)\), we first define

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

and

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

Now, from (5.2) and (5.5), using (5.8), 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 (5.3) and (5.6) 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 (5.4) and (5.7) we obtain the \(n-1\) equations

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

It is 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 (5.11) 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 (5.8)) to find a system of linear equations for the \(b\)-coefficients. From (5.11) and (5.12), 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 (5.9), we can solve for the \(n\) c-coefficients as follows:

\[\begin{aligned} c_{i} &=\frac{1}{h_{i}}\left(f_{i}-a_{i} h_{i}^{3}-b_{i} h_{i}^{2}\right) \\ &=\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{aligned} \nonumber \]

We can now find an equation for the \(b\)-coefficients by substituting (5.8), (5.13) and \((5.14)\) into \((5.10)\) :

\[\begin{gathered} 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) \end{gathered} \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, (5.15) represent \(n-1\) equations for the \(n+1\) unknown \(b\)-coefficients. Accordingly, we write the matrix equa- tion for the \(b\)-coefficients, leaving the first and last row absent, as

\[\left(\begin{array}{cccccccc} \cdots & \cdots & \cdots & \cdots & \text { missing } & \cdots & \cdots \\ \frac{1}{3} h_{0} & \frac{2}{3}\left(h_{0}+h_{1}\right) & \frac{1}{3} h_{1} & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & \frac{1}{3} h_{n-2} & \frac{2}{3}\left(h_{n-2}+h_{n-1}\right) & \frac{1}{3} h_{n-1} \\ \cdots & \cdots & \cdots & \cdots & \text { missing } & \cdots & \cdots \end{array}\right)\left(\begin{array}{c} b_{0} \\ b_{1} \\ \vdots \\ b_{n-1} \\ b_{n} \end{array}\right)\left(\begin{array}{c} \operatorname{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}} \\ \operatorname{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 (5.13) and (5.14), with the \(d\)-coefficients already known from (5.8). 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 \((5.6)\) and \((5.14)\), 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 \((5.6),(5.13)\), and \((5.14)\), 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 (5.5) we have

\[\begin{aligned} a_{0}\left(x-x_{0}\right)^{3}+b_{0}\left(x-x_{0}\right)^{2}+c_{0}(x&\left.-x_{0}\right)+d_{0} \\ &=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} \end{aligned} \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 (5.13), 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.