Spline interpolation

GUIDING QUESTION: How can I more accurately interpolate a given set of data points?

Motivation


Polynomial interpolation can be efficient (recall the barycentric formula), but using high degree polynomials can lead to large errors due to erratic oscillations, especially near the interval endpoints.

Piecewise linear interpolation

For starters, let's consider the simplest case.

Insight

We have constructed a linear interpolant $s_i(x)$ on each interval $[x_i, x_{i+1}]$, for $i=0,1, \ldots, n-1$.
On each subinterval $[x_i, x_{i+1}]$, the piecewise polynomial interpolant $s$ coincides with $s_i$. When $x_i \leq x \leq x_{i+1}$, $$s(x)= s_i(x)= \color{var(--emphColor)}{a_i} + \color{var(--emphColor)}{b_i} (x-x_i)$$
From the point-slope formula it follows that $$\color{var(--emphColor)}{a_i} = f(x_i) \text{ and } \color{var(--emphColor)}{b_i} = \frac{f(x_{i+1}) - f(x_i)}{x_{i+1} - x_i}.$$
We break up the interval $[a, b] = [x_0, x_n]$ into $n$ pieces $[x_i, x_{i+1}]$, for $i = 0, 1, \ldots, n-1$.
So the piecewise interpolant $s(x)$ on $[a, b]$ is defined by $n$ pieces $s_i(x)$, for $i = 0, 1, \ldots, n-1$.
For convenience, we define $s_n(x_n) = s_{n-1}(x_n)$.

Conditions

Two sets of conditions define the $s_i$.
  1. (Interpolation) $s(x_i) = s_i(x_i) = f(x_i)$, for $i = 0, 1, \ldots, n$.
  2. (Continuity) $s_i(x_{i+1}) = s_{i+1}(x_{i+1})$, for $i = 0, 1, \ldots, n-2$.
Each piece $s_i(x)= \color{var(--emphColor)}{a_i} + \color{var(--emphColor)}{b_i} (x-x_i)$ depends on two coefficients, so we must solve for $2n$ coefficients in total.
We have imposed $2n$ conditions to solve for theses $2n$ coefficients.
MATLAB's $\texttt{plot}$ function uses piecewise linear interpolation to display continous graphs when given only a finite set of points.

Cubic spline

Let's use the insight we extracted when building piecewise linear interpolants to construct piecewise cubic interpolants, or cubic splines.
We're skipping to cubic interpolants because you'll develop quadratic splines in HW5.
In this case, we require that the piecewise interpolant $$ s(x)= s_i(x) = \color{var(--emphColor)}{a_i} + \color{var(--emphColor)}{b_i} (x-x_i) + \color{var(--emphColor)}{c_i} (x-x_i)^2 + \color{var(--emphColor)}{d_i}(x-x_i)^3$$ on each subinterval $[x_i,x_{i+1}]$, where $\color{var(--emphColor)}{a_i}, \color{var(--emphColor)}{b_i}, \color{var(--emphColor)}{c_i}$, and $\color{var(--emphColor)}{d_i}$ are coefficients to be determined.

Cubic spline

$s(x)= s_i(x) = \color{var(--emphColor)}{a_i} + \color{var(--emphColor)}{b_i} (x-x_i) + \color{var(--emphColor)}{c_i} (x-x_i)^2 + \color{var(--emphColor)}{d_i}(x-x_i)^3$

As in the piecewise linear case, we require
  1. (Interpolation) $s(x_i) = s_i(x_i) = f(x_i)$, for $i = 0, 1, \ldots, n$.
  2. (Continuity) $s_i(x_\color{var(--emphColor)}{i+1}) = s_{i+1}(x_\color{var(--emphColor)}{i+1})$, for $i = 0, 1, \ldots, n-2$.
In this case there are $4n$ coefficients in total, so our best bet is to specify $4n$ constraints.

Smoothness

