Cubic splines - Fundamentals of Numerical Computation

A piecewise linear interpolant is continuous but has discontinuities in its derivative. We often desire a smoother interpolant, i.e., one that has some continuous derivatives. By far the most popular choice is piecewise cubic.

We use $S(x)$ to denote the cubic spline interpolant. As before, suppose that distinct nodes $t_0 < t_1 < \cdots < t_n$ (not necessarily equally spaced) and data $y_0,\ldots,y_n$ are given. For any $k=1,\ldots,n$ , the spline $S(x)$ on the interval $[t_{k-1},t_k]$ is by definition a cubic polynomial $S_k(x)$ , which we express as

S_k(x) = a_k + b_k(x-t_{k-1}) + c_k(x-t_{k-1})^2 + d_k(x-t_{k-1})^3, \qquad k=1,\ldots,n,

(5.3.1)

where $a_k,b_k,c_k,d_k$ are values to be determined. Overall there are $4n$ such undetermined coefficients.

5.3.1Smoothness conditions¶

We are able to ensure that $S$ has at least two continuous derivatives everywhere by means of the following constraints.

1. Interpolation by $S_k$ at both of its endpoints.

Algebraically we require $S_k(t_{k-1})=y_{k-1}$ and $S_k(t_k)=y_k$ for every $k=1,\dots,n$ . In terms of (5.3.1), these conditions are

a_k = y_{k-1},

(5.3.2)

a_k + b_k h_k + c_k h_k^2 + d_k h_k^3 = y_{k}, \qquad k=1,\ldots,n,

(5.3.3)

where we have used the definition

h_k = t_{k}-t_{k-1}, \qquad k=1,\ldots,n.

(5.3.4)

The values of $h_k$ are derived from the nodes. Crucially, the unknown coefficients appear only linearly in the constraint equations. So we will express the constraints using linear algebra. The left endpoint interpolation constraints (5.3.2) are, in matrix form,

\begin{bmatrix} \mathbf{I} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} \end{bmatrix} \begin{bmatrix} \mathbf{a} \\ \mathbf{b} \\ \mathbf{c} \\ \mathbf{d} \end{bmatrix} = \begin{bmatrix} y_0 \\ \vdots \\ y_{n-1} \end{bmatrix},

(5.3.5)

with $\mathbf{I}$ being an $n\times n$ identity. The right endpoint interpolation constraints, given by (5.3.3), become

\begin{bmatrix} \mathbf{I} & \mathbf{H} & \mathbf{H}^2 & \mathbf{H}^3 \end{bmatrix} \begin{bmatrix} \mathbf{a} \\ \mathbf{b} \\ \mathbf{c} \\ \mathbf{d} \end{bmatrix} = \begin{bmatrix} y_1 \\ \vdots \\ y_{n} \end{bmatrix},

(5.3.6)

where we have defined the diagonal matrix

\mathbf{H} = \begin{bmatrix} h_1 & & & \\ & h_2 & & \\ & & \ddots & \\ & & & h_n \end{bmatrix}.

(5.3.7)

Collectively, (5.3.5) and (5.3.6) express $2n$ scalar constraints on the unknowns.

2. Continuity of $S'(x)$ at interior nodes.

We do not know what the slope of the interpolant should be at the nodes, but we do want the same slope whether a node is approached from the left or the right. Thus we obtain constraints at the nodes that sit between two neighboring piecewise definitions, so that $S_1'(t_1)=S_2'(t_1)$ , and so on. Altogether these are

b_k + 2 c_k h_k + 3 d_k h_k^2 = b_{k+1}, \qquad k=1,\dots,n-1.

(5.3.8)

Moving the unknowns to the left side, as a system these become

\mathbf{E} \begin{bmatrix} \boldsymbol{0} & \mathbf{J} & 2\mathbf{H} & 3\mathbf{H}^2 \end{bmatrix} \begin{bmatrix} \mathbf{a} \\ \mathbf{b} \\ \mathbf{c} \\ \mathbf{d} \end{bmatrix} = \boldsymbol{0},

(5.3.9)

