% Except otherwise noted,
% this work is © 2009 Dmitri Fedorov, under CC-BY-SA license, version 3+.

\chapter{Interpolation}

\section{Introduction}
In practice one often meets a situation where the function of interest,
$f(x)$, is only represented by a discrete set of tabulated points,
$$\{x_i, y_i=f(x_i) \mid i=1\dots n\},$$ obtained, for example, by sampling,
experimentation, or extensive numerical calculations.

{\em Interpolation} means constructing a (smooth) function, called {\em
interpolating function} or {\em interpolant}\index{interpolant}, which
passes exactly through the given points and hopefully approximates the
tabulated function in between the tabulated points.  Interpolation is a
specific case of {\em curve fitting}\index{curve fitting} in which the
fitting function must go exactly through the data points.

The interpolating function can be used for different practical needs
like estimating the tabulated function between the tabulated points and
estimating the derivatives and integrals involving the tabulated function.

\section{Polynomial interpolation}

Polynomial interpolation uses a polynomial as the interpolating function.
Given a table of $n$ points, $\{x_i,y_i\}$, where no two $x_i$ are the
same, one can construct a polynomial $P(x)$ of the order $n-1$ which
passes exactly through the points: $P(x_i)=y_i$.  This polynomial can
be intuitively written in the {\em Lagrange form}\index{Lagrange
interpolating polynomial},
	\begin{equation}\label{}
P(x)=\sum_{i=1}^{n}y_i\prod_{k\ne i}^{n}\frac{x-x_k}{x_i-x_k}\;.
	\end{equation}
The interpolating polynomial always exists and is unique.

\begin{table}[htb]\
\caption{Polynomial interpolation in C}
%\lstset{frame=single,tabsize=2,language=C,basicstyle=\footnotesize}
\lstinputlisting{c/interpolation/polinterp.c}
\end{table}

