drnxxx

# Compact Finite Differences and Cubic Splines

Robert M. Corless
The Ontario Research Center for Computer Algebra and The School of
Mathematical and Statistical Sciences The University of Western Ontario London
Email<EMAIL_ADDRESS>Ontario  Canada

(15 October 2019)

###### Abstract

In this paper I uncover and explain—using contour integrals and residues—a
connection between cubic splines and a popular compact finite difference
formula. The connection is that on a uniform mesh the simplest Padé scheme for
generating fourth-order accurate compact finite differences gives exactly the
derivatives at the interior nodes needed to guarantee twice-continuous
differentiability for cubic splines. I also introduce an apparently new
spline-like interpolant that I call a compact cubic interpolant; this is
similar to one introduced in 1972 by Swartz and Varga, but has higher order
accuracy at the edges. I argue that for mildly nonuniform meshes the compact
cubic approach offers some potential advantages, and even for uniform meshes
offers a simple way to treat the edge conditions, relieving the user of the
burden of deciding to use one of the three standard options: free (natural),
complete (clamped), or “not-a-knot” conditions. Finally, I establish that the
matrices defining the compact cubic splines (equivalently, the fourth-order
compact finite difference formulæ) are positive definite, and in fact totally
nonnegative, if all mesh widths are the same sign, for instance if the mesh is
real and nodes are increasing only. compact finite differences; cubic splines;
barycentric form; compact cubic splines; contour integral methods; totally
nonnegative matrices

##  1 Introduction

“The most popular choice continues to be a piecewise cubic approximating
function.”

—Carl de Boor (de Boor, 1978, p. 49).

Cubic splines and the similar piecewise interpolant known as pchip are both
widely used in piecewise polynomial interpolation. Cubic splines give a twice
continuously differentiable interpolant through given data, while pchip tries
instead to preserve monotonicity and convexity. Both are useful not just
because they fit the data, but also because their derivatives are generally
also good approximations to the derivative of the underlying function that one
wants to approximate. Loosely speaking, on a uniform mesh of width $h$ the fit
of a cubic spline to a smooth function is accurate to $O(h^{4})$ and the
derivative is accurate to $O(h^{3})$. For more careful error bounds see Swartz
& Varga (1972) and de Boor (1978).

In contrast to both cubic splines and pchip, a true cubic Hermite interpolant
fits not only function values at the nodes but also true derivative values at
the nodes: both cubic splines and pchip use substitutes computed from the
function values. The similarity of names (pchip stands for Piecewise Cubic
Hermite Interpolant) does cause confusion. The relative rarity of the case
when true derivatives $\rho_{k,1}={{f}^{\,\prime}}(\tau_{k})$ are specified
makes this confusion bearable. This paper will concentrate on cubic splines,
and not on pchip or on true cubic Hermite interpolants. The classical
reference for splines is de Boor (1978), but see also chapter 3 of Moler
(2004).

Compact finite differences are an efficient and accurate way of approximating
the derivatives of known data. There are several formulæ in use. See for
instance Lele (1992) or Collatz (1966).

In this paper I detail a connection between cubic splines on a uniform mesh
and a popular compact finite difference formula that the traditional (monomial
basis) method of computing the cubic splines had obscured. Namely, the
simplest Padé scheme for generating fourth-order compact finite differences
gives exactly the derivatives at the interior nodes needed to guarantee twice-
continuous differentiability; that is, a spline. The literature both of
splines and of finite differences is vast. Nonetheless I believe that this
connection is new. In this paper I will give the two derivations, and show
that, in the case of an equally-spaced mesh, equation (11) from the compact
scheme and equation (9) from the cubic spline are identical after
normalization (each could be multiplied by an arbitrary constant without
altering the solution). I then offer an explanation in section 5 of why this
is so (or why it might have been expected to be so).

I also comment on the various choices for treating the degrees of freedom at
the endpoints, and introduce a new choice, which I call a “compact cubic
spline”, namely the use of a modified compact finite difference formula to
give the derivatives at the endpoints. I prove that the matrices involved are
positive definite (indeed totally nonnegative) if the mesh widths have the
same sign111Normally, the mesh will be real with the nodes increasing, so the
mesh widths will all be positive. Sometimes it is convenient to number real
nodes from the rightmost endpoint, so the mesh widths will all be negative. If
the nodes are complex and not in a straight line, then indeed the mesh widths
will not have the same (complex) sign.. This is similar in spirit to the
approach of Swartz & Varga (1972), who use a four-node explicit finite
difference at each end to provide $O(h^{3})$ accurate approximate derivatives
at the edges of a uniform grid. Here, we use a compact four-node formula that
gives $O(h^{4})$ accurate approximate derivatives at the edges, uniform grid
or not, if the function being interpolated has at least five continuous
derivatives.

##  2 Notation

We study two basic problems and their connection in this paper. The first
basic problem under consideration is construction of a piecewise polynomial
interpolant $p(t)$ from function value data $\rho_{i,0}$ on a possibly
nonuniform mesh made up of distinct nodes $\tau_{i}$, numbered $0\leq i\leq
n$. This will lead to $n+1$ by $n+1$ matrices. We will always take $n\geq 4$
because one can fit a single degree $n$ polynomial on $n+1$ points, and if
$n\leq 3$ this can be done with a single cubic.

We postpone discussion of the second basic problem, which is to take
derivatives of data given on a possibly nonuniform mesh, until section 4. For
simplicity in both problems, we restrict to the real case and suppose that
values of a smooth function $f:[a,b]\to\mathbb{R}$ have been sampled on some
partition of a finite interval $[a,b]$, with
$a=\tau_{0}<\tau_{1}<\tau_{2}<\cdots<\tau_{n-1}<\tau_{n}=b$. We write
$\rho_{i,0}=f(\tau_{i})$ for the known function values, for $0\leq i\leq n$.
On each subinterval $[\tau_{k},\tau_{k+1}]$ the piecewise interpolant will be
given by $p(t)=p_{k}(t)$, a polynomial of degree at most $3$ satisfying
$p_{k}(\tau_{k})=\rho_{k,0}$ and $p_{k}(\tau_{k+1})=\rho_{k+1,0}$, and
satisfying some other conditions that we will specify later. We will refer to
this as a piecewise cubic polynomial even if some or all of its degrees are
less than $3$.

As is very well known, both cubic splines and the shape-preserving interpolant
known as pchip construct values $\rho_{i,1}$ at the interior nodes $\tau_{i}$
for $1\leq i\leq n-1$ of the derivative of each $p_{k}(t)$ in order to ensure
continuous differentiability there: we wish
$p_{i-1}^{\prime}(\tau_{i})=p_{i}^{\prime}(\tau_{i})$ for all interior nodes,
$1\leq i\leq n-1$. A spline goes further than pchip and chooses these values
in order to ensure twice continuous differentiability at the interior nodes:
$p_{i-1}^{\prime\prime}(\tau_{i})=p_{i}^{\prime\prime}(\tau_{i})$. By not
imposing this last condition, a pchip can use some degrees of freedom to
preserve monotonicity instead, usually by taking the derivatives $\rho_{i,1}$
to be a harmonic mean of the slopes of the secants on either side (see Moler
(2004)).

##  3 Cubic Splines

What follows in this section is a new derivation, or at least a derivation
likely new to the reader, of the piecewise cubic Hermite spline. This
derivation was apparently first published in Chapter 8 of Corless & Fillion
(2013). The novelty is that it uses the barycentric forms for a cubic Hermite
interpolant on each piece: that is, instead of trying to fit unknown cubic
polynomials in a local monomial basis
$p_{k}(t)=a_{k}+b_{k}(t-\tau_{k})+c_{k}(t-\tau_{k})^{2}+d_{k}(t-\tau_{k})^{3}$
to the data and trying to find reasonable ways to determine the $4n$ unknowns
(there are $n$ subintervals), one instead works directly with

$\displaystyle
p_{k}(t)=\frac{\displaystyle\sum_{i=k}^{k+1}\sum_{j=0}^{1}\sum_{\ell=0}^{j}\beta_{i,j}\rho_{i\ell}(t-\tau_{i})^{\ell-j-1}}{\displaystyle\sum_{i=k}^{k+1}\sum_{j=0}^{1}\beta_{i,j}(t-\tau_{i})^{-j-1}}\>,$
(1)

which is the second barycentric form of the cubic Hermite interpolant (note
that only two nodes are used in this form and thus this determines only one
piece of the interpolant). There are four unknowns in this formula:
$\rho_{k,0}$, $\rho_{k,1}$, $\rho_{k+1,0}$, and $\rho_{k+1,1}$. The
$\beta_{i,j}$ can be found from the partial fraction decomposition of the
reciprocal of the node polynomial, and thus can be regarded as known once the
nodes $\tau_{k}$ are specified:

$\frac{1}{(t-\tau_{k+1})^{2}(t-\tau_{k})^{2}}=\sum_{i=k}^{k+1}\sum_{j=0}^{1}\beta_{i,j}(t-\tau_{i})^{-j-1}\\\
=\frac{\displaystyle\frac{-2}{(\tau_{k+1}-\tau_{k})^{3}}}{t-\tau_{k}}+\frac{\displaystyle\frac{1}{(\tau_{k+1}-\tau_{k})^{2}}}{(t-\tau_{k})^{2}}+\frac{\displaystyle\frac{2}{(\tau_{k+1}-\tau_{k})^{3}}}{t-\tau_{k+1}}+\frac{\displaystyle\frac{1}{(\tau_{k+1}-\tau_{k})^{2}}}{(t-\tau_{k+1})^{2}}\>.$

There are only four $\beta_{i,j}$ for each interval, and we see them written
above explicitly in terms of the given nodes $\tau_{k}$. For convenience, one
can simplify the second barycentric form to the usual cubic Hermite polynomial
basis (here $h_{k+1}=\tau_{k+1}-\tau_{k}$ is the width of the interval):

$p(t)={\frac{\left(t-\tau_{{k+1}}\right)^{2}\left(2\,t-2\,\tau_{{k}}+h_{{k+1}}\right)\rho_{{k,0}}}{{h_{{k+1}}}^{3}}}+{\frac{\left(t-\tau_{{k+1}}\right)^{2}\left(t-\tau_{{k}}\right)\rho_{{k,1}}}{{h_{{k+1}}}^{2}}}\\\
+{\frac{\left(t-\tau_{{k}}\right)^{2}\left(2\,\tau_{{k+1}}+h_{{k+1}}-2\,t\right)\rho_{{k+1,0}}}{{h_{{k+1}}}^{3}}}+{\frac{\left(t-\tau_{{k}}\right)^{2}\left(t-\tau_{{k+1}}\right)\rho_{{k+1,1}}}{{h_{{k+1}}}^{2}}}\>.$
(2)

######  Remark 3.1.