where now we have defined

\mathbf{J} = \begin{bmatrix} 1 & -1 & & & \\ & 1 & -1 & & \\ & & \ddots & \ddots & \\ & & &1 & -1 \\ & & & & 1 \end{bmatrix},

(5.3.10)

and $\mathbf{E}$ is the $(n-1)\times n$ matrix resulting from deleting the last row of the identity:

\mathbf{E} = \begin{bmatrix} 1 & 0 & & & \\ & 1 & 0 & & \\ & & \ddots & \ddots & \\ & & & 1& 0 \end{bmatrix}.

(5.3.11)

Left-multiplying by $\mathbf{E}$ deletes the last row of any matrix or vector. Hence (5.3.9) represents $n-1$ constraints on the unknowns. (Remember, there are only $n-1$ interior nodes.)

3. Continuity of $S''(x)$ at interior nodes.

These again apply only at the interior nodes $t_1,\dots,t_{n-1}$ , in the form $S_1''(t_1)=S_2''(t_1)$ and so on. Using (5.3.1) once more, we obtain

2 c_k + 6 d_k h_k = 2c_{k+1}, \qquad k=1,\dots,n-1.

(5.3.12)

In system form (after canceling a factor of 2 from each side) we get

\mathbf{E} \begin{bmatrix} \boldsymbol{0} & \boldsymbol{0} & \mathbf{J} & 3\mathbf{H} \end{bmatrix} \begin{bmatrix} \mathbf{a} \\ \mathbf{b} \\ \mathbf{c} \\ \mathbf{d} \end{bmatrix} = \boldsymbol{0}.

(5.3.13)

5.3.2End conditions¶

So far the equations (5.3.5), (5.3.6), (5.3.9), and (5.3.13) form $2n+(n-1)+(n-1)=4n-2$ linear conditions on the $4n$ unknowns in the piecewise definition (5.3.1). In order to obtain a square system, we must add two more constraints. If the application prescribes values for $S'$ or $S''$ at the endpoints, those may be applied. Otherwise there are two major alternatives:

Natural spline: $\quad S_1''(t_0)=S_n''(t_n)=0$
Not-a-knot spline: $\quad S_1'''(t_1)=S_2'''(t_1), \; S_{n-1}'''(t_{n-1})=S_n'''(t_{n-1})$

While natural splines have important theoretical properties, not-a-knot splines give better pointwise accuracy, and they are the only type we consider further.

In the not-a-knot spline, the values and first three derivatives of the cubic polynomials $S_1$ and $S_2$ agree at the node $t_1$ . Hence they must be the same cubic polynomial! The same is true of $S_{n-1}$ and $S_n$ .^[1] We could use these facts to eliminate some of the undetermined coefficients from our linear system of constraints. However, rather than rework the algebra we just append two more rows to the system, expressing the conditions

d_1=d_2, \quad d_{n-1}=d_n.

(5.3.14)

Collectively, (5.3.5), (5.3.6), (5.3.9), (5.3.13), and (5.3.14) comprise a square linear system of size $4n$ which can be solved for the coefficients defining the piecewise cubics in (5.3.1). This is a major difference from the piecewise linear interpolant, for which there is no linear system to solve. Indeed, while it is possible to find a basis for the cubic spline interpolant analogous to the hat functions, it is not possible in closed form to construct a cardinal basis, so the solution of a linear system cannot be avoided.

5.3.3Implementation¶

Algorithm 5.3.1 (spinterp)

Julia

MATLAB

Python

Cubic spline interpolation

spinterp.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
function S = spinterp(t,y)
% SPINTERP   Cubic not-a-knot spline interpolation.
% Input:
%   t     interpolation nodes (vector, length n+1)
%   y     interpolation values (vector, length n+1)
% Output:
%   S     not-a-knot cubic spline (function)

t = t(:);  y = y(:);  % ensure column vectors
n = length(t)-1;
h = diff(t);          % differences of all adjacent pairs

% Preliminary definitions.
Z = zeros(n);
I = eye(n);  E = I(1:n-1,:);
J = I - diag(ones(n-1,1),1);
H = diag(h);

