Splines
Contents
Splines¶
It should be clear at this point that if we want to approximate a function with an interpolating polynomial then selecting the nodes as something like Chebyshev points is essential. However, data from many sources comes to us at node intervals that we do not control. Most often, this is at regularly spaced intervals. How should we interpolate between such data? In most cases, the answer to this is that you should use splines. Interpolation using splines is a generalization of “connecting the dots”. The most simple form, a linear interpolating spline, involves just connecting the data points with straight lines to find intermediate values. This is illustrated below using a set of uniformly spaced points with the underlying function \(f(x)=e^{-x^2}\).
This is also a common technique for interpolating between values in a table. An obvious question is: Is this the best we can do with using different polynomials between each set of points (i.e. a function defined piecewise)? Aside from the issue of absolute errors in the interpolation (i.e. \(|f(x)-s_i(x)|\) where \(s_i(x)\) is the polynomial defined on the interval between points \(i\) and \(i+1\)), a linear interpolatory spline is not differentiable at the node points, typically called knots when discussing splines, which could be problematic for many applications. The kinks are also visually unappealing in the plot. It turns out that it is not that difficult to construct an interpolating spline that is continuous in both first and second derivatives at the knots by using piecewise cubic polynomial functions.
Cubic interpolatory splines¶
To be more specific, we are looking for a set of cubic polynomials \(\{ S_j(x) \}_{j=0}^{n-1}\) where our interpolatory function
This is illustrated below:
We want \(S(x)\) to not only be a continuous interpolation function, but continuously differentiable at least two times. This imposes the following four conditions at each interior node:
\(S_j(x_j)=f(x_j),\)
\(S_j(x_{j+1})=S_{j+1}(x_{j+1}),\) which means that \(S_j(x_{j+1})=f(x_{j+1})\) making use of the first condition,
\(S_j'(x_{j+1})=S_{j+1}'(x_{j+1}),\)
\(S_j''(x_{j+1})=S_{j+1}''(x_{j+1}).\)
Each of the \(S_j(x)\) is a cubic polynomial which we will write in the form
Each of the \(n-1\) different \(S_j\) requires four parameters so we will need \(4(n-1)\) conditions to determine them. As the last two conditions above only apply to the interior points, we are short two conditions overall. These are filled in with two boundary conditions. There are a number of possibilities for the boundary conditions. A few common choices are
\(S''(x_0)=0\) and \(S''(x_n)=0.\) This choice is called a natural cubic interpolatory spline.
\(S'(0) = A\) and \(S'(x_n)=B\) with \(A\) and \(B\) specified. This choice is called a clamped cubic interpolatory spline.
\(S_0'''(x_{1})=S_1'''(x_{1})\) and \(S_{n-2}'''(x_{n-1})=S_{n-1}'''(x_{n-1})\). This choice is called not-a-knot as it, together with the continuity of the lower derivatives already required at \(x_1\) and \(x_{n-1}\), forces \(S_0(x)\equiv S_1(x)\) and \(S_{n-2}(x)\equiv S_{n-1}(x)\).
Constructing the interpolatory spline requires using the above conditions to find the parameters in Equation (8) for all \(j\).
Condition 1 can be immediately exploited to find \(a_j\) by noting that evaluating Equation (8) at \(x_j\) gives \(S_j(x_j)=a_j\) and combining this with the first condition gives
We now evaluate \(S_j\) and \(S_{j+1}\) at \(x_{j+1}\) and set them equal, as required by condition 2, to get
where to simplify the notation going forward, we define
We have not yet applied condition 1 and 2 to \(S_{n-1}(x_{n+1})\). However, if we define a “dummy” polynomial \(S_n\) with \(a_n=f(x_n)\) then (10) can be extended to \(j=n-1\).
Before making use of the conditions on the derivatives of \(S(x)\) note that
As as this implies that \(S_{0}'(x)=b_0\) it should be clear that a clamped spline boundary condition explicitly determines \(b_0\) as the (given) derivative at \(x_0\). Similarly, a natural spline boundary condition explicity determines \(c_0=0\).
Continuity of the derivative, condition 3, then implies
To extend this equation to \(j=n-1\) we must make use of a boundary condition at \(x_{j+1}\). If we are constructing a clamped spline this would determine the derivative, and hence \(b_{n}\) for the “dummy” polynomial \(S_n\) if we extend this equation to \(j=n-1\). A natural spline instead just extends this equation to \(j=n-1\) by allowing it to effectively determine the derivative \(S_j'(x_{j+1})\) once we have determined the other parameters on the right hand side.
Similarly, continuity of the second derivative, condition 4, then implies
We can again extend this equation to \(j=n-1\) for a natural spline by replacing \(c_{j+1}\) with zero (the second derivative at \(x_{j+1}\)). The extension for the other boundary conditions are slightly more complicated, but not much. You will construct these in the problems.
Equations (10), (11), and (12) form a coupled linear set of equations for the parameters \(b_j\), \(c_j\), and \(d_j\) given that we know \(a_j\) (from (9)) and \(h_j\). We can simplify these equations further by eliminating \(b_j\) and \(d_j\) to get a set of equations for the \(c_j\) alone (and then once these are known they can be used to deterime the other parameters). We first eliminate the \(d_j\) by rearranging Equation (12),
Substituting this into (10) then gives
which can be solved to give the \(b_j\) in terms of the \(c_j\) to get
Substituting (14) and (13) into (11) gives
This is a set of equations (for different \(j\)) and it is more useful to center it around \(c_j\) than \(c_{j+1}\). i.e. use the equation for \(j\rightarrow j-1\)
\[ h_j c_{j+1} +2(h_j + h_{j-1}) c_j+h_{j-1} c_{j-1} = \frac{3}{h_j}(a_{j+1}-a_j)-\frac{3}{h_{j-1}}(a_j-a_{j-1}), \qquad j=1,2,...n-1\]
and using a boundary condition for \(j=0\). We can rewrite these equation and the boundary contions into a matrix equation as
for the natural cubic spline (the first equation is the boundary condition that sets \(c_0=0\) and the last equation uses the “dummy” polynomial extension of the equations where \(c_n=0\)). This is a diagonally dominant tridiagonal system which means that solving this system via Gaussian elimination/LU factorization without pivoting will provide a unique solution and will be numerically stable. In addition, this can be solved in \(\mathcal{O}(n)\) flops using the Thomas algorithm discussed in the chapter on linear systems.
So to construct our interpolatory spline \(S(x)\) (or, in other words to construct the set \(\{ S_j(x) \}_{j=0}^{n-1}\)) we
Use Eq.(9) to obtain the \(a_j\).
Solve the linear system in Eq.(15) to obtain the \(c_j\).
Use Eq.(14) and Eq.(13), potentially making use of boundary condtions, to obtain the \(b_j\) and \(d_j\).
Once we have all the paramaters, to evaluate \(S(x)\) at a specific point \(x\) we
Determine which interval \([x_j,x_j+1)\) the value \(x\) falls in.
Evaluate and return \(S_j(x)\).
Errors in Cubic Interpolatory Splines¶
As usual, we wish to ensure our errors are bounded and have some form of this error bound so it is clear what measures may, or may not, improve (or worsen) these errors. This is provided by the following theorem which summarizes the results of Hall and Meyer which also provides bounds on the derivatives.
Theorem Given \(f\in C^4[a,b]\), \(\max_{a\leq x \leq b} |f^{(4)}(x)|\leq M\), nodes at
then if \(S(x)\) is the cubic interpolatory spline with clamped end conditions
As we can see, splines can be reasonably good for obtaining approximations to derivatives of the underlying function as well.
Using Splines¶
Implementing splines requires solving a tridiagonal system. While not hard to do, it is usually not worth the effort to implement all of this yourself. Instead, you should use a library. SciPy contains an implementation for natural cubic interpolatory splines which we illustrate the use of below.
We first import the necessary libraries:
from scipy.interpolate import interp1d
import numpy as np
import matplotlib.pyplot as plt
The first call to interp1d creates an instance of the class which sets up the interpolatory spline based on the provided data points. The instance of this class defines a call method which can later be treated like a function which interpolates between known data values to obtain unknown value. Behavior at the boundary can be specified at instantiation time. The default is actually a linear interpolatory spline but contains options to get a cubic.
from scipy.interpolate import interp1d
import numpy as np
x = np.linspace(-2, 2, num=9, endpoint=True)
y = np.exp(-x**2)
f = interp1d(x, y)
f2 = interp1d(x, y, kind='cubic')
Now that we have the splines set up, we can make use of them to, for example, plot out the interpolatory function:
xnew = np.linspace(-2, 2, num=41, endpoint=True)
import matplotlib.pyplot as plt
plt.plot(x, y, 'o', xnew, f(xnew), '-', xnew, f2(xnew), '--')
plt.legend(['data', 'linear', 'cubic'], loc='best')
plt.show()
