28-2 Splines

A pratical method for interpolating a set of points with a curve is to use cubic splines. We are given a set $\{(x_i, y_i): i = 0, 1, \ldots, n\}$ of $n + 1$ point-value pairs, where $x_0 < x_1 < \cdots < x_n$. We wish to fit a piecewise-cubic curve (spline) $f(x)$ to the points. That is, the curve $f(x)$ is made up of $n$ cubic polynomials $f_i(x) = a_i + b_ix + c_ix^2 + d_ix^3$ for $i = 0, 1, \ldots, n - 1$, where if $x$ falls in the range $x_i \le x \le x_{i + 1}$, then the value of the curve is given by $f(x) = f_i(x - x_i)$. The points $x_i$ at which the cubic polynomials are "pasted" together are called knots. For simplicity, we shall assume that $x_i = i$ for $i = 0, 1, \ldots, n$.

To ensure continuity of $f(x)$, we require that

$$ \begin{aligned} f(x_i) & = f_i(0) = y_i, \\ f(x_{i + 1}) & = f_i(1) = y_{i + 1} \end{aligned} $$

for $i = 0, 1, \ldots, n - 1$. To ensure that $f(x)$ is sufficiently smooth, we also insist that the first derivative be continuous at each knot:

$$f'(x_{i + 1}) = f'_i(1) = f'_{i + 1}(0)$$

for $i = 0, 1, \ldots, n - 2$.

a. Suppose that for $i = 0, 1, \ldots, n$, we are given not only the point-value pairs $\{(x_i, y_i)\}$ but also the first derivatives $D_i = f'(x_i)$ at each knot. Express each coefficient $a_i$, $b_i$, $c_i$ and $d_i$ in terms of the values $y_i$, $y_{i + 1}$, $D_i$, and $D_{i + 1}$. (Remember that $x_i = i$.) How quickly can we compute the $4n$ coefficients from the point-value pairs and first derivatives?

The question remains of how to choose the first derivatives of $f(x)$ at the knots. One method is to require the second derivatives to be continuous at the knots:

$$f''(x_{i + 1}) = f''_i(1) = f''_{i + 1}(0)$$

for $i = 0, 1, \ldots, n - 2$. At the first and last knots, we assume that $f''(x_0) = f''_0(0) = 0$ and $f''(x_n) = f''_{n - 1}(1) = 0$; these assumptions make $f(x)$ a natural cubic spline.

b. Use the continuity constraints on the second derivative to show that for $i = 1, 2, \ldots, n - 1$,

$$D_{i - 1} + 4D_i + D_{i + 1} = 3(y_{i + 1} - y_{i - 1}). \tag{23.21}$$

c. Show that

$$ \begin{aligned} 2D_0 + D_1 & = 3(y_1 - y_0), & \text{(28.22)} \\ D_{n - 1} + 2D_n & = 3(y_n - y_{n - 1}). & \text{(28.23)} \end{aligned} $$

d. Rewrite equations $\text{(28.21)}$–$\text{(28.23)}$ as a matrix equation involving the vector $D = \langle D_0, D_1, \ldots, D_n \rangle$ or unknowns. What attributes does the matrix in your equation have?

e. Argue that a natural cubic spline can interpolate a set of $n + 1$ point-value pairs in $O(n)$ time (see Problem 28-1).

f. Show how to determine a natural cubic spline that interpolates a set of $n + 1$ points $(x_i, y_i)$ satisfying $x_0 < x_1 < \cdots < x_n$, even when $x_i$ is not necessarily equal to $i$. What matrix equation must your method solve, and how quickly does your algorithm run?

a. We have $a_i = f_i(0) = y_i$ and $b_i = f_i'(0) = f'(x_i) = D_i$. Since $f_i(1) = a_i + b_i + c_i + d_i$ and $f_i'(1) = b_i + 2c_i + 3d_i$, we have $d_i = D_{i + 1} - 2y_{i + 1} + 2y_i + D_i$ which implies $c_i = 3y_{i + 1} - 3y_i - D_{i + 1} - 2D_i$. Since each coefficient can be computed in constant time from the known values, we can compute the $4n$ coefficients in linear time.

b. By the continuity constraints, we have $f_i''(1) = f_{i + 1}''(0)$ which implies that $2c_i + 6d_i = 2c_{i + 1}$, or $c_i + 3d_i = c_{i + 1}$. Using our equations from above, this is equivalent to

$$D_i + 2D_{i + 1} + 3y_i - 3y_{i + 1} = 3y_{i + 2} - 3y_{i + 1} - D_{i + 2} - 2D_{i + 1}.$$

Rearranging gives the desired equation

$$D_i + 4D_{i + 1} + D_{i + 2} = 3(y_{i + 2} - y_i).$$

c. The condition on the left endpoint tells us that $f_0''(0) = 0$, which implies $2c_0 = 0$. By part (a), this means $3(y_1 − y_0) = 2D_0 + D_1$. The condition on the right endpoint tells us that $f_{n - 1}''(1) = 0$, which implies $c_{n - 1} + 3d_{n - 1} = 0$. By part (a), this means $3(y_n - y_{n - 1}) = D_{n - 1} + 2D_n$.

d. The matrix equation has the form $AD = Y$, where $A$ is symmetric and tridiagonal. It looks like this:

$$ \begin{pmatrix} 2 & 1 & 0 & 0 & \cdots & 0 \\ 1 & 4 & 1 & 0 & \cdots & 0 \\ 0 & \ddots & \ddots & \ddots & \cdots & \vdots \\ \vdots & \cdots & 1 & 4 & 1 & 0 \\ 0 & \cdots & 0 & 1 & 4 & 1 \\ 0 & \cdots & 0 & 0 & 1 & 2 \\ \end{pmatrix} \begin{pmatrix} D_0 \\ D_1 \\ D_2 \\ \vdots \\ D_{n - 1} \\ D_n \end{pmatrix} = \begin{pmatrix} 3(y_1 - y_0) \\ 3(y_2 - y_0) \\ 3(y_3 - y_1) \\ \vdots \\ 3(y_n - y_{n - 2}) \\ 3(y_n - y_{n - 1}) \end{pmatrix} . $$

e. Since the matrix is symmetric and tridiagonal, Problem 28-1 (e) tells us that we can solve the equation in $O(n)$ time by performing an LUP decomposition. By part (a), once we know each $D_i$ we can compute each $f_i$ in $O(n)$ time.

f. For the general case of solving the nonuniform natural cubic spline problem, we require that $f(x_{i + 1}) = f_i(x_{i + 1} − x_i) = y_{i + 1}$, $f'(x_{i + 1}) = f_i'(x_{i + 1} - x_i) = f_{i + 1}'(0)$ and $f''(x_{i + 1}) = f_i''(x_{i + 1} - x_i) = f_{i + 1}''(0)$. We can still solve for each of $a_i$, $b_i$, $c_i$ and $d_i$ in terms of $y_i$, $y_{i + 1}$, $D_i$ and $D_{i + 1}$, so we still get a tridiagonal matrix equation. The solution will be slightly messier, but ultimately it is solved just like the simpler case, in $O(n)$ time.