7.1 Cubic spline interpolation
The simplest piecewise interpolant joins consecutive data points by straight lines. It is continuous, but has corners at the nodes. Cubic pieces have enough freedom to match values, slopes and curvatures across the nodes.
Definition 7.1.1. Let \(a = x_0<x_1<\dots <x_n = b\). A cubic spline interpolant \(S\) for \(f\) is a function that is a cubic polynomial \(S_j\) on each \([x_j,x_{j+1}]\), \(j = 0,\dots ,n-1\), such that
- (a)
- \(S_j(x_j) = f(x_j)\) and \(S_j(x_{j+1}) = f(x_{j+1})\) for each \(j\);
- (b)
- \(S'_{j+1}(x_{j+1}) = S'_j(x_{j+1})\) and \(S''_{j+1}(x_{j+1}) = S''_j(x_{j+1})\) for \(j = 0,\dots ,n-2\);
together with one of the boundary conditions
- (c)
- \(S''(x_0) = S''(x_n) = 0\) (a natural spline), or
- (d)
- \(S'(x_0) = f'(x_0)\) and \(S'(x_n) = f'(x_n)\) (a clamped spline).
So \(S\) interpolates \(f\) and has continuous first and second derivatives. The natural condition models a flexible strip of wood (the draughtsman’s “spline”) pinned at the data points and free at its ends; the clamped condition uses extra information about \(f\), and is more accurate when that information is available.
Constructing the spline
Write \(S_j(x) = a_j+b_j(x-x_j)+c_j(x-x_j)^2+d_j(x-x_j)^3\) and \(h_j = x_{j+1}-x_j\). Condition (a) gives \(a_j = f(x_j)\). Since \(S_j''(x_j) = 2c_j\), continuity of \(S''\) gives \(c_{j+1} = c_j+3d_jh_j\), so \[d_j = \frac {c_{j+1}-c_j}{3h_j} .\] Substituting into \(S_j(x_{j+1}) = a_{j+1}\) and solving for \(b_j\), \[b_j = \frac {a_{j+1}-a_j}{h_j}-\frac {h_j}{3}(2c_j+c_{j+1}),\] and then continuity of \(S'\), \(b_j = b_{j-1}+2c_{j-1}h_{j-1}+3d_{j-1}h_{j-1}^2\), becomes the linear system \begin {equation} \label {eq:splinesys} h_{j-1}c_{j-1}+2(h_{j-1}+h_j)c_j+h_jc_{j+1} = \frac {3}{h_j}(a_{j+1}-a_j)-\frac {3}{h_{j-1}}(a_j-a_{j-1}),\qquad j = 1,\dots ,n-1 . \end {equation} For a natural spline, \(c_0 = c_n = 0\); for a clamped spline the two boundary conditions give \[2h_0c_0+h_0c_1 = \frac {3}{h_0}(a_1-a_0)-3f'(a),\qquad h_{n-1}c_{n-1}+2h_{n-1}c_n = 3f'(b)-\frac {3}{h_{n-1}}(a_n-a_{n-1}) .\] Once the \(c_j\) are known, the \(b_j\) and \(d_j\) follow.
Theorem 7.1.2. A function defined at \(a = x_0<\dots <x_n = b\) has a unique natural spline interpolant and, if \(f'(a)\) and \(f'(b)\) are given, a unique clamped spline interpolant.
Proof. In both cases the \(c_j\) satisfy a tridiagonal system whose matrix is strictly diagonally dominant: in equation (7.1) the diagonal entry is \(2(h_{j-1}+h_j)>h_{j-1}+h_j\), the sum of the off-diagonal entries, and the same holds for the boundary rows (\(1>0\) in the natural case, \(2h_0>h_0\) in the clamped case). By Theorem 4.5.2 the system has a unique solution, which determines the spline. It can be computed by the tridiagonal Algorithm 4.5.6 in \(O(n)\) operations. □
Solution. Here \(h_0 = h_1 = 1\), \(c_0 = c_2 = 0\), and (7.1) with \(j = 1\) is \[4c_1 = 3(5-3)-3(3-2) = 3,\qquad c_1 = \tfrac 34 .\] Then \(b_0 = (3-2)-\frac 13\left (0+\frac 34\right ) = \frac 34\), \(b_1 = (5-3)-\frac 13\left (\frac 32+0\right ) = \frac 32\), \(d_0 = \frac {3/4}{3} = \frac 14\), \(d_1 = \frac {-3/4}{3} = -\frac 14\). So \[S(x) = \begin {cases}2+\frac 34(x-1)+\frac 14(x-1)^3, & 1\leq x\leq 2,\\[2pt] 3+\frac 32(x-2)+\frac 34(x-2)^2-\frac 14(x-2)^3, & 2\leq x\leq 3 .\end {cases}\] One checks that both pieces have value \(3\), slope \(\frac 32\) and second derivative \(\frac 32\) at \(x = 2\).
Example 7.1.4. Construct the natural and the clamped cubic splines for \(f(x) = e^x\) at \(x = 0,1,2,3\), and use each to approximate \(\int _0^3e^x\,dx = e^3-1 = 19.08554\).
Solution. With \(h_j = 1\), the natural spline system is \(4c_1+c_2 = 3(e^2-2e+1)\), \(c_1+4c_2 = 3(e^3-2e^2+e)\), and the coefficients are
| \(j\) | \(a_j\) | \(b_j\) | \(c_j\) | \(d_j\) |
| \(0\) | \(1.00000\) | \(1.46600\) | \(0\) | \(0.25228\) |
| \(1\) | \(2.71828\) | \(2.22285\) | \(0.75685\) | \(1.69107\) |
| \(2\) | \(7.38906\) | \(8.80977\) | \(5.83007\) | \(-1.94336\) |
For the clamped spline, with \(f'(0) = 1\) and \(f'(3) = e^3\):
| \(j\) | \(a_j\) | \(b_j\) | \(c_j\) | \(d_j\) |
| \(0\) | \(1.00000\) | \(1.00000\) | \(0.44468\) | \(0.27360\) |
| \(1\) | \(2.71828\) | \(2.71016\) | \(1.26548\) | \(0.69513\) |
| \(2\) | \(7.38906\) | \(7.32652\) | \(3.35087\) | \(2.01909\) |
Integrating each piece, \(\int _{x_j}^{x_{j+1}}S_j = a_jh_j+\frac {b_j}2h_j^2+\frac {c_j}3h_j^3+\frac {d_j}4h_j^4\), gives \(19.55229\) for the natural spline and \(19.05965\) for the clamped one: errors \(0.47\) and \(0.026\). The clamped spline is far better because it uses the true slopes at the ends, while the natural condition \(S'' = 0\) is false for \(e^x\). The same is seen pointwise: at \(x = 2.5\) the natural spline gives \(13.009\), the clamped \(12.142\), and \(e^{2.5} = 12.182\).
Note. For the clamped spline there is the error bound \[\max _{a\leq x\leq b}|f(x)-S(x)|\leq \frac {5M}{384}\max _jh_j^4,\qquad M = \max _{a\leq x\leq b}|f^{(4)}(x)|,\] so halving the spacing reduces the error by a factor of about \(16\); the natural spline has the same \(O(h^4)\) accuracy in the interior but only \(O(h^2)\) near the ends unless \(f''(a) = f''(b) = 0\). Splines do not suffer from Runge’s phenomenon: refining the nodes always improves the fit.
Questions on this section
Stuck on something here? Ask below and it stays attached to this topic.