Lecture 13: Interpolation and splines
Jamie Haddock
In several of our recent tasks, and in many scientific problems, the solution is a function. Representing functions nuerically is a more complicated and difficult task than that of representing real numbers. The process of converting functions into numberical represesentations of finite length is known as discretization.
The interpolation problem
Definition: Interpolation problem
Given \(n+1\) distinct points \((t_0, y_0), (t_1,y_1), \cdots, (t_n,y_n)\), with \(t_0 < t_1 < \cdots < t_n\) called nodes, the interpolation problem is to find a function \(p(x)\), called the interpolant, such that \(p(t_k) = y_k\) for \(k = 0, \cdots, n\).
Polynomials
Polynomials are a natural first candidate for an interpolant, but it is easy to find examples of data where a polynomial interpolant is unusable for the application or task.
The points are unremarkable but the interpolant is quite unusual.
Fact:
Interpolation by a polynomial at equally spaced nodes is ill-conditioned as the degree of the polynomial grows.
Piecewise polynomials
To keep polynomial degrees low while interpolating large data sets, we use piecewise polynomials as interpolants – that is \(p\) is a polynomial on each subinterval \([t_{k-1},t_k]\) for \(k = 1, \cdots, n\).
We can even choose a higher degree polynomial for a smoother interpolant.
Conditioning of interpolation
We consider the nodes to be fixed and let \(a = t_0, b = t_n\), so the data for the interpolant problem consists of the vector of \(y\) values, \(\mathbf{y}\), and the result is a function on \([a,b]\).
Suppose \(\mathcal{I}\) is a method for producing the interpolant from a data vector – that is \(\mathcal{I}(\mathbf{y}) = p\), where \(p(t_k) = y_k\) for all \(k\).
We consider methods that are linear, which in particular, satisfy \[\mathcal{I}(\mathbf{y}) = \mathcal{I}\left(\sum_{k=0}^n y_k \mathbf{e}_k\right) = \sum_{k=0}^n y_k \mathcal{I}(\mathbf{e}_k).\]
Definition: Cardinal function
A cardinal function \(\phi_k\) for a node set \(t_0, \cdots, t_n\) is the function that interpolates the value \((t_k,1)\) and \((t_j,0)\) for all \(j \not= k\)
For any set of \(n+1\) nodes, there are \(n+1\) cardinal functions \(\phi_0, \cdots, \phi_n\) which each single out a different interpolation node, and we then have \[\mathcal{I}(\mathbf{y}) = \sum_{k=0}^n y_k \phi_k.\]
We’ll now measure the conditioning of the interpolation problem with respect to the function \(\infty\)-norm defined as \[\|f\|_\infty = \max_{x \in [a,b]} |f(x)|.\]
Theorem: Conditioning of interpolation
Suppose that \(\mathcal{I}\) is a linear interpolation method on nodes \(t_0, \cdots, t_n\). Then with respect to \(\|\cdot\|_\infty\), the absolute condition number of \(\mathcal{I}\) satisfies \[\max_{0\le k \le n} \|\phi_k\|_\infty \le \kappa(\mathbf{y}) \le \sum_{k=0}^n \|\phi_k\|_\infty,\] where the \(\phi_k\) are cardinal interpolating functions.
Proof:
Suppose the data is perturbed from \(\mathbf{y}\) to \(\mathbf{y} + \mathbf{d}\), then \[\frac{\|\mathcal{I}(\mathbf{y}+\mathbf{d}) - \mathcal{I}(\mathbf{y})\|_\infty}{\|\mathbf{d}\|_\infty} = \left\| \sum_{k=0}^n \frac{d_k}{\|\mathbf{d}\|_\infty}\phi_k\right\|_\infty.\] The absolute condition number maximizes this quantity over all \(\mathbf{d}\). Suppose \(j\) is such that \(\|\phi_j\|_\infty\) is maximal. Letting \(\mathbf{d} = \mathbf{e}_j\), we have the lower bound for \(\kappa(\mathbf{y})\). The upper bound follows from the triangle inequality.Code
The piecewise cubic cardinal function has \(\|\phi\|_\infty \le 1\) which ensures a good condition number for any interpolation with these functions. This is very different than the case for global polynomial interpolation!
Code
The condition number for global polynomial interpolation here is at least 500.
Piecewise linear interpolation
Here, the data points are joined pair-wise by line segments.
Definition: Piecewise linear interpolant
Given nodes \(t_0 < t_1 < \cdots < t_n\), the piecewise linear interpolant \(p(x)\) is given by \[p(x) = y_k + \frac{y_{k+1} - y_k}{t_{k+1}-t_k}(x - t_k) \quad\quad \text{ for } x \in [t_k, t_{k+1}].\]
On each interval \([t_k,t_{k+1}]\), \(p(x)\) is a linear function passing through \((t_k,y_k)\) and \((t_{k+1},y_{k+1})\).
Hat functions
The interpolant is chosen from among linear combinations of a preselected finite set of functions known as hat functions; for \(k = 0, \cdots, n\), \[H_k(x) = \begin{cases} \frac{x - t_{k-1}}{t_k - t_{k-1}} & \text{if } x \in [t_{k-1},t_k], \\ \frac{t_{k+1}-x}{t_{k+1}-t_k} & \text{if } x \in [t_k,t_{k+1}], \\ 0 &\text{otherwise}. \end{cases}\]
The hat functions form a basis for the set of functions that are continuous and piecewise linear relative to the node points \(\mathbf{t}\).
Code
"""
hatfun(t,k)
Create a piecewise linear hat function, where `t` is a vector of n+1
interpolation nodes and `k` is an integer in 0:n giving the index of
the node where the hat function equals one.
"""
function hatfun(t,k)
n = length(t)-1
return function(x)
if k > 0 && t[k] ≤ x ≤ t[k+1]
return (x-t[k])/(t[k+1]-t[k])
elseif k < n && t[k+1] ≤ x ≤ t[k+2]
return (t[k+2] - x)/(t[k+2]-t[k+1])
else
return 0
end
end
endhatfunCardinality conditions
Notice that hat functions are the cardinal functions for piecewise linear interpolation. We can then say \[p(x) = \sum_{k=0}^n y_k H_k(x).\]
Conditioning and convergence
The previous condition number bounds give that the condition number for PL interpolation satisfies \(1 \le \kappa \le n+1\). However, an even stronger results is:
Theorem: Conditioning of PL interpolation
The absolute condition number of piecewise linear interpolation in the infinity norm equals 1. More specifically, if \(\mathcal{I}\) is the piecewise linear interpolation operator, then \[\|\mathcal{I}(\mathbf{y}+\mathbf{z}) - \mathcal{I}(\mathbf{y})\|_\infty = \|\mathbf{z}\|_\infty.\]
This result follows from linearity of \(\mathcal{I}\) and the triangle inequality.
We also ask how close the PL interpolant \(p(x)\) is to function \(f\) that generated \(\mathbf{y}\).
Theorem: Convergence of PL interpolation
Suppose that \(f(x)\) has a continuous second derivative in \([a,b]\) (\(f \in C^2[a,b]\)). Let \(p_n(x)\) be the piecewise linear interpolant of \((t_i,f(t_i))\) for \(i = 0, \cdots, n\), where \(t_i = a + ih\) and \(h = (b-a)/n\). Then \[\|f - p_n\|_\infty = \max_{x\in[a,b]}|f(x) - p_n(x)| \le M h^2,\] where \(M = \|f''\|_\infty\).
Definition: Algebraic convergence
If an approximation has error that is \(O(h^m)\) as \(h \rightarrow 0\) for an integer \(m\) and a discretization size parameter \(h\), then we say the approximation has algebraic convergence. if the error is not also \(O(h^{m+1})\), then \(m\) is the order of accuracy.
PL interpolation is second-order accurate.
Code
┌───────┬────────────┐
│ n │ err │
│ Int64 │ Float64 │
├───────┼────────────┤
│ 10 │ 0.150471 │
│ 100 │ 0.00166421 │
│ 1000 │ 1.66494e-5 │
└───────┴────────────┘Cubic splines
Piecewise cubic interpolation provides a smoother interpolant than PL interpolation, and has continuous derivatives.
Definition: Cubic spline
A cubic spline is a piecewise cubic function that has two continuous derivatives everywhere.
The spline \((x)\) is defined on the interval \([t_{k-1},t_k]\) by cubic polynomial \[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, \quad\quad k=1,\cdots,n,\] where \(a_k, b_k, c_k, d_k\) are values to be determined. Overall there are \(4n\) such undetermined coefficients.
Smoothness conditions
The following constraints guarantee that \(S\) has at least two continuous derivative everywhere.
- \(S_k(t_{k-1}) = y_{k-1}\) and \(S_k(t_k) = y_k\) for every \(k = 1, \cdots, n\).
This condition can be rewritten as \(a_k = y_{k-1}\) and \(a_k + b_k h_k + c_k h_k^2 + d_k h_k^3 = y_k\) for \(k = 1, \cdots, n\). These linear constraints can be written as \[\begin{bmatrix} \mathbf{I} & \mathbf{0} & \mathbf{0} & \mathbf{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},\] and \[\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}\] where \(\mathbf{H}\) is a diagonal matrix with \(h_1, \cdots, h_n\) on the diagonal.
- Continuity of \(S'(x)\) at interior nodes.
This can be written as \(b_k + 2 c_k h_k + 3 d_k h_k^2 = b_{k+1}\) for \(k = 1,, \cdots, n-1\).
This can be rewritten as \[\mathbf{E} \begin{bmatrix} \mathbf{0} & \mathbf{J} & 2\mathbf{H} & 3\mathbf{H}^2 \end{bmatrix} \begin{bmatrix} \mathbf{a} \\ \mathbf{b} \\ \mathbf{c} \\ \mathbf{d} \end{bmatrix} = \mathbf{0},\] where \(\mathbf{J}\) has ones on the diagonal and \(-1\)s on the first diagonal, and zeros elsewhere, and \(\mathbf{E}\) is the \((n-1) \times n\) matrix obtained by deleting the last row from the identity matrix.
- Continuity of \(S''(x)\) at interior nodes.
This can be written as \(2c_k + 6 d_k h_k = 2 c_{k+1}\) for \(k = 1, \cdots, n-1\).
This can be rewritten as \[\mathbf{E} \begin{bmatrix} \mathbf{0} & \mathbf{0} & \mathbf{J} & 3\mathbf{H} \end{bmatrix} \begin{bmatrix} \mathbf{a} \\ \mathbf{b} \\ \mathbf{c} \\ \mathbf{d} \end{bmatrix} = \mathbf{0}.\]
End conditions
So far there are \(4n-2\) conditions on the \(4n\) unknowns. To get a well-defined problem, we add two more constraints. There are two typical approaches:
- Natural spline: \(S_1''(t_0) = S''_n(t_n) = 0\)
- Not-a-knot spline: \(S_1'''(t_1) = S_2'''(t_1)\), \(S_{n-1}'''(t_{n-1}) = S_n'''(t_{n-1})\)
We consider the not-a-knot spline as they give better pointwise accuracy.
Algebraically, these conditions are \(d_1 = d_2\), \(d_{n-1} = d_n\). We append these to the systems we built above.
Note that there is no closed form expression for a cardinal basis for the cubic spline, so the process of computing it requires solving the constructed linear system.
Code
┌───────┬─────────────┐
│ n │ error │
│ Int64 │ Float64 │
├───────┼─────────────┤
│ 8 │ 0.0305634 │
│ 11 │ 0.0207562 │
│ 16 │ 0.00590761 │
│ 23 │ 0.00134587 │
│ 32 │ 0.000367049 │
│ 45 │ 9.17785e-5 │
│ 64 │ 2.15306e-5 │
│ 91 │ 5.04292e-6 │
│ 128 │ 1.24012e-6 │
└───────┴─────────────┘Code
Theorem:
Suppose that \(f(x)\) has four continuous derivatives in \([a,b]\). Let \(S_n(x)\) be the not-a-knot cubic spline interpolant of \((t_i,f(t_i))\) for \(i = 0, \cdots, n\), where \(t_i = a + ih\) and \(h = (b-a)/n\). Then for all sufficiently small \(h\), there is a constant \(C > 0\) such that \[\|f - S_n\|_\infty \le C h^4.\]