Evaluation of these formulæ and their derivatives does not entail significant
rounding error. A standard backward error analysis (not given here) shows that
if IEEE standard floating-point arithmetic is used, then (using the notation
of Higham (2002) ) $\mathrm{fl}\left(p(x)\right)$ is the exact value of a
polynomial going through the data
$[\rho_{k,0}(1+\theta_{6}),\rho_{k,1}(1+\theta_{6}),\rho_{k+1,0}(1+\theta_{5}),\rho_{k+1,1}(1+\theta_{5})]$;
that is the floating-point evaluation of $p(x)$ is the exact value of a
polynomial going through data that is at most six rounding errors different to
the original data.

Notice that the $\rho_{i,0}$ (not the $\rho_{i,1}$, which represent derivative
values) are the known data values. We want to choose the $n+1$ slopes
$\rho_{i,1}$ to make the resulting interpolant as smooth as possible. We will
see that we can make it $\mathscr{C}^{2}[\tau_{0},\tau_{n}]$; that is, the
second derivative will be continuous at each interior node. This makes the
piecewise interpolant a cubic spline, albeit represented in a different basis.
Notice also that we may choose the $\rho_{i,1}$ in such a way that we
automatically have $p(t)\in\mathscr{C}^{1}[\tau_{0},\tau_{n}]$: just take the
slope at the right end of one interval to be the same slope at the left end of
the next. This is very natural because $\rho_{i,1}$ is then interpreted as
‘the’ slope at the node $\tau_{i}$ (indeed it would be somewhat unnatural to
have different slopes on the left and right, though we could do that if we
wanted). Having made our interpolant continuously differentiable by this
device, then
$p_{k-1}^{\prime}(\tau_{k}^{-})=\rho_{k,1}=p^{\prime}_{k}(\tau_{k}^{+})$. To
further ensure $p(t)\in\mathscr{C}^{2}[\tau_{0},\tau_{n}]$, we want to make
the second derivatives equal, i.e.,

$\displaystyle
p^{\prime\prime}_{k-1}(\tau^{-}_{k})=p^{\prime\prime}_{k}(\tau_{k}^{+})\qquad
k=1,2,\ldots,n-1\>.$ (3)

I personally found this easier than the necessary algebra in the local
monomial case. We would have

$\displaystyle
p_{k}(t)=\rho_{k,0}+\rho_{k,1}(t-\tau_{k})+c_{k}(t-\tau_{k})^{2}+d_{k}(t-\tau_{k})^{3}$
(4)

and even to make the function just $\mathscr{C}^{1}$ we would have to impose
the condition

$\displaystyle
p_{k}^{\prime}(\tau^{-}_{k+1})=p^{\prime}_{k+1}(\tau_{k+1}^{+})\>,$

which isn’t automatic; we would have to enforce

$\displaystyle\rho_{k,1}+2c_{k}(\tau_{k+1}-\tau_{k})+3d_{k}(\tau_{k+1}-\tau_{k})^{2}=\rho_{k+1,1}\>.$

For $\mathscr{C}^{2}$, we would have to enforce yet another condition, namely

$\displaystyle 2c_{k}+6d_{k}(\tau_{k+1}-\tau_{k})=2c_{k+1}\>.$

Of course, this can be done, and the solution is even elegant. These equations
can be reduced to a tridiagonal system of equations for the slopes
$\rho_{k,1}$, and explicit formulæ for the $c_{k}$ and $d_{k}$ are known once
the slopes are known. The solution is shown in splinetx.m in Moler (2004).

But here, because we start with the Hermite interpolational basis, we have a
simpler (but equivalent) task: just enforce the second derivative conditions.
To do this, we need a formula for the second derivative of the cubic Hermite
interpolant at the nodes. A short computation by hand or with a computer
algebra language such as Maple shows that

$\displaystyle
p_{k-1}^{\prime\prime}(\tau_{k}^{-})=\frac{2}{\tau_{k}-\tau_{k-1}}(2\rho_{k,1}+\rho_{k-1,1})-\frac{6}{(\tau_{k}-\tau_{k-1})^{2}}(\rho_{k,0}-\rho_{k-1,0})$
(5)

and

$\displaystyle
p_{k}^{\prime\prime}(\tau_{k}^{+})=\frac{-2}{\tau_{k+1}-\tau_{k}}(2\rho_{k,1}+\rho_{k+1,1})+\frac{6}{(\tau_{k+1}-\tau_{k})^{2}}(\rho_{k+1,0}-\rho_{k,0})\>.$
(6)

Equating these at interior nodes $1\leq k\leq n-1$ gives $n-1$ equations
constraining the slopes. Explicitly, these equations are

$\frac{2}{\tau_{k}-\tau_{k-1}}\rho_{k-1,1}+4\left(\frac{1}{\tau_{k}-\tau_{k-1}}+\frac{1}{\tau_{k+1}-\tau_{k}}\right)\rho_{k,1}+\frac{2}{\tau_{k+1}-\tau_{k}}\rho_{k+1,1}\\\
=\frac{6}{(\tau_{k+1}-\tau_{k})^{2}}(\rho_{k+1,0}-\rho_{k,0})+\frac{6}{(\tau_{k}-\tau_{k-1})^{2}}(\rho_{k,0}-\rho_{k-1,0})$
(7)

for $k=1,2,3,\ldots,n-1$. There are $n-1$ equations. The structure of the
resulting matrix is tridiagonal. These are equivalent to those of equation (9)
of (de Boor, 1978, p. 54). It is convenient to introduce
$h_{k}=\tau_{k}-\tau_{k-1}$ and $h_{k+1}=\tau_{k+1}-\tau_{k}$ and multiply
equation (7) by $h_{k}h_{k+1}/(h_{k}+h_{k+1})$ to get

$\frac{2h_{k+1}}{h_{k}+h_{k+1}}\rho_{k-1,1}+4\rho_{k,1}+\frac{2h_{k}}{h_{k}+h_{k+1}}\rho_{k+1,1}\\\
=\frac{6h_{k}}{h_{k+1}(h_{k}+h_{k+1})}(\rho_{k+1,0}-\rho_{k,0})+\frac{6h_{k+1}}{h_{k}(h_{k}+h_{k+1})}(\rho_{k,0}-\rho_{k-1,0})\>.$
(8)

If the mesh is equally-spaced, that is $\tau_{k}=a+(b-a)k/n=a+kh$ for $0\leq
k\leq n$, with $h=(b-a)/n$, then these equations become for $1\leq i\leq n-1$

$\displaystyle\rho_{i-1,1}+4\rho_{i,1}+\rho_{i+1,1}$
$\displaystyle=\frac{3}{h}\left(\rho_{i+1,0}-\rho_{i-1,0}\right)\>.$ (9)

If we specify two edge conditions and then solve this tridiagonal system of
$n-1$ equations, we will get derivative values at the interior nodes that will
ensure a twice continuously differentiable piecewise interpolant, where each
piece $p_{k}(t)$ is given by equation (2).

This leaves two degrees of freedom, which (counterintuitively) may be
unwanted: the user will have to choose what to do with these two degrees of
freedom. Some common choices are, first, to ask for three times continuous
differentiability at $\tau_{1}$ and at $\tau_{n-1}$: this is called the “not-
a-knot” condition. Second, one can specify arbitrary values of the derivatives
at the ends, which is sometimes called a “clamped” spline; if the derivative
values are correct for the function being interpolated, it is called a
complete spline in de Boor (1978). Finally, one can ask for a second
derivative of zero at the ends, called the “natural” spline.

We will see a new choice, actually a new set of choices, similar to a fifth
choice due to Swartz & Varga (1972) that is mentioned in de Boor (1978), after
the next section.

##  4 Compact finite differences

This section is adapted from chapter 11 of Corless & Fillion (2013). The
problem being addressed is the problem of computing derivatives of a function
known only at discrete values; we are not necessarily interested in computing
an interpolant. One could use, for instance, simple finite differences and
approximate $f^{\prime}(\tau_{k})$ by
$(f(\tau_{k+1})-f(\tau_{k}))/(\tau_{k+1}-\tau_{k})$. The main idea of a
compact finite difference is that instead of using a single explicit finite
difference formula to evaluate a derivative at a point, we have a whole mesh
of function values and we wish to compute the derivatives at all the nodes.
This is quite like the case of a global interpolating polynomial where one
constructs a differentiation matrix (see Weideman & Reddy (2000)), and indeed
we will have the equivalent of a differentiation matrix here; but it will not
be explicitly formed. Instead we will solve a banded linear system.

The purpose of this paper is to show that the canonical example of a compact
finite difference happens to have a significant relationship with cubic
splines. This fact does not seem to have been noticed before, although it is
difficult to be sure, given the variation in nomenclature and the large number
of works on piecewise polynomials of one kind or another. Specifically, the
relationship is as follows.

Let $\Delta$ be the difference operator $\Delta(f)=f(x+h)-f(x)$. Then the so-
called operator approach to finite differences (which appears for instance in
Boole (1960)) gives a relationship between the differentiation operator $D$
and $\Delta$, namely $D=\ln(1+\Delta)/h$. Then the $(2,2)$ Padé approximant
for $\ln(1+\Delta)$, namely

$D\doteq{\frac{\Delta+{\Delta}^{2}/2}{1+\Delta+{\Delta}^{2}/6}}\>,$ (10)

gives rise to a compact finite difference scheme (also called a Padé scheme)
that happens to be fourth-order accurate: apply the denominator
$1+\Delta+{\Delta}^{2}/6$ to ${{f}^{\,\prime}}$ and the numerator
$\Delta+\Delta^{2}/2$ to $f$ and we find after shifting to center at $t$ that,
up to terms of $O(h^{4})$,

$\displaystyle\frac{1}{6}{{f}^{\,\prime}}(t-h)+\frac{2}{3}{{f}^{\,\prime}}(t)+\frac{1}{6}{{f}^{\,\prime}}(t+h)=\frac{1}{2h}\left(f(t+h)-f(t-h)\right)\>.$

We will discuss further the operator approach briefly in section 6.2. This
formula gives us a tridiagonal (whence “compact”) system of equations for the
unknown derivatives, at all interior nodes; again something special has to be
done at the edges. This tridiagonal system of equations turns out to be
identical up to scaling—when the nodes are equally-spaced—to equations (7).

The idea of solving a system of equations to find our finite difference
approximation to the derivatives may be unfamiliar. The simplest finite
difference formulæ make straightforward linear combinations such as
$(f(t+h)-f(t-h))/(2h)$ to approximate $f^{\prime}(t)$. But compact finite
difference formulæ are different. That is, instead of simply applying a
formula to a vector of function values to get a vector of derivative values,
we instead have to set up and solve a linear system of equations for the
unknown derivatives. Having to solve equations is more complicated than just
using a formula, but it has several advantages.