% Left endpoint interpolation:
AL = [ I, Z, Z, Z ];
vL = y(1:n);

% Right endpoint interpolation:
AR = [ I, H, H^2, H^3 ];
vR = y(2:n+1);

% Continuity of first derivative:
A1 = E*[ Z, J, 2*H, 3*H^2 ];
v1 = zeros(n-1,1);

% Continuity of second derivative:
A2 = E*[ Z, Z, J, 3*H ];
v2 = zeros(n-1,1);

% Not-a-knot conditions:
nakL = [ zeros(1,3*n), [1,-1, zeros(1,n-2)] ];
nakR = [ zeros(1,3*n), [zeros(1,n-2), 1,-1] ];

% Assemble and solve the full system.
A = [ AL; AR; A1; A2; nakL; nakR ];
v = [ vL; vR; v1; v2; 0 ;0 ];
z = A\v;

% Break the coefficients into separate vectors.
rows = 1:n;
a = z(rows);
b = z(n+rows);  c = z(2*n+rows);  d = z(3*n+rows);
S = @evaulate;

    % This function evaluates the spline when called with a value for x.
    function f = evaulate(x)
        f = zeros(size(x));
        for k = 1:n       % iterate over the pieces
            % Evalaute this piece's cubic at the points inside it.
            index = (x>=t(k)) & (x<=t(k+1));   
            f(index) = polyval( [d(k),c(k),b(k),a(k)], x(index)-t(k) );
        end
    end

end

Cubic spline interpolation

spinterp.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
def spinterp(t, y):
    """
    spinterp(t,y)

    Create a cubic not-a-knot spline interpolating function for data values in `y` given at nodes in `t`.
    """
    n = len(t) - 1
    h = [t[i + 1] - t[i] for i in range(n)]

    # Preliminary definitions.
    Z = np.zeros([n, n])
    I = np.eye(n)
    E = I[: n - 1, :]
    J = np.eye(n) + np.diag(-np.ones(n - 1), 1)
    H = np.diag(h)

    # Left endpoint interpolation:
    AL = np.hstack([I, Z, Z, Z])
    vL = y[:-1]

    # Right endpoint interpolation:
    AR = np.hstack([I, H, H**2, H**3])
    vR = y[1:]

    # Continuity of first derivative:
    A1 = E @ np.hstack([Z, J, 2 * H, 3 * H**2])
    v1 = np.zeros(n - 1)

    # Continuity of second derivative:
    A2 = E @ np.hstack([Z, Z, J, 3 * H])
    v2 = np.zeros(n - 1)

    # Not-a-knot conditions:
    nakL = np.hstack([np.zeros(3 * n), np.hstack([1, -1, np.zeros(n - 2)])])
    nakR = np.hstack([np.zeros(3 * n), np.hstack([np.zeros(n - 2), 1, -1])])

    # Assemble and solve the full system.
    A = np.vstack([AL, AR, A1, A2, nakL, nakR])
    v = np.hstack([vL, vR, v1, v2, 0, 0])
    z = solve(A, v)

    # Break the coefficients into separate vectors.
    rows = np.arange(n)
    a = z[rows]
    b = z[n + rows]
    c = z[2 * n + rows]
    d = z[3 * n + rows]
    S = [np.poly1d([d[k], c[k], b[k], a[k]]) for k in range(n)]

    # This function evaluates the spline when called with a value for x.
    def evaluate(x):
        f = np.zeros(x.shape)
        for k in range(n):
            # Evaluate this piece's cubic at the points inside it.
            index = (x >= t[k]) & (x <= t[k + 1])
            f[index] = S[k](x[index] - t[k])
        return f

    return evaluate

Function 5.3.1 gives an implementation of cubic not-a-knot spline interpolation. For clarity it stays very close to the description given above. There are some possible shortcuts—for example, one could avoid using $\mathbf{E}$ and instead directly delete the last row of any matrix it left-multiplies. Observe that the linear system is assembled and solved just once, and the returned evaluation function simply uses the resulting coefficients. This allows us to make multiple calls to evaluate $S$ without unnecessarily repeating the linear algebra.

