from numpy import *
from scipy import linalg
from scipy.linalg import norm
from matplotlib.pyplot import *
from prettytable import PrettyTable
from timeit import default_timer as timer
import sys
sys.path.append('fncbook/')
import fncbook as FNC

# This (optional) block is for improving the display of plots.
# from IPython.display import set_matplotlib_formats
# set_matplotlib_formats("svg","pdf")
# %config InlineBackend.figure_format = 'svg'
rcParams["figure.figsize"] = [7, 4]
rcParams["lines.linewidth"] = 2
rcParams["lines.markersize"] = 4
rcParams['animation.html'] = "jshtml"  # or try "html5"

5.3.Cubic splines¶

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

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

where we have used the definition

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,

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

where we have defined the diagonal matrix

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

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

where now we have defined

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

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

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

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

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¶

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¶

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 the values in $\mathbf{y}$ (see Exercise 5.3.5). The details may be found in more advanced texts.

5.3.5Exercises¶