To better understand where the system of equations for this formula comes
from, we give another derivation. Make the following simplifying assumptions,
for the moment. Suppose ${{f}^{\,\prime}}(\tau_{0})$ and
${{f}^{\,\prime}}(\tau_{n})$ are known (just to make it simple) and that
$\tau_{k+1}-\tau_{k}=h$ is constant. Then fix attention on one particular
node, say $\tau_{k}$. The formula above becomes, when multiplied by $6$ and
putting $t-h=\tau_{k-1}$, $t=\tau_{k}$, and $t+h=\tau_{k+1}$,

$\displaystyle{{f}^{\,\prime}}(\tau_{k-1})+4{{f}^{\,\prime}}(\tau_{k})+{{f}^{\,\prime}}(\tau_{k+1})=\frac{3}{h}(f(\tau_{k+1})-f(\tau_{k-1}))\>.$
(11)

This equation is exactly equation (9). Recall that $\rho_{k,1}$ represents
derivatives $f^{\prime}(\tau_{k})$ and $f(\tau_{k})=\rho_{k,0}$. Now we let
$k$ vary over all the indices of the interior nodes, $1\leq k\leq n-1$. Each
interior node gives us one equation. Each equation only contains at most three
of the unknown derivatives (and the equation for $k=1$ touches the known
derivative ${{f}^{\,\prime}}(\tau_{0})$, while the equation for $k=n-1$
touches the known derivative ${{f}^{\,\prime}}(\tau_{n})$). This gives us a
tridiagonal linear system of equations to solve for the unknown derivatives
${{f}^{\,\prime}}(\tau_{k})$, $1\leq k\leq n-1$. Call the tridiagonal matrix
$\bm{\mathbf{A}}$. Notice also that the right-hand side of the system involves
linear combinations of the values of $f(\tau_{k})$ at different nodes—these
are supposed to be known. Call the (also tridiagonal) matrix that forms that
combination, $\bm{\mathbf{B}}$. Note that $\bm{\mathbf{B}}$ has a zero
diagonal. The system $\bm{\mathbf{A}}\bm{\mathbf{v}}=\bm{\mathbf{B}}\bm{\rho}$
needs to be solved computationally to get the vector $\bm{\mathbf{v}}$ of
desired derivatives.

In effect this computes the differentiation matrix $\bm{\mathbf{D}}$ as
$\bm{\mathbf{A}}^{-1}\bm{\mathbf{B}}$, but in practice one never explicitly
computes $\bm{\mathbf{A}}^{-1}$ because it is a full matrix. Instead, of
course, to compute $\bm{\mathbf{D}}\bm{\mathbf{y}}$, one solves
$\bm{\mathbf{A}}\bm{\mathbf{v}}=\bm{\mathbf{B}}\bm{\mathbf{y}}$ for
$\bm{\mathbf{v}}$, which is formally
$\bm{\mathbf{A}}^{-1}\bm{\mathbf{B}}\bm{\mathbf{y}}$.

Figure 1:  Fourth-order uniform mesh compact finite difference derivative
graphed with the exact derivative (dashed line). Even with just five
subintervals, visual accuracy is achieved in this example.

I now give a derivation of a fourth-order compact finite difference method for
nonuniform meshes, using contour integrals because this also works for the
barycentric form. This provides a point of theoretical continuity with section
3, and will be useful for the explanation of the connection to cubic splines
in section 5. I will also give a similar derivation for a compact formula
valid at the endpoints, which gives us a new way to deal with the extra
degrees of freedom in a cubic spline.

Consider the following partial fraction expansion: here we put $h_{k-1}=rh$
and $h_{k}=sh$ and we assume both $r$ and $s$ are different from zero (we
could without loss of generality assume that one of $r$ or $s$ was $1$ but
this doesn’t help much, and I find the symmetry in the formulæ below useful
for understanding).

$\frac{1}{(z+rh)^{2}z^{2}(z-sh)^{2}}=\frac{1}{{r}^{2}{h}^{4}\left(s+r\right)^{2}\left(z+rh\right)^{2}}+\frac{4\,r+2\,s}{\left(z+rh\right){r}^{3}{h}^{5}\left(s+r\right)^{3}}\\\
+\frac{1}{{r}^{2}{h}^{4}{s}^{2}{z}^{2}}+\frac{2\,r-2\,s}{z{r}^{3}{h}^{5}{s}^{3}}+\frac{1}{{h}^{4}\left(s+r\right)^{2}{s}^{2}\left(z-sh\right)^{2}}+\frac{-4\,s-2\,r}{\left(z-sh\right){h}^{5}\left(s+r\right)^{3}{s}^{3}}\>,$

By a standard argument in complex variables, the contour integral of

$\frac{1}{2\pi i}\oint_{C}\frac{f(z)}{(z+rh)^{2}z^{2}(z-sh)^{2}}\,dz$ (12)

over the contour that encloses all the zeros of the denominator (the numerator
is a polynomial and has no poles) is zero for polynomials $f(x)$ of degree
$6-2=4$ or less. Once we expand this using the partial fraction above and use
the Cauchy Integral formula

$\frac{f^{(j)}(p)}{j!}=\frac{1}{2\pi i}\oint_{C}\frac{f(z)}{(z-p)^{j+1}}\,dz$
(13)

and multiply by $(s+r)^{2}$ for convenience, this gives us the following
formula, which is exact for polynomials of degree at most four222The fact that
this is valid for polynomials of degree four is an important point in
understanding why the equations are the same as for cubic splines, as we will
see in section 5..

${\frac{1}{{r}^{2}}}{{f}^{\,\prime}}(-rh)+{\frac{\left(s+r\right)^{2}}{{r}^{2}{s}^{2}}}{{f}^{\,\prime}}(0)+{\frac{1}{{s}^{2}}}{{f}^{\,\prime}}(sh)\\\
=-{\frac{4\,r+2\,s}{{r}^{3}{h}(r+s)}}f(-rh)-{\frac{2\left(r-s\right)\left(s+r\right)^{2}}{{r}^{3}{h}{s}^{3}}}f(0)+{\frac{4\,s+2\,r}{{h}(r+s){s}^{3}}}f(sh)\>.$
(14)

These equations, if $r\neq s$, are not the same as equations (8), whatever
constant we use to normalize them. Setting $r=s=1$, however, gives equation
(11). If we now have a mesh of distinct points
$\tau_{0}<\tau_{1}<\cdots<\tau_{n-1}<\tau_{n}$, we can lay this compact
formula first on $\tau_{0}$, $\tau_{1}$, and $\tau_{2}$, with
$rh=\tau_{1}-\tau_{0}$ and $sh=\tau_{2}-\tau_{1}$ which, if we have a
reference step width $h$, say $h=(\sum(\tau_{k+1}-\tau_{k}))/n$, gives us an
equation relating the derivative values on these three mesh points to the
function values on the mesh points. We then lay the formula over $\tau_{1}$,
$\tau_{2}$, and $\tau_{3}$, giving us another equation—one for each interior
point, $\tau_{1}$, $\tau_{2}$, $\ldots$, $\tau_{n-1}$. The linear system for
the unknown derivative values is, like the spline equations, tridiagonal; but
we only have $n-1$ equations and $n+1$ unknowns ${{f}^{\,\prime}}(\tau_{0})$,
${{f}^{\,\prime}}(\tau_{1})$, $\ldots$, ${{f}^{\,\prime}}(\tau_{n})$. To make
a square system, we need two more equations.

###  4.1 Truncation error in the compact formula

Once the formula has been found, one can do a more standard Taylor series
analysis on it, for instance in Maple. Taking the Taylor series in $h$ of the
difference $L-R$ where $L$ is the left-hand side of equation (14) and $R$ is
the right-hand side, the error term turns out to be

$L-R={\frac{\left(s+r\right)^{2}}{120}}{h}^{4}f^{(5)}(0)-{\frac{\left(r-s\right)\left(s+r\right)^{2}}{360}}{h}^{5}f^{(6)}(0)+O(h^{6})\>,$
(15)

which is $O(h^{4})$ regardless of the local mesh ratios $r$ and $s$. If
$r=s=1$ then this becomes equation (11) and the leading error term becomes
$h^{4}/30$ times the fifth derivative of $f$ evaluated at zero.

###  4.2 Truncation error in the spline formula, considered as a compact
finite difference

We can do the same kind of error analysis as we have just done for the compact
formula, but for equations (8). Put $\rho_{k,1}=f^{\prime}(\tau_{k})$ and
$\rho_{k,0}=f(\tau_{k})$ in equations (8). Take $h_{k-1}=rh$ and $h_{k}=sh$.
For a sufficiently differentiable function $f$ we may take the Taylor series
in $h$ of the difference $L-R$ between the left- and right-hand sides. Then
the result is

$L-R=-\frac{rs(r-s)}{12}f^{(4)}(0)h^{3}+\frac{rs(r^{2}-rs+s^{2})}{30}f^{(5)}(0)h^{4}+O(h^{5})\>.$
(16)

This is the truncation error that occurs if the spline continuity equations
(7) are used as a compact finite difference formula. We see that if the mesh
is not locally approximately uniform, that is $r\neq s+O(h)$, then this error
is $O(h^{3})$, not $O(h^{4})$. We may compare this with equation (15); we see
that if $r=s$ they are identical, as we will see that they must be.

###  4.3 What to do at the edges

We now return to computing fourth-order accurate compact finite differences.
We can use the same contour integral method to look for fourth-order formulae
at either end, giving equations involving $\tau_{0}$ and its nearest mesh
neighbours, and $\tau_{n}$ and its nearest neighbours. We will want the
formulæ likewise exact for polynomials of degree four or less. Since the
matrix is so far tridiagonal, we try to keep it that way and we thus look for
relations of the form

$a_{0}{{f}^{\,\prime}}(\tau_{0})+b_{0}{{f}^{\,\prime}}(\tau_{1})=c_{0}f(\tau_{0})+c_{1}f(\tau_{1})+c_{2}f(\tau_{2})+c_{3}f(\tau_{3})\>.$
(17)

The reader may verify that fourth order accuracy is not possible in general
without the function value at $\tau_{3}$. This still qualifies as ‘compact’
though because we use only two extra mesh points at the left end, and
similarly only two extra on the right, and these appear in the right-hand side
and do not change the tridiagonality of the matrix. We will need $n\geq 4$,
giving five by five matrices at the smallest. This ansatz suggests looking at
the partial fraction decomposition of

${\frac{1}{\left(z-\tau_{{0}}\right)^{2}\left(z-\tau_{{1}}\right)^{2}\left(z-\tau_{{2}}\right)\left(z-\tau_{{3}}\right)}}\>,$
(18)

from which we straightforwardly find (of course by using a computer algebra
system) that (with $h_{1}=\tau_{1}-\tau_{0}$, $h_{2}=\tau_{2}-\tau_{1}$, etc
and normalizing so $b_{0}=1$)

$\displaystyle
a_{{0}}={\frac{h_{2}\left(h_{2}+h_{3}\right)}{\left(h_{1}+h_{2}\right)\left(h_{1}+h_{2}+h_{3}\right)}}\quad\mathrm{and}\quad
b_{{0}}=1\>,$