5.3.4Conditioning and convergence¶

Example 5.3.1

Julia

MATLAB

Python

Cubic splines

For illustration, here is a spline interpolant using just a few nodes.

f = x -> exp(sin(7 * x))

plot(f, 0, 1, label="function", xlabel=L"x", ylabel=L"y")

t = [0, 0.075, 0.25, 0.55, 0.7, 1]  # nodes
y = f.(t)                           # values at nodes

scatter!(t, y, label="values at nodes")

S = FNC.spinterp(t, y)

plot!(S, 0, 1, label="spline")

Now we look at the convergence rate as the number of nodes increases.

x = (0:10000) / 1e4              # sample the difference at many points
n = @. round(Int, 2^(3:0.5:7))  # numbers of nodes
err = zeros(0)
for n in n
    t = (0:n) / n
    S = FNC.spinterp(t, f.(t))
    dif = @. f(x) - S(x)
    push!(err, norm(dif, Inf))
end

pretty_table((; n, err), header=["n", "error"])

Since we expect convergence that is $O(h^4)=O(n^{-4})$ , we use a log-log graph of error and expect a straight line of slope -4.

order4 = @. (n / n[1])^(-4)

plot(n, [err order4], m=[:o :none], l=[:solid :dash],
    label=["error" "4th order"],
    xaxis=(:log10, "n"), yaxis=(:log10, L"|| f-S\,||_\infty"),
    title="Convergence of spline interpolation")

Cubic splines

For illustration, here is a spline interpolant using just a few nodes.

clf
f = @(x) exp(sin(7 * x));
fplot(f, [0, 1], displayname="function")

t = [0, 0.075, 0.25, 0.55, 0.7, 1];    % nodes
y = f(t);                              % values at nodes
hold on, scatter(t, y, displayname="values at nodes")

S = spinterp(t, y);
fplot(S, [0, 1], displayname="spline")

xlabel("x");  ylabel("y")
legend()

Now we look at the convergence rate as the number of nodes increases.

x = (0:10000)' / 1e4;              % sample the difference at many points
n = round(2 .^ (3:0.5:7))';        % numbers of nodes
maxerr = zeros(size(n));
for i = 1:length(n)
    t = (0:n(i))' / n(i);
    S = spinterp(t, f(t));
    err = f(x) - S(x);
    maxerr(i) = norm(err, Inf);
end
disp(table(n(1:2:end), maxerr(1:2:end), variableNames=["n", "inf-norm error"]))

Since we expect convergence that is $O(h^4)=O(n^{-4})$ , we use a log-log graph of error and expect a straight line of slope -4.

clf
loglog(n, maxerr, "-o", displayname="error")
order4 = 0.5 * maxerr(end) * (n / n(end)) .^ (-4);
hold on
loglog(n, order4, "k--", displayname="O(n^{-4})")
xlabel("n");  ylabel("|| f-S ||_{\infty}")
title("Convergence of spline interpolation")

Cubic splines

For illustration, here is a spline interpolant using just a few nodes.

f = lambda x: exp(sin(7 * x))

x = linspace(0, 1, 500)
fig, ax = subplots()
ax.plot(x, f(x), label="function")

t = array([0, 0.075, 0.25, 0.55, 0.7, 1])  # nodes
y = f(t)  # values at nodes

xlabel("$x$")
ylabel("$y$")
ax.scatter(t, y, label="nodes")

S = FNC.spinterp(t, y)
ax.plot(x, S(x), label="spline")
ax.legend()
fig

Now we look at the convergence rate as the number of nodes increases.

N = floor(2 ** linspace(3, 8, 17)).astype(int)
err = zeros(N.size)
for i, n in enumerate(N):
    t = linspace(0, 1, n + 1)  # interpolation nodes
    p = FNC.spinterp(t, f(t))
    err[i] = max(abs(f(x) - p(x)))
print(err)

Since we expect convergence that is $O(h^4)=O(n^{-4})$ , we use a log-log graph of error and expect a straight line of slope -4.