Higher order interpolating polynomials are susceptible to the {\em Runge's
phenomenon}\index{Runge's phenomenon}: erratic oscillations close to
the end-points of the interval (see Figure~\ref{fig:runge}).

\begin{figure}[htb]
\centerline{\includegraphics{runge.eps}}
\caption{Lagrange interpolating polynomial, solid line, showing the
Runge's phenomenon: large oscillations at the edges. For comparison
the dashed line shows a cubic spline. }
	\label{fig:runge}
\end{figure}

\section{Spline interpolation}
Spline interpolation uses a {\em piecewise polynomial}, $S(x)$, called {\em
spline}\index{spline}, as the interpolating function,
	\begin{equation}\label{eq:spl-spl}
S(x)=S_i(x) \:,\:\mathrm{if}\; x \in [x_i,x_{i+1}]\:,
\: i=1,\dots,n-1
	\end{equation}
where $S_i(x)$ is a polynomial of a given order $k$.

The spline of the order $k\ge1$ can be made continuous at the tabulated
points,
	\begin{equation}\label{eq:spl-s}
S_i(x_i)=y_i \:,\: S_i(x_{i+1})=y_{i+1}\:,\: i=1,\dots,n-1\;,
	\end{equation}
together with its $k-1$ derivatives,
	\begin{equation}\label{eq:spl-sp}
\left. \begin{array}{c}
S^\prime_i(x_{i+1})=S^\prime_{i+1}(x_{i+1})\:,\\
S^{\prime\prime}_i(x_{i+1})=S^{\prime\prime}_{i+1}(x_{i+1})\:,\\
\vdots \\
S^{(k-1)}_i(x_{i+1})=S^{(k-1)}_{i+1}(x_{i+1})\:.
	\end{array} \right| {i=1,\dots,n-2}
	\end{equation}

Continuity conditions (\ref{eq:spl-s}) and (\ref{eq:spl-sp}) make
$kn+n-2k$ linear equations for the $(n-1)(k+1)=kn+n-k-1$ coefficients of
$n-1$ polynomials (\ref{eq:spl-spl}) of the order $k$. The missing $k-1$
conditions can be chosen (reasonably) arbitrarily.

The most popular is the cubic spline, where the polynomials $S_i(x)$ are
of third order. The cubic spline is a continuous function together with
its first and second derivatives. The cubic spline has a nice feature
that it (sort of) minimizes the total curvature of the interpolating
function. This makes the cubic splines look good.

Quadratic splines---continuous with the first derivative---are not nearly
as good as cubic splines in most respects.  In particular they might
oscillate unpleasantly when a quick change in the tabulated function
is followed by a period where the function is nearly a constant.
Cubic splines are somewhat less susceptible to such oscillations.

Linear spline is simply a {\em polygon} drawn through the tabulated
points.

\subsection{Linear interpolation}
{\em Linear interpolation} is a spline with linear polynomials.
The continuity conditions (\ref{eq:spl-s}) can be satisfied by choosing
the spline in the (intuitive) form
	\begin{equation}\label{eq:spl-ls}
S_i(x)=y_i+p_i (x-x_i)\;,
	\end{equation}
where
	\begin{equation}
p_i = {\Delta y_i\over \Delta x_i}\;,\;\; \Delta y_i\doteq
y_{i+1}-y_i\;,\;\; \Delta x_i\doteq x_{i+1}-x_i\;.
	\end{equation}

\begin{table}[htb]
\caption{Linear interpolation in C}
%\lstset{frame=single,tabsize=2,language=C,basicstyle=\footnotesize}
\lstinputlisting{c/interpolation/linterp.c}
\end{table}

Note that the search of the interval $[x_i\le x\le x_{i+1}]$ in an ordered
array $\{x_i\}$ should be done with the {\em binary search}\index{binary
search} algoritm (also called {\em half-interval search}\index{half-interval
search}): the point $x$ is compared to the middle element of the array,
if it is less than the middle element, the algorithm repeats its action
on the sub-array to the left of the middle element, if it is greater,
on the sub-array to the right.  When the remaining sub-array is reduced
to two elements, the interval is found.  The average number of operations
for a binary search is $O(\log n)$.

\subsection{Quadratic spline}
Quadratic spline is made of second order polynomials,
conveniently written in the form
	\begin{equation}\label{eq:spl-qs}
S_i(x)=y_i+p_i(x-x_i)+c_i(x-x_i)(x-x_{i+1}) \;\Big|_{i=1,\dots,n-1} \;,
	\end{equation}
which identically satisfies the continuity conditions
	\begin{equation}\label{}
S_i(x_i)=y_i \:,\: S_i(x_{i+1})=y_{i+1} \;\Big|_{i=1,\dots,n-1} \;.
	\end{equation}
Substituting (\ref{eq:spl-qs}) into the derivative continuity condition,
	\begin{equation}\label{}
S^\prime_i(x_{i+1})=S^\prime_{i+1}(x_{i+1}) \;\Big|_{i=1,\dots,n-2}\:,
	\end{equation}
gives $n-2$ equations for $n-1$ unknown coefficients $c_i$,
	\begin{eqnarray}\label{eq:spl-qa}
p_i+c_i\Delta x_i= p_{i+1}-c_{i+1}\Delta x_{i+1} \;\Big|_{i=1,\dots,n-2}\:.
	\end{eqnarray}

One coefficient can be chosen arbitrarily, for example $c_1=0$.  The other
coefficients can now be calculated recursively,
	\begin{equation}
c_{i+1}={1\over\Delta x_{i+1}}\left( p_{i+1} - p_{i}-c_i\Delta x_i \right)
\;\Big|_{i=1,\dots,n-2}\;.
	\end{equation}
Alternatively, one can choose $c_{n-1}=0$ and make the backward-recursion
	\begin{equation}
c_{i}={1\over\Delta x_{i}}\left( p_{i+1} - p_{i} - c_{i+1}\Delta x_{i+1}
\right) \;\Big|_{i=n-2,\dots,1}\;.
	\end{equation}

In practice, unless you know what your $c_1$ (or $c_{n-1}$) is, it is
better to run both recursions and then average the resulting $c$'s.
This amounts to first running the forward-recursion from $c_1=0$
and then the backward recursion from $\tfrac12 c_{n-1}$.

The optimized form (\ref{eq:spl-qs}) of the quadratic spline can aslo
be written in the ordinary form, suitable for differentiation and
integration
	\begin{equation}\label{}
S_i(x) = y_i + b_i(x-x_i)+c_i(x-x_i)^2 \:,\;\mathrm{where}\;
b_i = p_i-c_i\Delta x_i \:.
	\end{equation}

An implementation of quadratic spline in C is listed in
Table~\ref{interp:lst:qs}

\begin{table}[htb]\label{interp:lst:qs}
\caption{Quadratic spline in C}
\lstinputlisting{c/interpolation/qspline.c}
\end{table}

\subsection{Cubic spline}
Cubic splines are made of third order polynomials,
	\begin{equation}
S_i(x)=y_i+b_i(x-x_i)+c_i(x-x_i)^2+d_i(x-x_i)^3 \;.
	\end{equation}
This form automatically satisfies the first half of continuity
conditions~(\ref{eq:spl-s}): $S_i(x_i)=y_i$.  The second half of
continuity conditions (\ref{eq:spl-s}), $S_i(x_{i+1})=y_{i+1}$, and the
continuity of the first and second derivatives~(\ref{eq:spl-sp}) give a
set of equations,
	\begin{eqnarray}\label{eq:spl-c1}
y_i+b_ih_i+c_ih_i^2+d_ih_i^3 = y_{i+1} \;,\; i=1,\dots,n-1 \nonumber\\
b_i+2c_ih_i+3d_ih_i^2 = b_{i+1} \;,\; i=1,\dots,n-2 \nonumber\\
2c_i+6d_ih_i = 2c_{i+1}\;,\; i=1,\dots,n-2
	\end{eqnarray}
where
	\begin{equation}
h_i=x_{i+1}-x_i\;.
	\end{equation}

The set of equations (\ref{eq:spl-c1}) is a set of $3n-5$ linear equations
for the $3(n-1)$ unknown coefficients $\{a_i,b_i,c_i\mid i=1,\dots,n-1\}$.
Therefore two more equations should be added to the set to find the
coefficients. If the two extra equations are also linear, the total
system is linear and can be easily solved.

The spline is called {\em natural}\index{natural spline} if the extra
conditions are given as vanishing second derivative at the end-points,
	\begin{equation}
S^{\prime\prime}(x_1)= S^{\prime\prime}(x_n)= 0\;,
	\end{equation}
which gives
	\begin{eqnarray}\label{eq:spl-n}
c_1=0 \;, \nonumber \\ c_{n-1}+3d_{n-1}h_{n-1}=0 \;.
	\end{eqnarray}

Solving the first two equations in (\ref{eq:spl-c1}) for $c_i$ and $d_i$
gives\footnote{introducing an auxiliary coefficient $b_n$.}
	\begin{eqnarray}
c_ih_i &=& -2b_i-b_{i+1}+3p_i \nonumber \;, \\
d_ih_i^2 &=& b_i+b_{i+1}-2p_i \;,
	\end{eqnarray}
where $p_i\doteq\frac{\Delta y_i}{h_i}$. The natural conditions
(\ref{eq:spl-n}) and the third equation in
(\ref{eq:spl-c1}) then produce the following tridiagonal system of $n$
linear equations for the $n$ coefficients $b_i$,
	\begin{eqnarray}\label{eq:spl-t}
2b_1+b_2 &=& 3p_1 \;, \nonumber \\
b_i+(2\frac{h_i}{h_{i+1}}+2)b_{i+1}+\frac{h_i}{h_{i+1}}b_{i+2}
 &=& 3(p_i+p_{i+1}\frac{h_i}{h_{i+1}})\;\Big|_{i=1,\dots,n-2}\;,\nonumber \\
b_{n-1}+2b_n &=& 3p_{n-1} \;,
	\end{eqnarray}
or, in the matrix form,
	\begin{equation}
\left(
\begin{array}{cccccc}
D_1 & Q_1 & 0   & 0 & \dots \\
1   & D_2 & Q_2 &0 &\dots \\
0 & 1 & D_3 & Q_3 &\dots \\
\vdots & \vdots & \ddots & \ddots & \ddots \\
\dots & \dots & 0 & 1 & D_n
\end{array}
\right)
\left(\begin{array}{c}
b_1 \\ \vdots \\ \vdots \\ b_n
\end{array} \right)
=
\left(\begin{array}{c}
B_1 \\ \vdots \\ \vdots \\ B_n
\end{array} \right)
	\end{equation}
where the elements $D_i$ at the main diagonal are
	\begin{equation}
D_1 = 2\:,\;
D_{i+1} = 2\frac{h_i}{h_{i+1}}+2\;\Big|_{i=1,\dots,n-2}\:,\;
D_n = 2 \;,
	\end{equation}
the elements $Q_i$ at the above-main diagonal are
	\begin{equation}
Q_1 = 1\:,\;
Q_{i+1} = \frac{h_i}{h_{i+1}} \;\Big|_{i=1,\dots,n-2}\;,
	\end{equation}
and the right-hand side terms $B_i$ are
	\begin{equation}
B_1 = 3p_1\:,\;
B_{i+1} = 3\left(p_i+p_{i+1}\frac{h_i}{h_{i+1}}\right)
\;\Big|_{i=1,\dots,n-2}\:,
B_n = 3p_{n-1} \:.
	\end{equation}

This system can be solved by one run of Gauss elimination and then
a run of back-substitution.  After a run of Gaussian elimination the
system becomes

	\begin{equation}\label{eq:spl-ge}
\left(
\begin{array}{cccccc}
\tilde D_1 & Q_1 & 0   & 0 & \dots \\
0   & \tilde D_2 & Q_2 &0 &\dots \\
0 & 0 & \tilde D_3 & Q_3 &\dots \\
\vdots & \vdots & \ddots & \ddots & \ddots \\
\dots & \dots & 0 & 0 & \tilde D_n
\end{array}
\right)
\left(\begin{array}{c}
b_1 \\ \vdots \\ \vdots \\ b_n
\end{array} \right)
=
\left(\begin{array}{c}
\tilde B_1 \\ \vdots \\ \vdots \\ \tilde B_n
\end{array} \right) \;,
	\end{equation}
where
	\begin{equation}
\tilde D_1=D_1\:,\; \tilde{D}_i=D_i-Q_{i-1}/\tilde D_{i-1}
\;\Big|_{i=2,\dots,n} \;,
	\end{equation}
and
	\begin{equation}
\tilde B_1=B_1 \:,\; \tilde{B}_i=B_i-\tilde B_{i-1}/\tilde D_{i-1}
\;\Big|_{i=2,\dots,n}\;.
	\end{equation}

The triangular system (\ref{eq:spl-ge}) can be solved by a run of
back-substitution,
	\begin{equation}
b_n=\frac{1}{\tilde D_n}\tilde B_n \:,\;
b_i=\frac{1}{\tilde D_i}(\tilde B_i-Q_ib_{i+1})
\;\Big|_{i=n-1,\dots,1}\;.
	\end{equation}

An implementation of quadratic spline in C is listed in
Table~\ref{interp:lst:cs}

\begin{table}[htb]\label{interp:lst:cs}
\caption{Cubic spline in C}
\lstinputlisting{c/interpolation/cspline.c}
\end{table}

\subsection{Akima sub-spline interpolation}
{\em Akima sub-spline}\index{Akima sub-spline} is an interpolating
function in the form of a piecewise cubic polynomial~\cite{akima},
similar to a cubic spline,
	\begin{equation}\label{interp:eq:ss}
\mathcal{A}(x)\Big|_{x\in [x_i,x_{i+1}]}
=a_i+b_i(x-x_i)+c_i(x-x_i)^2+d_i(x-x_i)^3 \doteq A_i(x) \:.
	\end{equation}
However, unlike the cubic spline, Akima sub-spline dispenses with the
demand of maximal differentiability of the spline---in this case, the
continuity of the second derivative---hence the name {\em sub-spline}.
	\begin{figure}[htb]
\centerline{\includegraphics{akima.eps}} \caption{A cubic spline (solid
line) showing the typical wiggles, compared to the Akima sub-spline
(dashed line) where the wiggles are essentially removed.}
	\label{interp:fig:akima}
	\end{figure}
Instead of achieving maximal differentiability Akima sub-splines try
to reduce the wiggling the ordinary splines are typically prone to
(see figure~\ref{interp:fig:akima}).

First let us note that the coefficients $\{a_i,b_i,c_i,d_i\}$ in
eq.~(\ref{interp:eq:ss}) are determined by the values of the derivatives
$\mathcal{A}'_i\doteq \mathcal{A}'(x_i)$ of the sub-spline through the
continuity conditions for the sub-spline and its first derivative,
	\begin{equation}\label{interp:eq:ac}
A_i(x_i)=y_i ,\: A'_i(x_i)=\mathcal{A}'_i ,\: A_i(x_{i+1})=y_{i+1} ,\:
A'_i(x_{i+1})=\mathcal{A}'_{i+1} .
	\end{equation}
Inserting~(\ref{interp:eq:ss}) into~(\ref{interp:eq:ac}) and solving for
the coefficients gives
	\begin{equation}\label{}
a_i=y_i,\; b_i=\mathcal{A}'_i ,\:
c_i=\frac{3p_i -2\mathcal{A}'_i-\mathcal{A}'_{i+1}}{\Delta x_i} ,\:
d_i=\frac{\mathcal{A}'_i+\mathcal{A}'_{i+1}-2p_i}{(\Delta x_i)^2} ,
	\end{equation}
where $p_i\doteq\Delta y_i/\Delta x_i$, $\Delta y_i\doteq y_{i+1}-y_i$,
$\Delta x_i\doteq x_{i+1}-x_i$.

In the ordinary qubic spline the derivatives $\mathcal{A}'_i$ are
determined by the continuity condition of the second derivative of
the spline.  Sub-splines do without this continuity condition and can
instead use the derivatives as free parameters to be chosen to satisfy
some other condition.

Akima suggested to minimize the wiggling by choosing the
derivatives as linear combinations of the nearest slopes,
	\begin{eqnarray}\label{}
\mathcal{A}'_i = \frac{w_{i+1}p_{i-1}+w_{i-1}p_{i}}{w_{i+1}+w_{i-1}} \:,
&\mathrm{if}& w_{i+1}+w_{i-1}\ne 0 \:,\\
\mathcal{A}'_i = \frac{p_{i-1}+p_{i}}{2} \:,
&\mathrm{if}& w_{i+1}+w_{i-1}= 0 \:,
	\end{eqnarray}
where the weights $w_i$ are given as
	\begin{equation}\label{}
w_i = |p_{i}-p_{i-1}| \:.
	\end{equation}
The idea is that if three points lie close to a line, the sub-spline
in this vicinity has to be close to this line.  In other words, if
$|p_{i}-p_{i-1}|$ is small, the nearby derivatives must be close to
$p_i$.

The first two and the last two points need a special prescription, for
example (naively) one can simply use
	\begin{equation}\label{}
\mathcal{A}'_1 = p_1 ,\:
\mathcal{A}'_1 = \frac12 p_1 + \frac12 p2,\:
\mathcal{A}'_n = p_{n-1} ,\:
\mathcal{A}'_{n-1} = \frac12 p_{n-1} + \frac12 p_{n-2}.
	\end{equation}

Table~(\ref{interp:t:a}) shows a javascript implementation of this
algorithm.
	\begin{table}[htb]
	\caption{Javascript implementation of Akima sub-spline.}
	\label{interp:t:a}
\lstinputlisting{javascript/interpolation/akima.js}
	\end{table}

\section{Other forms of interpolation}
Other forms of interpolation can be constructed by choosing different
classes of interpolating functions, for example, rational function
interpolation, trigonometric interpolation, wavelet interpolation etc.

Sometimes not only the values of the function are tabulated but also the
values of its derivative. This extra information can be taken advantage
of when constructing the interpolation function.

\section{Multivariate interpolation}
Interpolation of a function in more than one varibale is called {\em
multivariate interpolation}.  The function of interest is represented
as a set of discrete points in a multidimentional space.  The points
may or may not lie on a regular grid.

\subsection{Nearest-neighbor interpolation}
Nearest-neighbor interpolation approximates the value of the function
at a non-tabulated point by the value at the nearest tabulated point,
yielding a piecewise-constant interpolating function.  It can be used
for both regular and irregular grids.

\subsection{Piecewise-linear interpolation}

Piecewise-linear interpolation is used to interpolate functions of
two variables tabulated on irregular grids.  The tabulated 2D region is
triangulated -- subdivided into a set of non-intersecting triangles whose
union is the original region.  Inside each triangle the interpolating
function $S(x,y)$ is taken in the linear form,
	\begin{equation}
S(x,y)=a + bx + cy \;,
	\end{equation}
where the three constants are determined by the three conditions that
the interpolating function is equal the tabulated values at the three
vertices of the triangle.

\subsection{Bi-linear interpolation}
{\em Bi-linear interpolation} is used to interpolate functions of
two variables tabulated on regular rectilinear 2D grids.

The interpolating function $B(x,y)$ inside each of the grid
rectangles is taken as a product of two linear functions of $x$ and $y$
respectively,
	\begin{equation}
B(x,y)=(\alpha + \beta x)(\gamma + \delta y)
\doteq a + bx + cy +dxy \;,
	\end{equation}
where the four constants $a, b, c, d$ are obtained from the four
conditions that the interpolating function is equal the tabulated values
at the four nearest tabulated points (which are the vertices of the
given grid rectangle).