and that the $c_{k}$s are

$\displaystyle c_{0}$
$\displaystyle=-{\frac{h_{2}\left(h_{2}+h_{3}\right)\left(4\,{h_{1}}^{2}+6\,h_{1}h_{2}+3\,h_{1}h_{3}+2\,{h_{2}}^{2}+2\,h_{{1}}h_{3}\right)}{h_{1}\left(h_{1}+h_{2}\right)^{2}\left(h_{1}+h_{2}+h_{3}\right)^{2}}}$
$\displaystyle c_{1}$
$\displaystyle={\frac{2\,h_{2}(h_{2}-h_{1})+h_{3}(2h_{2}-h_{1})}{h_{1}h_{2}\left(h_{2}+h_{3}\right)}}$
$\displaystyle c_{2}$
$\displaystyle={\frac{{h_{1}}^{2}\left(h_{2}+h_{3}\right)}{h_{2}\left(h_{1}+h_{2}\right)^{2}h_{3}}}$
$\displaystyle c_{3}$
$\displaystyle=-{\frac{{h_{1}}^{2}h_{2}}{h_{3}\left(h_{2}+h_{3}\right)\left(h_{1}+h_{2}+h_{3}\right)^{2}}}\>.$

It turns out that the residual error in (17) is, as desired, $O(h^{4})$. In
detail, if $\tau_{k+1}=\tau_{k}+r_{k+1}h$, for $k=0$, $1$, and $2$, and
$h_{1}=\tau_{1}-\tau_{0}=h$, then the residual error is seen by a Taylor
series computation to be

$L-R=\frac{(r_{1}+r_{2})r_{1}}{120}h^{4}f^{(5)}(0)+O(h^{5})\>.$ (19)

In the uniform mesh case this reduces to $h^{4}/60$ times the fifth
derivative. A similar formula holds for the other end (indeed, simply reverse
the labels, $\tau_{k}\leftrightarrow\tau_{n-k}$). This gives us closure in our
search for a compact, variable mesh fourth-order finite difference formula. We
are left with a tridiagonal matrix $\bm{\mathbf{A}}$ with entries depending on
the $h_{k}=\tau_{k}-\tau_{k-1}$.

Apart from the first and last rows, this matrix is diagonally dominant because
$4-4(h_{k}^{2}+h_{k-1}^{2})/(h_{k}+h_{k-1})^{2}=8h_{k}h_{k-1}/(h_{k}+h_{k-1})^{2}$.
We will see that it is indeed positive definite if the mesh widths have the
same sign. For meshes with widely varying widths, however, the matrix can be
ill-conditioned. See figure 2.

By direct computation with Maple for small examples, we see that the
determinant is a rational polynomial in all the $h_{k}$ with positive
coefficients in the numerator and a squared denominator, and indeed is
positive definite. Monitoring the sign of the determinant of $\bm{\mathbf{A}}$
as computed numerically by Matlab’s det function, we find that for $n\geq 4$
the determinant is always reported as being positive. In fact, more is true.

######  Theorem 4.1.

If all mesh widths $h_{k}$ have the same sign, then for $n\geq 4$ the $n+1$ by
$n+1$ matrix $\bm{\mathbf{A}}$ of the compact finite difference formula with
these fourth order formulæ at the edges is totally nonnegative: that is, the
determinants of all minors are nonnegative. More, the determinants of all
principal minors are strictly positive, so the matrix is positive definite.

######  Proof 4.2.

Without loss of generality we will assume all $h_{k}>0$. By a theorem of
Gantmacher and Krein (see (Hogben, 2006, Chapter 29)), an irreducible
tridiagonal matrix is totally nonnegative if and only if its entries are
nonnegative and its leading principal minors are nonnegative. Here all entries
are nonnegative. To prove the theorem we must then verify only that all
leading principal minors are positive. We do so by induction, for all but the
last, and then handle that specially.

For simplicity, write
$\bm{\mathbf{A}}=\bm{\mathbf{Z}}\bm{\mathbf{T}}^{n+1,n+1}$ where
$\bm{\mathbf{Z}}=\mathrm{diag}(1,4/(h_{1}+h_{2})^{2},4/(h_{2}+h_{3})^{2},\ldots,4/(h_{n-1}+h_{n})^{2},1)$
to remove the common denominator of each of the interior rows. Multiplying by
positive diagonal matrices leaves the properties of total positivity
(nonnegativity) undisturbed. This leaves each interior row as
$h_{k}^{2},(h_{k}^{2}+h_{k-1})^{2},h_{k-1}^{2}$, for easier manipulation.
Explicitly writing the six-by-six case, $\bm{\mathbf{T}}^{6,6}=$

$\left[\begin{array}[]{cccccc}{\frac{h_{{2}}\left(h_{{3}}+h_{{2}}\right)}{\left(h_{{2}}+h_{{1}}\right)\left(h_{{3}}+h_{{2}}+h_{{1}}\right)}}&1&0&0&0&0\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr{h_{{2}}}^{2}&\left(h_{{2}}+h_{{1}}\right)^{2}&{h_{{1}}}^{2}&0&0&0\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr
0&{h_{{3}}}^{2}&\left(h_{{3}}+h_{{2}}\right)^{2}&{h_{{2}}}^{2}&0&0\\\ \vskip
6.0pt plus 2.0pt minus 2.0pt\cr
0&0&{h_{{4}}}^{2}&\left(h_{{3}}+h_{{4}}\right)^{2}&{h_{{3}}}^{2}&0\\\ \vskip
6.0pt plus 2.0pt minus 2.0pt\cr
0&0&0&{h_{{5}}}^{2}&\left(h_{{4}}+h_{{5}}\right)^{2}&{h_{{4}}}^{2}\\\ \vskip
6.0pt plus 2.0pt minus 2.0pt\cr
0&0&0&0&1&{\frac{h_{{4}}\left(h_{{3}}+h_{{4}}\right)}{\left(h_{{4}}+h_{{5}}\right)\left(h_{{3}}+h_{{4}}+h_{{5}}\right)}}\end{array}\right]\>.$

Let $\bm{\mathbf{T}}^{k,n+1}=\bm{\mathbf{T}}[1{:}k,1{:}k]$ denote the leading
principal submatrices. Let $D^{k,n+1}=\det(\bm{\mathbf{T}}^{k,n+1})$ for
$1\leq k\leq n+1$. By direct computation,

$\displaystyle D^{1,n+1}$
$\displaystyle={\frac{h_{{2}}\left(h_{{3}}+h_{{2}}\right)}{\left(h_{{2}}+h_{{1}}\right)\left(h_{{3}}+h_{{2}}+h_{{1}}\right)}}$
$\displaystyle D^{2,n+1}$
$\displaystyle={\frac{h_{{1}}h_{{2}}h_{{3}}\left(h_{{2}}+h_{{1}}\right)}{\left(h_{{2}}+h_{{1}}\right)\left(h_{{3}}+h_{{2}}+h_{{1}}\right)}}$
$\displaystyle D^{3,n+1}$
$\displaystyle={\frac{h_{{1}}h_{{3}}{h_{{2}}}^{2}\left(h_{{3}}+h_{{2}}\right)\left(h_{{3}}+h_{{2}}+h_{{1}}\right)}{\left(h_{{2}}+h_{{1}}\right)\left(h_{{3}}+h_{{2}}+h_{{1}}\right)}}$
(20)

######  Remark 4.3.

If $n=3$, then $D^{4,4}=0$. This is excluded from the theorem; as discussed
previously, in that case, just fit a polynomial to all the data.

Returning to the proof, we see that $D^{1,n+1}>0$, $D^{2,n+1}>0$, and
$D^{3,n+1}>0$, for any choice of positive $h_{k}$.

Laplace expansion about the last row of $\bm{\mathbf{T}}^{k,n+1}$ for $3\leq
k\leq n$ gives a recursive formula for the determinant:

$D^{k,n+1}=(h_{k-1}+h_{k})^{2}D^{k-1,n+1}-h_{k-2}^{2}h_{k}^{2}D^{k-2,n+1}\>.$
(21)

For the last minor, namely the determinant of the full matrix, we have

$D^{n+1,n+1}=\frac{h_{n-1}(h_{n-1}+h_{n-2})}{(h_{n-1}+h_{n})(h_{n}+h_{n-1}+h_{n-2})}D^{n,n+1}-h_{n-1}^{2}D^{n-1,n+1}\>.$
(22)

By direct computation,

$\displaystyle D^{3,n+1}$ $\displaystyle=h_{2}^{2}D^{2,n+1}+q_{2}$ (23)

where $q_{2}$ is the positive function

$q_{2}=\frac{h_{1}h_{2}^{2}h_{3}^{2}(h_{1}+2h_{2}+h_{3})}{(h_{1}+h_{2})(h_{1}+h_{2}+h_{3})}\>.$
(24)

Now assume inductively that $D^{k-1,n+1}=h_{k-2}^{2}D^{k-2,n+1}+q_{k-2}$ for a
positive function $q_{k-2}$. We see above that this is true for $k-1=3$, or
$k=4$. We will simultaneously establish by induction that $D^{k,n+1}>0$ and
that $D^{k,n+1}=h_{k-1}^{2}D^{k-1,n+1}+q_{k-1}$ for another positive function
$q_{k-1}$, for $k\leq n$.

######  Remark 4.4.

The assumption $n\geq 4$ of the theorem is used here: we cannot start the
induction otherwise. If $n=3$ then $D^{4,4}$ is not given by the interior
formula, but by the end formula.

Returning to the proof, by the interior recurrence relation above,

$\displaystyle D^{k,n+1}$
$\displaystyle=(h_{k}+h_{k-1})^{2}D^{k-1,n+1}-h_{k}^{2}h_{k-2}^{2}D^{k-2,n+1}$
$\displaystyle=h_{k-1}^{2}D^{k-1,n+1}+q_{k-1}$ (25)

where, by the inductive assumption and expanding $(h_{k}+h_{k-1})^{2}$,

$q_{k-1}=2h_{k}h_{k-1}D^{k-1,n+1}+h_{k}^{2}q_{k-2}$ (26)

is clearly positive, being the sum of positive terms.

This establishes by induction that all leading minors $D^{k,n+1}$ up to, but
not including, the final determinant $D^{n+1,n+1}$, are positive. We verify
that the final determinant is positive by a different method.

If the $(n+1,n+1)$ entry in $\bm{\mathbf{T}}^{n+1,n+1}$ were $\alpha$, then
the determinant $D^{n+1,n+1}$ would be positive if

$\alpha D^{n,n+1}-h_{n-1}^{2}D^{n-1,n+1}>0\>;$ (27)

or, on rearranging,