order4 = (N / N[0]) ** (-4)
loglog(N, err, "-o", label="observed error")
loglog(N, order4, "--", label="4th order")
xlabel("$n$")
ylabel("$\|f-S\|_\infty$")
legend()

Besides having more smoothness than a piecewise linear interpolant, the not-a-knot cubic spline improves the order of accuracy to 4.

The conditioning of spline interpolation is much more complicated than for the piecewise linear case. First, the fact that the coefficients of all the cubics must be solved for simultaneously implies that each data value in $\mathbf{y}$ has an influence on $S$ over the entire interval. Second, $S$ can take on values larger in magnitude than all of the values in $\mathbf{y}$ (see Exercise 5). The details may be found in more advanced texts.

5.3.5Exercises¶

✍ In each case, write out the entries of the matrix and right-hand side of the linear system that determines the coefficients for the cubic not-a-knot spline interpolant of the given function and node vector.
(a) $\cos (\pi^2 x^2 ), \: \mathbf{t} = [-1,1,4]$
(b) $\cos (\pi^2 x^2), \: \mathbf{t} = [0,\tfrac{1}{2},\tfrac{3}{4},1]$
(c) $\ln(x), \: \mathbf{t} = [1,2,3]$
(d) $\sin(x^2),\: \mathbf{t} = [-1,0,1]$
only - Unknown Directive
⌨ (continuation) For each case in the preceding problem, use Julia to solve the linear system you wrote down. Then plot the resulting cubic spline over the interval between the second and third nodes.
⌨ For each given function, interval, and value of $n$ , define $n+1$ evenly spaced nodes. Then use Function 5.3.1 to plot the cubic spline interpolant using those nodes, together with the original function over the given interval.
(a) $\cos(\pi x^2)$ on $[0,4]$ , $n=18$
(b) $\ln(x)$ on $[1,20]$ , $n=4$
(c) $\sin\left(\frac{1}{x}\right)$ on $\left[\frac{1}{2},7\right]$ , $n=9$
⌨ For each given function and interval, perform piecewise linear interpolation using Function 5.3.1 for $n+1$ equispaced nodes with $n=10,20,40,80,160,320$ . For each $n$ , estimate the error
$E(n) = \| f-p \|_\infty = \max_x | f(x) - p(x) |$
(5.3.16)
by evaluating the function and interpolant at 1600 points in the interval. Make a log-log plot of $E$ as a function of $n$ and add the line $E=Cn^{-4}$ for a constant $C$ of your choosing.
(a) $\cos(\pi x^2)$ on $[0,4]$
(b) $\ln(x)$ on $[1,20]$
(c) $\sin\left(\frac{1}{x}\right)$ on $\left[\frac{1}{2},7\right]$
⌨ Although the cardinal cubic splines are intractable in closed form, they can be found numerically. Each cardinal spline interpolates the data from one column of an identity matrix. Define the nodes $\mathbf{t} = \bigl[0,\, 0.075,\, 0.25,\, 0.55,\, 1]$ . Plot over $[0,1]$ the five cardinal functions for this node set over the interval $[0,1]$ .
✍ Suppose you were to define a piecewise quadratic spline that interpolates $n+1$ given values and has a continuous first derivative. Follow the derivation of this section to express all of the interpolation and continuity conditions. How many additional conditions are required to make a square system for the coefficients?
(a) ✍ If $y_0=y_n$ , another possibility for cubic spline end conditions is to make $S(x)$ a periodic function. This implies that $S'$ and $S''$ are also periodic. Write out the two new algebraic equations for these constraints in terms of the piecewise coefficients.
(b) ⌨ Modify Function 5.3.1 to compute a periodic spline interpolant. Test by making a plot of the interpolant for $f(x) =\exp(\sin(3x))$ over the interval $[0,2\pi/3]$ with equally spaced nodes and $n=8$ .

Footnotes¶

This explains the name of the not-a-knot spline—for splines, “knots” are the points at which different piecewise definitions meet.
↩

Preface

Piecewise linear interpolation

Preface

Finite differences