The piecewise linear interpolant can have sharp corners, making it non-differeniable.
This is undesirable, so we require smoothness for the cubic spline.
  1. (Continuity of $s'$) $s_i'(x_\color{var(--emphColor)}{i+1}) = s_{i+1}'(x_\color{var(--emphColor)}{i+1})$, for $i = 0, 1, \ldots, n-2$.
  2. (Continuity of $s''$) $s_i''(x_\color{var(--emphColor)}{i+1}) = s_{i+1}''(x_\color{var(--emphColor)}{i+1})$, for $i = 0, 1, \ldots, n-2$.

Insufficient requirements

The interpolation and coninuity constraints amount to $$n+ 1 + 3(n-1) = 4n - 2$$ conditions, so we must specify two more conditions in order to hope we have defined the cubic spline uniquely.
Before discussing these additional constraints, let's work with what we've got!

Manipulation

We can manipulate the given constraints to obtain formulas for $a_i$, formulas for $b_i$ and $d_i$ in terms of the $c_i$, and we can set up a linear system to solve for the $c_i$.

Continuity of $s''$

Let $h_i = x_{i+1} - x_i$. Then continuity of $s''$ means $$2c_ih_i +3 \cdot 2\color{var(--emphColor)}{d_i}h_i=s''_i(x_{i+1}) = s''_{i+1}(x_{i+1})=2 c_{i+1},$$ for $i=0,1, \ldots, n-2$.
We may solve for $\color{var(--emphColor)}{d_i}$ to obtain $$\color{var(--emphColor)}{d_i} = \frac{c_{i+1} - c_i}{3h_i}.$$

Continuity of $s'$

Substituting $\color{var(--emphColor)}{d_i}$ into the constraints requiring continuity of $s'$, we obtain \begin{align} b_{i+1} &= b_i + 2c_ih_i + (c_{i+1} - c_i) h_i \\ &= b_i+(c_{i+1} + c_i) h_i. \end{align}

Continuity of $s$

Similarly, substituting $\color{var(--emphColor)}{d_i}$ into the constraints requiring continuity of $s$, we obtain \begin{align} a_{i+1}&= a_i + \color{var(--emphColor)}{b_i} h_i + c_ih_i^2 + \frac{c_{i+1} - c_i}{3}h_i^2\\ &= a_i + \color{var(--emphColor)}{b_i} h_i + \frac{c_{i+1} +2c_i}{3}h_i^2. \end{align}
We may solve the last equation for $\color{var(--emphColor)}{b_i}$ to conclude that $$\color{var(--emphColor)}{b_i} = \frac{a_{i+1} - a_i}{h_i} -\frac{2c_i +c_{i+1}}{3}h_i.$$

Interpolation

The interpolation conditions imply \begin{align} f(x_i) &= s_i (x_i) \\ &= \color{var(--emphColor)}{a_i} + b_i(x_i-x_i) + c_i (x_i-x_i)^2+ d_i (x_i-x_i)^3\\ &= \color{var(--emphColor)}{a_i} \end{align} for $i=0, 1, \ldots, n$.
For convenience, we write $a_n = f(x_n)$.
At this point, all the coefficients are either known or they are given in terms of the $c_i$!

Return of the linear system

After doing some more algebra we find simple linear relations defining the $c_i$'s: $$c_{i-1} + 4 c_i + c_{i+1} = \frac{3}{h^2} (a_{i+1} - 2a_i + a_{i-1}),$$ for $i = 1, \ldots, n-2$, assuming the nodes $x_i$ are evenly spaced, so that $h_0 = h_1 = \cdots = h_{n-1} = h$.
Please refer to the supplementary notes or to Bradie's textbook for more details.

Boundary conditions

At this point, we're ready to introduce the 2 additional conditions needed to obtain with a square system.

Not-a-knot

We require continuity of $s'''$ at $x = x_1$, and at $x = x_{n-1}$.
Since $s'''_i(x) = 6d_i$, the not-a-knot boundary conditions require $$d_0 = d_1, \quad \text{and} \quad d_{n-2} = d_{n-1}.$$
We may obtain equations for $c_i$ using $d_i = \frac{c_{i+1}-c_i}{3h_i}$.
We may then solve for $c_0$ and $c_n$ to preserve the tridiagonal structure of the system.

Not-a-knot

When all is said and done, we obtain an $(n-1) \times (n-1)$ tridiagonal system defining the $c_i$: \begin{align} 6c_1 \quad \quad \quad \quad &=\frac{3}{h^2}(a_2 -2a_1 +a_0)\\ c_{i-1} + 4 c_i +c_{i+1} &= \frac{3}{h^2} (a_{i+1} - 2a_i + a_{i-1}), \quad i=2, \ldots, n-2\\ 6c_{n-1} &= \frac{3}{h^2} (a_n -2a_{n-1} +a_{n-2}) \end{align}
We have again assumed the nodes are equidistant. See the supplementary notes for more details.
When no additional information about $f$ is given, the not-a-knot boundary conditions are recommended.

Clamped spline

Imagine placing a clamp at $x = a$ and at $x = b$. This would force your piecewise polynomial interpolant to have the same initial and final slope as $f$.
If you know $f'(a)$ and $f'(b)$, use them!

Clamped spline

In this case we set $$s'(a) = f'(a), \quad s'(b)=f'(b).$$
Note that $$f'(a) = s'_0(a) = b_0 \quad \text{and} \quad f'(b) = s'_{n-1}(b)= b_n,$$
so we can use \begin{equation} b_j = \frac{a_{j+1}-a_j}{h_j} - \frac{2c_j+c_{j+1}}{3}h_j \label{eqn3star} \end{equation} with $j=0$ and $j = n-1$ to re-write the boundary conditions in terms of the $c_j$.
I'll let you work out the rest of the algebra on this one!
See the supplementary notes for more details.

Conclusion

To sum up, computing the coefficients defining a cubic spline is a two-step process.
  1. Use your chosen boundary conditions to solve a tridiagonal linear system describing the $c_i$.
  2. Compute $b_i$ and $d_i$ using your $c_i$.
Recall that the $a_i = f(x_i)$ are known!

Results

Evaluating your cubic spline is also a two-step process.
  1. Given $x \in [a, b]$, find $i$ such that $x_i \leq x \leq x_{i+1}$.
  2. Compute $$s_i(x) = \color{var(--emphColor)}{a_i} + \color{var(--emphColor)}{b_i} (x-x_i) + \color{var(--emphColor)}{c_i} (x-x_i)^2 + \color{var(--emphColor)}{d_i}(x-x_i)^3.$$
Congratulations! You reached the end of this lecture.
Take me home