$\displaystyle\alpha$ $\displaystyle>\frac{h_{n-1}^{2}D^{n-1,n+1}}{D^{n,n+1}}$
$\displaystyle>\frac{h_{n-1}^{2}}{(h_{n}+h_{n-1})^{2}-h_{n}^{2}h_{n-2}^{2}D^{n-2,n+1}/D^{n-1,n+1}}$
(28)

by equation (21). Examining the ratio $D^{n-2,n+1}/D^{n-1,n+1}$ will lead to
the sufficient condition

$\alpha>\frac{h_{n-1}^{2}}{(h_{n-1}+h_{n})^{2}-h_{n}^{2}h_{n-2}/(h_{n-2}+2h_{n-1})}$
(29)

for $D^{n+1,n+1}$ to be positive. In detail, because $D^{n-2,n+1}$ and
$D^{n-1,n+1}$ are interior determinants, they satisfy

$\frac{D^{n-1,n+1}}{D^{n-2,n+1}}=h_{n-2}^{2}+2h_{n-2}h_{n-1}+\frac{h_{n}^{2}q_{n-2}}{D^{n-2,n+1}}$
(30)

where that last term is positive. Dropping the last term and replacing $=$ by
$>$ then following the manipulations (inversion, negation, inversion) through
the algebra of inequalities gives the upper bound

$\frac{h_{n-1}^{2}}{(h_{n}+h_{n-1})^{2}-h_{n}^{2}h_{n-2}^{2}D^{n-2,n+1}/D^{n-1,n+1}}\\\
<\frac{h_{n-1}^{2}}{(h_{n}+h_{n-1})^{2}-h_{n}^{2}h_{n-2}^{2}/(h_{n-2}^{2}+2h_{n-2}h_{n-1})}.$
(31)

Verifying that
$\alpha=\frac{h_{n-1}(h_{n-1}+h_{n-2})}{(h_{n-1}+h_{n})(h_{n}+h_{n-1}+h_{n-2})}$
is larger than the right hand side, and thus larger than the left, is
straightforward: after some algebra we have

$\alpha-\mathrm{rhs}={\frac{h_{{n-1}}h_{{n-2}}h_{{n}}\left(h_{{n-2}}+2\,h_{{n-1}}+h_{{n}}\right)}{\left(h_{{n-1}}+h_{{n}}\right)\left(h_{{n}}+h_{{n-1}}+h_{{n-2}}\right)A}}$
(32)

where

$A=2\,{h_{{n}}}^{2}+2\,h_{{n}}h_{{n-2}}+4\,h_{{n}}h_{{n-1}}+h_{{n-2}}h_{{n-1}}+2\,{h_{{n-1}}}^{2}$
(33)

which, being positive, completes the theorem. $\natural$

Figure 2: Condition numbers of the tridiagonal compact finite difference
matrix for random meshes. We sampled $10,000$ meshes of dimension $n=10,946$
(a Fibonacci number). The meshes were generated by h = rand(1,n). The vertical
line is at the dimension, $n$, suggesting that the mean condition number grows
like $n$. Similar graphs were generated for other $n$, up to $n=46,368$. Plot
generated by the command  histogram(
log(condsR(:,1))/log(10),’DisplayStyle’,’stairs’).

####  4.3.1 Some details of the equally-spaced case

In the equally-spaced case the $(1,1)$ and $(n+1,n+1)$ entries each become
$1/3$, the diagonal entries are $4$ otherwise, and all sub and superdiagonal
entries are $1$. LU factoring the matrix gives bidiagonal $\bm{\mathbf{L}}$
and $\bm{\mathbf{U}}$. The diagonal entries of $\bm{\mathbf{U}}$ are $1/3$,
$1$, $3$, $11/3$, $41/11$, and so on up until the penultimate entry, for
$2<k<n-2$ with the form $a_{k+1}/a_{k}$ where $a_{k+1}=4a_{k}-a_{k-1}$, and
initial conditions $a_{0}=a_{1}=1$ not reflected in the first diagonal entry.
This is easily proved by induction. The final entry in $\bm{\mathbf{U}}$, call
it $u_{n+1,n+1}$ must have $a_{n-1}/a_{n}\cdot 1+u_{n+1,n+1}=1/3$ (the
subdiagonal entries in $L$ are the reciprocals of the diagonal entries of
$\bm{\mathbf{U}}$) and so $u_{n+1,n+1}=1/3-a_{n-1}/a_{n}$. Solving the
recurrence relation for the $a_{k}$ we have

$a_{k}=\left(\frac{1}{2}-\frac{\sqrt{3}}{6}\right)\left(2+\sqrt{3}\right)^{k}+\left(\frac{1}{2}+\frac{\sqrt{3}}{6}\right)\left(2-\sqrt{3}\right)^{k}$
(34)

so $1/3-u_{n+1,n+1}$ very rapidly approaches
$2-\sqrt{3}=1/(2+\sqrt{3})=0.2677949$; indeed already $11/41=0.26829$ and
thereafter $a_{k}/a_{k+1}=1/(2+\sqrt{3})+O(c^{-k})$ where $c\approx 13.9$.
This is important because the final entry $u_{n+1,n+1}$ is then very close to
$1/3-1/(2+\sqrt{3})\approx 0.06538$.

These details can be used to establish the following, which I first observed
as an experimental fact:

######  Theorem 4.5.

The $1$-norm condition number of the $n+1$ by $n+1$ matrix $\bm{\mathbf{A}}$
for uniform mesh is asymptotically $63+36\sqrt{3}\approx 125.35$,
exponentially quickly independent of $n$.

######  Proof 4.6.

Omitted for brevity.

####  4.3.2 Yet another alternative treatment of the endpoints

In Zhao et al. (2007) we didn’t use equation (19). Instead we added an extra
node $\tau_{4}$, which requires $n>3$, in order to specify the ratio
$a_{0}/b_{0}=c:=2+\sqrt{3}$. This allows analytical factoring of the matrix,
speeding the process up slightly (it lowers the constant in the $O(n)$ cost).
In detail, we found the edge formula

$\displaystyle cf^{\prime}(0)+f^{\prime}(h)=-\frac{1}{h}\bigg{(}$
$\displaystyle\left(-25c/12-1/4\right)f\left(0\right)$
$\displaystyle+\left(4\,c-5/6\right)f\left(h\right)+\left(-3\,c+3/2\right)f\left(2\,h\right)$
$\displaystyle+\left(4c/3-1/2\right)f\left(3\,h\right)+\left(-c/4+1/12\right)f\left(4\,h\right){\bigg{)}}$
(35)

(and a similar formula at the right edge with the sign of $h$ reversed but
where $c=4$) which has truncation error

$\frac{(4c-1)h^{4}}{20}f^{(v)}(0)+\frac{(40c-11)h^{5}}{120}f^{(vi)}(0)+\cdots\>.$
(36)

This gives compact finite differences, again fourth order accurate. Note that
the truncation error above is larger than the truncation error for the four-
node formula; this seems to be compensated for by the fact that the condition
number is smaller, by a similar amount.

A similar formula is available for a nonuniform mesh. However, it does not
have the same speed advantage that it does for a uniform mesh, because the
matrix cannot (so far as I know) be analytically factored. The proof of
Theorem 4.1 is easily modified to show that the resulting matrix (for
nonuniform meshes) is also totally nonnegative for positive mesh widths; in
fact, this follows from a general fact about totally nonnegative matrices that
says that the $(1,1)$ entries and the $(n+1,n+1)$ entries can always be
increased without altering the total nonnegativity (see (Hogben, 2006, Chapter
29) for a statement and reference).

As mentioned, the advantage of this formula for a uniform mesh is that the
tridiagonal matrix $\bm{\mathbf{A}}$ can be analytically factored into
$\bm{\mathbf{L}}$ and $\bm{\mathbf{U}}$, where independent of dimension the
subdiagonal elements of $\bm{\mathbf{L}}$ are all $1/c=1/(2+\sqrt{3})$, and
the diagonal elements of $\bm{\mathbf{U}}$ are all $c$. This saves the $O(n)$
cost of factoring the matrix, both in terms of operation count and in terms of
storage. The relatively trivial disadvantage of this formula is that it needs
one more function value ($f(4h)$) on the right hand side. This makes no
practical difference because in the only case where this matters, $n=3$, one
can instead use a single cubic. For the uniform mesh formula, the $1$-norm
condition number is again asymptotically constant, but now being less than $3$
irrespective of the dimension; however, the truncation error at the edges is
larger, and thus this advantage does not seem to matter.

###  4.4 Nonuniform mesh compact derivatives/spline in practice

So, how well do these methods work in practice? Actually, pretty well. As
expected, if the mesh ratios are not “too large,” i.e., adjacent subintervals
are not too different in width (so that the $r_{k}$ factors in the edge
formulæ and the $r$ and $s$ factors in the interior formula are not too
large), then the formula is very similar to the cubic spline and its
derivative, for smooth functions. As the mesh ratios depart more seriously
from uniformity, the compact method begins to be more accurate than the cubic
spline.

Note, however, that the formulæ as initially written in terms of the
$\tau_{k}$ are quite likely to produce rounding errors, especially in the
formulæ at the edges, and should be rewritten using
$h_{i+1}=\tau_{i+1}-\tau_{i}$ and factored, as done in the six-by-six example
above, wherever possible. When this is done the influence of rounding errors,
while still felt, is significantly reduced. These formulæ have been
implemented in Matlab in a program called vcompact4, available on the code
repository for Corless & Fillion (2013). For this paper, I have written a
program called compactcubic, which I discuss later.

Figure 3:  Fourth-order nonuniform mesh compact finite difference maximum
error computing $(\nicefrac{{1}}{{\Gamma}})^{\prime}(x)$ on $1\leq x\leq 3$.
The spatial mesh was Chebyshev points $x_{j}=2+\cos(\nicefrac{{\pi j}}{{n}})$
for $0\leq j\leq n$, for various $n$. Theoretical fourth-order behaviour is
shown by the dashed line. Figure 4:  Error in fourth-order compact cubic
spline and its derivative for the Runge function on random meshes on $[-1,1]$.
The random meshes were generated by h = rand(1,n) and then using a cumulative
sum to generate grid points. The mesh widths varied in each individual mesh by
factors up to about $10^{5}$. Figure 5:  Error in fourth-order compact cubic
spline and its derivative for the Runge function on uniform meshes on
$[-1,1]$.

######  Remark 4.7.

The tridiagonal system for splines or compact finite differences is always
nonsingular for positive mesh widths, as we proved. Moreover, it is often
well-conditioned (as we will see in the next section). As always with finite-
differences, however, rounding errors will make their presence felt for small
enough interval widths. This is because differentiation is infinitely ill-
conditioned (a fact which is well-known, but rarely discussed as such: see the
exposition in chapter 11 of Corless & Fillion (2013) for one instance).

##  5 Explanation

Why are the equations for a cubic spline, which guarantee twice-continuous
differentiability for the piecewise cubic interpolation of the data, the same
equations as those for (equally-spaced) fourth-order compact finite
differences? The problems being solved by the two algorithms are different!

Saying it again: in the cubic spline case, we are constructing separate cubic
polynomials on each subinterval $\tau_{k}\leq t\leq\tau_{k+1}$. To guarantee
twice continuous differentiability, we specify that the interior derivatives
$\rho_{i,1}$ for $1\leq i\leq n-1$ should satisfy certain equations. But why
should those equations give us fourth-order accurate derivatives in the
equally-spaced case, but only third-order accurate derivatives for locally
nonuniform meshes?

The most satisfactory explanation I have found so far is that the partial
fraction decomposition of

$\frac{1}{(z-\tau_{k-1})^{2}(z-\tau_{k})^{3}(z-\tau_{k+1})^{2}}$ (37)

has zero coefficient333This is similar to the contour integral explanation of
the well-known extra order of accuracy of Simpson’s rule: we fit three points,
but in the equally-spaced case we get exact fit for a cubic, not just a
quadratic. of $1/(z-\tau_{k})^{2}$ if (and only if)
$\tau_{k+1}-\tau_{k}=\tau_{k}-\tau_{k-1}$. Explicitly, here, if
$\tau_{k+1}-\tau_{k}=\tau_{k}-\tau_{k-1}=h$ (as we will see without loss of
generality we may take $h=1$) we have the following

######  Theorem 5.1.

Let $t$ and $h$ be fixed complex numbers, $h\neq 0$, and $p(z)$ be a
polynomial in the variable $z$. If the polynomial $p(z)$ has degree at most
$5$, then the values of $p(t-h)$, $p^{\prime}(t-h)$, $p(t)$, $p(t+h)$, and
$p^{\prime}(t+h)$ determine $p^{\prime\prime}(t)$ uniquely. Simultaneously, if
$p(z)$ has degree less than or equal to $4$, those same values determine
$p^{\prime}(t)$ uniquely.

######  Proof 5.2.

For intelligibility we may without loss of generality replace $z$ with
$\theta$ where $z=t+\theta h$, so $-1\leq\theta\leq 1$. Then we have the
partial fraction decomposition

$\frac{1}{(\theta+1)^{2}\theta^{3}(\theta-1)^{2}}=-\frac{\nicefrac{{1}}{{4}}}{(\theta+1)^{2}}-\frac{1}{\theta+1}+\frac{1}{\theta^{3}}+\frac{2}{\theta}+\frac{1}{\theta-1}+\frac{\nicefrac{{1}}{{4}}}{(\theta-1)^{2}}\>,$
(38)

which (as noted in the preamble to the theorem) has no $1/\theta^{2}$ term. By
a standard argument, the contour integral

$\frac{1}{2\pi
i}\oint_{C}\frac{p(\theta)}{(\theta+1)^{2}\theta^{3}(\theta-1)^{2}}\,d\theta=0$
(39)

for all polynomials of degree $5$ or less if the contour $C$ contains $-1$,
$0$, and $1$. Therefore the partial fraction decomposition and Cauchy’s
integral formula gives

$0=-p(-1)+2p(0)-p(1)+\frac{1}{4}\left(p^{\prime}(1)-p^{\prime}(-1)\right)+\frac{1}{2}p^{\prime\prime}(0)\>.$
(40)

This means that if we know the value of the function at three consecutive
points (which for a spline would take two cubic polynomials to fit), with
equal width intervals on either side of the midpoint $0$, and we know the
first derivatives at the two endpoints, then the value of the second
derivative at the midpoint $0$ is determined irrespective of the value of the
first derivative $p^{\prime}(0)$. Since these same quantities determine
$p^{\prime}(0)$ by equation (11) for polynomials of degree at most four, the
proof is complete. $\natural$

Undoing the nondimensionalization and translating this into a finite-
difference formula on the three points $\tau_{k-1}$, $\tau_{k}$, $\tau_{k+1}$
where $h=\tau_{k+1}-\tau_{k}=\tau_{k}-\tau_{k+1}$, we have

$f^{\prime\prime}(\tau_{k})=2\frac{f(\tau_{k+1})-2f(\tau_{k})+f(\tau_{k-1})}{h^{2}}-\frac{f^{\prime}(\tau_{k+1})-f^{\prime}(\tau_{k-1})}{2h}+\frac{1}{360}f^{(6)}(0)h^{4}+O(h^{8})\>.$
(41)

Notice the familiar second difference formula and the centered difference
formula (using $f^{\prime}$) appearing. That this combination is fourth-order
accurate was new to me but hardly seems surprising; doubtless the formula is
in one book of finite differences or another, perhaps Milne-Thomson (1951),
Collatz (1966), Jordán (1965), or even Boole (1960). The error term in
equation (41) was computed in Maple, and indicates that the formula is exact
if the sixth derivative is identically zero, i.e. if $f(z)$ is a polynomial of
degree at most five. Since the value of the derivative at $\tau_{k}$ is not
used, we can see that imposing a fourth-order accurate value of the derivative
there must respect this constraint on the second derivative.

An alternative view is that in the equally-spaced case we may force an
approximate value of $f^{\prime}(\tau_{k})$ that is in $O(1)$ error without
altering this finite-difference value of $f^{\prime\prime}(\tau_{k})$.

That this symmetry forces the value of the second derivative to be a
particular value for all fourth-order accurate derivatives still seems
somewhat surprising. The key seems to be the extra accuracy, allowing the
second derivative to be determined, that is permitted by the zero residue in
the (locally) equally-spaced case.

##  6 Numerical tests

I wrote a small didactic (and therefore relatively inefficient) Matlab program
that I called compactcubic, based on the codes pchiptx and splinetx from Moler
(2004), and which I will make available, either in the code repository for
Corless & Fillion (2013) or some other convenient place.

Rather than leave the interpolant in the Hermite interpolational basis, which
I prefer, for ease of comparison I computed the local monomial basis
coefficients from the derivatives, much as is done in splinetx. Indeed the
formulæ are very similar. I added the ability to output the first derivative.
I also examine the $2$nd derivative, which turns out to be informative, though
normally clearly departing from the $2$nd derivative of the underlying
function. That this is not continuous for nonuniform meshes is demonstrated in
figure 8.

I ran this on several functions, notably the Runge example
$f(x)=1/(1+25x^{2})$ and the signum function (sign in Matlab), using uniform
meshes and nonuniform meshes of random widths. The nonuniform meshes of random
widths were constructed by first generating the random widths with rand(1,n)
and then constructing the nodes by use of cumsum.

Some of the results are given in figures 4 and 5. In these tests, the programs
appear to have behaved satisfactorily. The condition numbers of the
tridiagonal matrices depend on the node family used. Equally-spaced nodes
generate a condition number essentially constant; random meshes generate
condition numbers that depend on the mesh spacing, but not the dimension. The
$O(1/n^{4})$ behaviour of the errors on uniform grids and the $O(1/n^{3})$
behaviour of the errors of the derivative are shown.

Figure 6:  Plot of the error in the fourth-order uniform mesh ($n=64$) compact
finite difference derivative, on the Runge example $f(x)=1/(1+25x^{2})$, on
the interval $[-1,1]$. Solid line is the error between mesh points. The error
is largest in the middle, and believably $O(h^{3})$. The error in the
interpolant itself (not the derivative as we have it in this figure) is
$O(h^{4})$.  Figure 7:  Plot of the error in second derivative of the cubic
Hermite interpolant. We use $n=7$ subintervals on $[-1,1]$ for the Runge
example $f(x)=1/(1+25x^{2})$. The graph shows clearly that the second
derivatives of the spline are continuous at the nodes, although as expected
they do not match the second derivatives at the nodes computed by the repeated
application of vcompact4.m. Figure 8:  Plot of the error in second derivative
of the compact cubic interpolant. We use $n=7$ subintervals with nonuniform
(in fact, Chebyshev) nodes on $[-1,1]$ for the Runge example
$f(x)=1/(1+25x^{2})$. The graph shows clearly that the (linear) second
derivatives of the cubic spline are not continuous at the nodes, as they would
have been for equally-spaced nodes. Circles are the exact second derivatives
at the nodes. Accuracy is about as expected with only seven nodes and errors
are clearly visible.

######  Remark 6.1.

Somewhat annoyingly, the plain vanilla Matlab built-in spline function and
piecewise polynomial evaluation routines do not have automatic access to
evaluating the derivatives of the piecewise functions. Of course one could
write one’s own spline differentiator using unmkpp and mkpp (it only takes two
lines of code), but bookkeeping other people’s choices of piecewise polynomial
representation is not always trivial, and the user might not notice the
solution given in the help page for mkpp. I believe that the built-in ability
to take derivatives is always helpful444As a colleague points out, piecewise
polynomial differentiation functionality is available in the CurveFitting
package, which can be downloaded and installed (at my institution, at no extra
cost; on the Mathworks page it offers the possibility of a free trial).. But
for this paper, to make comparisons, instead of using mkpp and unmkpp I
modified the splinetx.m function to return derivatives.

As predicted, for equally-spaced nodes the derivatives were the same. For
Chebyshev nodes the derivatives were nearly indistinguishable ($r-s$ is
typically quite small); for Fibonacci nodes there was a difference, but not
much.

###  6.1 Contour Integrals and Higher Derivatives

The contour integral technique used above is of course classical. I learned it
from John C. Butcher, who used the following elegant variation in Butcher
(1967) and later in Butcher et al. (2010) to compute higher derivatives. The
idea is to multiply by another factor, which I call the Butcher factor and
denote by $B(z)$, that zeros out the unneeded and interfering residues. As an
example, suppose that we wish to find a formula for the second derivative of a
cubic spline. We then wish a rational function with denominator
$(z-\tau_{n})^{2}(z-t)^{3}(z-\tau_{n+1})^{3}$, so that we will pick up a term
of the form $1/(z-t)^{3}$ which will supply us with $p^{\prime\prime}(t)/2!$
by the Cauchy Integral Formula. But this approach will also generate a term
for $p(t)$ and another for $p^{\prime}(t)$, all mixed in together; so we
introduce the Butcher factor

$B(z)=b_{0}+b_{1}z+b_{2}z^{2}$ (42)

and choose the coefficients $b_{0}$, $b_{1}$, and $b_{2}$ (not all zero) so
that the residues of $1/(z-t)$ and $1/(z-t)^{2}$ are zero. This gives us two
linear equations in three variables to solve, which we can do
straightforwardly. One reasonably presentable solution is to put
$t=\tau_{n}+h\phi$, $z=\tau_{n}+h\theta$ and $\tau_{n+1}=\tau_{n}+h$ for
neatness, and then

$B(\theta)=\left(\left(3\,\phi-2\right){\phi}^{3}\right)-\left(2\,{\phi}^{2}\left(4\,\phi-3\right)\right)\,\theta+\left(6\,{\phi}^{2}-6\,\phi+1\right)\,{\theta}^{2}$
(43)

which gives in turn

$\frac{B(\theta)}{\theta^{2}(\theta-\phi)^{3}(\theta-1)^{2}}={\frac{-3\,\phi+1}{\left(\theta-1\right)^{2}}}+{\frac{-6\,\phi+3}{\theta}}+{\frac{6\,\phi-3}{\theta-1}}+{\frac{-3\,\phi+2}{{\theta}^{2}}}+\frac{1}{\left(\theta-\phi\right)^{3}}\>.$
(44)

Notice that as intended the coefficients of $1/(z-t)$ and $1/(z-t)^{2}$
(equivalently $1/(\theta-\phi)$ and $1/(\theta-\phi)^{2}$) are zero. This in
turn gives (after the same contour argument as before) that

$\frac{2-3\phi}{h}p^{\prime}(\tau_{n})+\frac{1-3\phi}{h}p^{\prime}(\tau_{n+1})+\frac{p^{\prime\prime}(t)}{2!}+\frac{6\phi-3}{h^{2}}p(\tau_{n+1})+\frac{3-6\phi}{h}p(\tau_{n})=0\>,$
(45)

if the degree of $p(z)$ is at most three. Rearranging, we get

$p^{\prime\prime}(\tau_{n}+\phi
h)=\frac{2(3\phi-2)}{h}p^{\prime}(\tau_{n})+\frac{2(3\phi-1)}{h}p^{\prime}(t_{n+1})+\frac{6(1-2\phi)}{h}\left(\frac{p(t_{n+1})-p(\tau_{n})}{h}\right)\>.$
(46)

To use the same idea to find a compact formula for the second derivative given
only function values on equally-spaced nodes, the idea is to search for a
Butcher factor $B(z)$ that makes the coefficients of $(\theta+1)^{-2}$,
$\theta^{-2}$, and $(\theta-1)^{-2}$ zero in the partial fraction expansion of

$\frac{b_{0}+b_{1}\theta+b_{2}\theta^{2}}{(\theta+1)^{3}\theta^{3}(\theta-1)^{3}}\>.$
(47)

Here we have taken advantage of the equal spacing of the mesh to put
$z=\tau_{k}+\theta h$ again. A short computation shows that
$B(\theta)=9\theta^{2}-5$ gives

$\frac{9\theta^{2}-5}{(\theta+1)^{3}\theta^{3}(\theta-1)^{3}}=\frac{-3}{\theta+1}+\frac{6}{\theta}-\frac{3}{\theta-1}+\frac{\nicefrac{{1}}{{2}}}{(\theta+1)^{3}}+\frac{5}{\theta^{3}}+\frac{\nicefrac{{1}}{{2}}}{(\theta-1)^{3}}\>.$
(48)

After using the Cauchy Integral Formula again, we find the following fourth-
order formula

$p^{\prime\prime}(\tau_{n}-h)+10p^{\prime\prime}(\tau_{n})+p^{\prime\prime}(\tau_{n}+h)=\frac{12}{h^{2}}\left(p(\tau_{n}-h)-2p(\tau_{n})+p(\tau_{n}+h)\right)\>.$
(49)

The Cauchy formula already tells us that this is exact for polynomials of
degree at most $5$. A separate Taylor series expansion (trivial in a computer
algebra system) gives that the error term for functions not polynomials of
degree at most five is asymptotically, as $h\to 0$,

$\frac{h^{4}}{20}p^{(6)}(\tau_{n})+\frac{11h^{6}}{5040}p^{(8)}(\tau_{n})+O(h^{8})\>.$
(50)

###  6.2 Notes and further reading

The books de Boor (1978) and Jordán (1965) are each rich sources of the
history of their respective subjects. Much of this section comes from them.
The word “spline” comes from the same root as the English word “splinter”,
meaning a thin fragment of wood; such wooden pieces gave rise to a flexible
mechanical device that was used, apparently in shipbuilding, to draw smooth
curves between fixed points. I. J. Schoenberg in Schoenberg (1946) chose the
name “spline” for the piecewise cubic polynomials computed in this fashion
because they minimize an approximation of the strain energy that the original
flexible instrument minimized. The modern terminology, and in particular the
difference between a cubic spline, a pchip, and a true cubic Hermite
interpolant (where the derivative values $f^{\prime}(\tau_{k})$ are also known
exactly and not chosen by the algorithm), is somewhat confusing, sometimes
even to experts.

I was unable to trace the first use of the word “compact” for finite
differences that are defined as a relation amongst function and derivative
values on a small set of mesh points. Some authors use instead “optimal”,
conveying that the maximum accuracy possible is attained given the constraints
on the number of mesh points. Collatz used the term “Mehrstellenverfahren”,
which is still used.

This paper is concerned with the connection between these two concepts. To
explain this connection I used the contour integral approach. As is very well
known, this is not the only approach to finite differences, compact or
otherwise.

One standard approach is simply to use brute Taylor expansion of an ansatz
with undetermined coefficients, then set all the desired coefficients of
powers of the step size to zero, and brutally solve for the unknowns. Doing
this once and for all for a fixed (possibly symbolic) grid in a computer
algebra system is perfectly straightforward and perfectly useful, and I have
used this brutal method in a few papers such as Zhao et al. (2007), and before
that in Corless & Rokicki (1996). As pointed out in Fornberg (1998) this
approach cannot easily be used in a purely numerical environment because the
linear systems are often quite ill-conditioned, and more so for larger
formulæ; but in a computer algebra system the arithmetic can be done exactly
(if rational), or at whatever precision desired, and the expense of doing so
is amortized over the life of the use of the formulæ.

Also as pointed out in Fornberg (1998), and in Corless & Fillion (2013), the
use of the formulæ (however obtained) is subject to numerical difficulty; this
is because differentiation is infinitely ill-conditioned. But the brute-force
approach of generating formulæ by computer algebra systems works and gets
formulæ as accurate as any method does555Well, maybe there is still room for
blunders. The paper Keller & Pereyra (1978) used Macsyma to generate weights
for finite difference formulæ. The paper Fornberg (1988) points out that the
tables in Keller & Pereyra (1978) contained “both isolated and systematic
errors”. Certainly while I have been as careful as I can for this paper, and I
have taken out a lot of typos, I have probably introduced some in transferring
formulæ from Maple to LaTeX..

I believe that the contour integral approach is better at explanation,
however. Here, of course, we are trying to understand the connection of cubic
splines to compact finite differences, and the contour integral approach
supplied a satisfactory explanation. Contour methods are often used in this
way. According to Weideman (2005) contour integrals were used by the great
number theorist K. Mahler to explore Padé approximations to the logarithms of
algebraic numbers in Mahler (1953). The connection to Padé approximation (and
indeed compact finite difference formulæ are often called “Padé methods”)
arises because of another, perhaps even more elegant method for constructing
finite difference formulæ, the operator method. One version of the operator
method, derived by applying difference operations to the exponential function,
is used to great effect in Fornberg (1998) to compute weights in many finite-
difference formulæ and to give general recurrence relations for doing so on
arbitrary grids. In older works, such as Fox (1957), we see a more “formal”
treatment of operator methods; here ‘formal’ means that rigour is not
maintained throughout the computation (a curious usage of this word, that
often causes puzzlement, meaning only that computations are just carried out
by ‘form’ and not worrying until later whether they were correct). The
operator method was recently blogged about in Scientific American by Lamb
(2019). I gave a talk on this method at a CAIMS meeting in Winnipeg in 2005.
It seems to be common in early works on finite differences, such as Milne-
Thomson (1951), and is at least a hundred and sixty years old, being present
in Boole (1960); in fact, according to Jordán (1965) the method is due to
Lagrange himself.

The operator method rests on an analogy: the shift operator
$E_{h}(f)(t)=f(t+h)$ can be expressed in terms of the derivative operator $hD$
by Taylor’s theorem for analytic functions as

$E_{h}(f)(t)=f(t+h)=\sum_{k\geq 0}\frac{h^{k}}{k!}D^{(k)}(f)(t),$ (51)

where the symbol $D^{(k)}$ means repeated application of the differentiation
operator. By analogy with multiplication, this is written as
$E_{h}=\exp(hD)(f)(t)$, which though likely familiar to the reader from
various contexts is actually a remarkable leap to a definition of an
exponential of an operator. Formally solving for $D$ gives
$D=\ln(E_{h})(f)(t)/h=\ln(1+\Delta)(f)(t)/h$, motivating Padé approximations
to the logarithm near $1$ (or the inverse sinh function near $0$ because
analogously $\exp(hD/2)-\exp(-hD/2)=2\sinh(hD/2)$ applied to $f(t)$ gives
$f(t+h/2)-f(t-h/2)$ and so the inverse sinh function expresses $D$ in terms of
compositions of central differences). The resulting Padé approximations then
give us compact finite difference formulæ.

The paper Weideman (2005) extends this idea to Hermite-Padé approximations,
and solves several general problems associated with these. Further, that paper
details an open problem about explicit descriptions of arbitrary width compact
formulæ for the second derivative; I do not know if the Butcher factor idea
for contour integration can be used to attack that open problem, but John C.
Butcher did use his approach to solve another general degree interpolation
problem in quadrature (we detailed his solution in Butcher et al. (2010)), so
it might be possible. I note that one can establish that a necessary condition
for the unwanted residues to vanish if the degree of $B(z)$ is to be no more
than $n$ is that at each node $\tau_{k}$, $0\leq k\leq n$, the Butcher factor
must have

$B^{\prime}(\tau_{k})-3\left(\sum_{\begin{subarray}{c}j=0\\\ j\neq
k\end{subarray}}^{n}\frac{1}{\tau_{k}-\tau_{j}}\right)B(\tau_{k})=0\>,$ (52)

which implies that $B(z)$ can be identified in the Lagrange basis by its
values on the nodes by finding a vector in the nullspace of
$\bm{\mathbf{D}}-3\mathrm{diag}(s_{0},s_{1},\ldots s_{n})$ where
$s_{k}=\sum_{\begin{subarray}{c}j=0\\\ j\neq
k\end{subarray}}^{n}\frac{1}{\tau_{k}-\tau_{j}}$ and $\bm{\mathbf{D}}$ is the
differentiation matrix on the nodes (for differentiation matrices, see e.g.
Amiraslani et al. (2018)). I leave this to future work.

###  6.3 Higher dimension

Compact finite differences are widely used in problems having more than one
dimension. An example on a regular grid is Rokicki & Floryan (1995), and on
irregular grids is Rokicki & Floryan (1999). Automatic computation of such
formulæ for scattered nodes in more than one dimension are considered in more
detail in the University of Western Ontario Ph.D. thesis of Jichao Zhao (2006)
and a subsequent paper Zhao & Corless (2006) and a Maple Share Library package
released in that same year updating the original 1994 package. Other
researchers have also studied the problem, for example Wright & Fornberg
(2006) who use radial basis functions and were likely the first to do this in
full generality.

##  7 Concluding remarks, and a marriage of convenience

The usual methods of dealing with the two extra degrees of freedom of a cubic
spline have always seemed somewhat unnatural to me. As has been demonstrated
with many psychological studies, the extra choice seems unwelcome. What I
recommend here is to use compact finite differences to not only compute the
derivatives at the interior nodes, but also at the endpoints. This is very
similar in spirit to the recommendation of Swartz & Varga (1972), except here
we are using a compact finite difference formula to give fourth-order
accuracy, where they recommended a third-order accurate direct formula. The
method recommended here requires essentially no more effort, being simply a
modification of the tridiagonal system needed to find the spline derivatives,
and has higher order accuracy for the derivatives at the edges. This new
compact method, with its positive definite (totally nonnegative) matrices,
therefore seems as though it will be as good a method as any, and might be
psychologically more satisfying in that the derivatives at the end are
determined by the data. Of course, they are also determined by the data, in a
sense, if instead the “not-a-knot” condition is used, but that makes a
qualitative difference with the two subintervals at each end. Whether it makes
a practical difference is another matter.

I have not made extensive comparisons with Matlab ’s built-in spline, for
example (but the few experiments I have done show that the behaviour is at
least similar). The fact that the Not-A-Knot condition generates an $O(h^{3})$
error in the derivative at the end may or may not be significant. Nor have I
made a choice between the $2+\sqrt{3}$ method—which needs at least $n=4$ so
therefore $5$ nodes and has an analytic factoring for uniform meshes—and the
$h_{2}(h_{2}+h_{3})/((h_{1}+h_{2})(h_{1}+h_{2}+h_{3}))$ method, which also
needs at least $5$ nodes because its four-by-four matrix is singular; in the
first case the condition number of the matrix is smaller, but the truncation
error is larger (by a similar factor). My only experience with this in
practice is with the uniform mesh case, where the efficiency of the analytic
factoring is noticeably useful. More experimentation is necessary to learn if
the variable mesh formula is truly useful.

The method has been coded in Matlab as the program compactcubic.m and will be
made available at the code repository for Corless & Fillion (2013), namely
nfillion.com/coderepository.

It seems that this approach gives a potentially reasonable alternative to
splines, in that its derivative can be expected to be more accurate at the
nodes for nonuniform meshes, and one does not have to make a choice about what
to do at the endpoints. That first caveat, “at the nodes” is important. The
error in piecewise cubic interpolation means that the derivatives between the
nodes can only be $O(h^{3})$ accurate. This weakens the case for the compact
cubic interpolant proposed here, it’s true; away from the nodes, one doesn’t
see more accurate derivatives.

This approach needs a name to distinguish it from “spline”, “cubic Hermite”,
and “pchip” (itself a misleading name for the interpolant whose derivatives
are chosen to preserve qualitative features, at the cost of even less accuracy
in the derivative). I propose the name “compact cubic interpolant” for this
construction. It also needs an efficient and stable implementation; probably
conversion to the local monomial basis is the most efficient, as is done here,
but there are likely other tricks needed to ensure that the computation is as
stable as possible.

I realize that I have not made a strong case for using this method—cubic
splines are pretty good, after all, even if they occasionally wiggle too much
(which is what pchip is for)—but the exploration of this alternative has
yielded some advantages for the compact cubic method: file-it-and-forget-it
treatment of the edges, and in theory more accurate derivatives at the nodes.
If one is using compact finite differences for some other purpose anyway
(perhaps as part of solving a PDE by the method of lines) then the analysis in
this paper may help you to choose to use the compact cubic interpolant.

I remark that totally positive (nonnegative) matrices occur in various
interpolation problems—see Marco et al. (2017) for an example.

The other result of this paper, namely an explanation, using contour integrals
and residues, of the apparent coincidence that cubic splines give the same
derivatives that the fourth order Padé compact formula does, is only one such
explanation. This coincidence can also be explained by noting first (as can be
found in Appendix A of Fornberg & Flyer (2015) for instance) that a cubic
spline over $x_{0}-h,x_{0},x_{0}+h$ is the sum of a pure cubic and a multiple
of $|x-x_{0}|^{3}$. Then since the compact formula is exact for pure cubics,
one only needs to think about its behaviour on $|x-x_{0}|^{3}$; but this is an
even function, and a moment’s reflection shows that (9) is exact for even
functions, giving zero derivative at the centre; so since the equation is
exact for each, it is exact for their sum, and hence for splines. Thus this
coincidence might not be so surprising for one “skilled in the art”. This
argument works only in the case of real nodes, of course, whereas the residue
argument works over $\mathbb{C}$, but the overwhelming majority of spline
interpolation takes place over $\mathbb{R}$, so that doesn’t matter much.
Still, I think that this coincidence has hardly been noticed in the
literature; that it can be explained several ways is perhaps not a surprise. I
do like the residue explanation in preference to the $|x-x_{0}|^{3}$
explanation because zero residues are connected to the extra degree of
accuracy obtained in Simpson’s Rule, and connected to the Butcher method of
finding higher derivative formulæ and solving the Birkhoff interpolation
problem. This, however, is clearly a matter of taste: I am certain that there
are many who would prefer the explanation using $|x-x_{0}|^{3}$.

A referee points out that I have not reviewed the spline literature thoroughly
in this paper, omitting reference even to the classic book Schumaker (2007).
In my defence I point out that this paper is already too long.

## References

  * Amiraslani et al. (2018) Amiraslani, A., Corless, R. M. & Gunasingam, M. (2018)  Differentiation matrices for univariate polynomials.  Numerical Algorithms, 1–31.
  * Boole (1960) Boole, G. (1960)  A treatise on the calculus of finite differences.  Dover.  (Orig. pub 1860. 2nd edition 1872 edited by J. F. Moulton).
  * Butcher (1967) Butcher, J. C. (1967)  A multistep generalization of Runge-Kutta methods with four or five stages.  J. ACM, 14, 84–99.
  * Butcher et al. (2010) Butcher, J. C., Corless, R. M., Gonzalez-Vega, L. & Shakoori, A. (2010)  Polynomial algebra for Birkhoff interpolants.  Numerical Algorithms, 1–29.
  * Collatz (1966) Collatz, L. (1966)  The numerical treatment of differential equations, 3rd edn.  Springer-Verlag.
  * Corless & Fillion (2013) Corless, R. M. & Fillion, N. (2013)  A graduate introduction to numerical methods, vol. 10.  Springer, p. 12.
  * Corless & Rokicki (1996) Corless, R. M. & Rokicki, J. (1996)  The symbolic generation of finite-difference formulae.  Z.A.M.M., 76, 381–382.
  * de Boor (1978) de Boor, C. (1978)  A practical guide to splines.  Springer-Verlag.
  * Fornberg (1988) Fornberg, B. (1988)  Generation of finite difference formulas on arbitrarily spaced grids.  Mathematics of computation, 51, 699–706.
  * Fornberg (1998) Fornberg, B. (1998)  Classroom note: Calculation of weights in finite difference formulas.  SIAM Review, 40, 685–691.
  * Fornberg & Flyer (2015) Fornberg, B. & Flyer, N. (2015)  A primer on radial basis functions with applications to the geosciences.  SIAM.
  * Fox (1957) Fox, L. (1957)  The numerical solution of two-point boundary problems in ordinary differential equations.  Oxford Clarendon.
  * Higham (2002) Higham, N. J. (2002)  Accuracy and Stability of Numerical Algorithms, 2nd edn.  Philadelphia: SIAM.
  * Hogben (2006) Hogben, L. (ed.) (2006)  Handbook of Linear Algebra.  Chapman and Hall/CRC.
  * Jordán (1965) Jordán, K. (1965)  Calculus of finite differences, vol. 33.  American Mathematical Soc.
  * Keller & Pereyra (1978) Keller, H. & Pereyra, V. (1978)  Symbolic generation of finite difference formulas.  Mathematics of Computation, 32, 955–971.
  * Lamb (2019) Lamb, E. (2019)  The theorem that unites different kinds of calculus.  Scientific American (blogs).
  * Lele (1992) Lele, S. K. (1992)  Compact finite difference schemes with spectral-like resolution.  Journal of Computational Physics, 103, 16–42.
  * Mahler (1953) Mahler, K. (1953)  On the approximation of logarithms of algebraic numbers.  Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences, 245, 371–398.
  * Marco et al. (2017) Marco, A., Martínez, J.-J. & Peña, J. M. (2017)  Accurate bidiagonal decomposition of totally positive Cauchy–Vandermonde matrices and applications.  Linear Algebra and its Applications, 517, 63–84.
  * Milne-Thomson (1951) Milne-Thomson, L. M. (1951)  The Calculus of Finite Differences, 2nd edn.  MacMillan and Company.  (1st edition 1933).
  * Moler (2004) Moler, C. (2004)  Numerical Computing with MATLAB, Electronic edn.  Philadelphia: SIAM.
  * Rokicki & Floryan (1995) Rokicki, J. & Floryan, J. (1995)  Domain decomposition and the compact fourth-order algorithm for the Navier-Stokes equations.  Journal of Computational Physics, 116, 79–96.
  * Rokicki & Floryan (1999) Rokicki, J. & Floryan, J. (1999)  Higher-order unstructured domain decomposition method for Navier–Stokes equations.  Computers & fluids, 28, 87–120.
  * Schoenberg (1946) Schoenberg, I. J. (1946)  Contributions to the problem of approximation of equidistant data by analytic functions.  Quart. Appl. Math., 4, 45–99 and 112–141.
  * Schumaker (2007) Schumaker, L. (2007)  Spline functions: basic theory.  Cambridge University Press.
  * Swartz & Varga (1972) Swartz, B. K. & Varga, R. S. (1972)  Error bounds for spline and L-spline interpolation.  Journal Of Approximation Theory, 6, 6–49.
  * Weideman (2005) Weideman, J. A. C. (2005)  Padé approximations to the logarithm I: Derivation via differential equations.  Quaestiones Mathematicae, 28, 375–390.
  * Weideman & Reddy (2000) Weideman, J. A. C. & Reddy, S. C. (2000)  A matlab differentiation matrix suite.  ACM Trans. Math. Softw., 26, 465–519.
  * Wright & Fornberg (2006) Wright, G. B. & Fornberg, B. (2006)  Scattered node compact finite difference-type formulas generated from radial basis functions.  Journal of Computational Physics, 212, 99–123.
  * Zhao et al. (2007) Zhao, J., Davison, M. & Corless, R. M. (2007)  Compact finite difference method for American option pricing.  J. Comput. Appl. Math., 206, 306–321.
  * Zhao & Corless (2006) Zhao, J. & Corless, R. M. (2006)  Compact finite difference method for integro-differential equations.  Applied Mathematics and Computation, 177, 271 – 288.