% This file was created by the WP2LaTeX program version: 2.95 
\documentclass[11pt]{article}
\usepackage{wp2latex}


\begin{document}
\baselineskip=2.50ex \nwln
\begin{center}
More Accurate Linear Interpolation For Table Lookup Algorithms
\end{center}
\nwln
\begin{center}
John Forkosh\\
285 Stegman Parkway \#309\\
Jersey City, NJ 07305\\
(201) 860$-$9677
\end{center}

\nwln
\begin{center}
Resume\\
$-$$-$$-$$-$$-$$-$$-$$-$
\end{center}
\hspace*{1.28cm}John Forkosh has an ms in physics, and for the past ten years has been an
independent consultant working primarily under VAX/VMS in "C".  He developed the
algorithm discussed in this article while at the Goddard Institute for Space
Studies in New York City, working on the radiative transfer code for a general
circulation model of the earth's atmosphere (used to study climatic change,
particularly due to increased carbon dioxide resulting from burning fossil
fuels).

\nwln
\begin{center}
I. Introduction\\
$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$
\end{center}
\hspace*{1.28cm}Scientific applications often require the very frequent evaluation of
complicated functions.  To reduce computational overhead, such a function can
be defined by a table of its values at various fixed points.  Some form of
interpolation is then required to estimate values of the function at
intermediate points.  Linear interpolation is the quickest procedure; however,
it is usually the least accurate.

\hspace*{1.28cm}A lookup table for a function f(x) is a sequence of numbers y$_{{\rm i}}$=f(x$_{{\rm i}}$),
i=1...n.  Using linear interpolation, f(x) is approximated by straight line
segments connecting the y$_{{\rm i}}$'s.  Consider the very simple function f(x)=x$^{{\rm ^{2}}{}}$. 
Because this function is concave (its second derviative f"(x)=2 is always
positive) the straight line segments will always lie above the curve, i.e.,
linearly interpolated values will always be too large.
\newpage

 \bigskip
\hspace*{1.28cm}To avoid such systematic errors, we suggest line segments connecting a
sequence of "pseudo$-$values" y$_{{\rm i}}$* that are determined by minimizing the sum of the
integrated square errors attributable to each segment.  A least$-$squares
technique is used to determine the appropriate y$_{{\rm i}}$*'s.  The result is that the
table does not contain exact values for the function f(x) at the points x$_{{\rm i}}$. 
However, the mean square error due to sampling random points in the function's
domain is minimized.  For instance, for f(x)=x$^{{\rm ^{2}}{}}$ between $-$10$\leq$x$\leq$10 with table
values at each integer, the mean square error is reduced by a factor of 6 (the
rms error by a factor of 2.45).

\hspace*{1.28cm}Our discussion is pitched at about the level of Numerical Recipes in C,
Press et.al., Cambridge University Press 1988.  It illustrates not only
interpolation (Chapter 3), but also least squares (Chapter 14) and the solution
of linear algebraic equations (Chapter 2).  The latter arises because, for a
table with n values, our analysis expresses the y$_{{\rm i}}$*'s as n simultaneous
equations in terms of the y$_{{\rm i}}$'s.  Complete code is provided so that the reader
can apply the algorithm to any function f(x) whatsoever.

\nwln
\begin{center}
II. Functions of a Single Variable\\
$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$
\end{center}
\hspace*{1.28cm}Interpolation techniques represent a function y=f(x) by a table of values
y$_{{\rm i}}$=f(x$_{{\rm i}}$), i=1...n.  Between y$_{{\rm i}}$ and y$_{{\rm i+1}}$, standard linear interpolation
approximates f(x) by a line segment, which we'll denote by u$_{{\rm i}}$(x)=a$_{{\rm i}}$x+b$_{{\rm i}}$,
x$_{{\rm i}}$$\leq$x$\leq$x$_{{\rm i+1}}$, with a$_{{\rm i}}$ and b$_{{\rm i}}$ chosen so that u$_{{\rm i}}$(x$_{{\rm i}}$)=y$_{{\rm i}}$ and u$_{{\rm i}}$(x$_{{\rm i+1}}$)=y$_{{\rm i+1}}$.  This
provides exact values of f(x) at the x$_{{\rm i}}$'s, but does not constrain the mean
square errors

\begin{equation}
\epsilon_{i} ={\frac{1}{x_{i+1} -x_{i}}} \int_{x_{i}}^{x_{i+1} } {(u_{i} (x)-f(x))^{2} dx}.
\end{equation}
In addition, it gives rise to systematic errors when the sign of f"(x) doesn't
change between points.  If f"$>$0 (concave function) then u$-$f$>$0 for all x in the
interval (and similarly for a convex function).  This situation is illustrated
in figure 1a.
\newpage

 \bigskip
\hspace*{1.28cm}Figure 1b illustrates the procedure to be developed.  In this case the line
segments u$_{{\rm i}}$(x) connect points u$_{{\rm i}}$(x$_{{\rm i}}$)=y$_{{\rm i}}$* and u$_{{\rm i}}$(x$_{{\rm i+1}}$)=y$_{{\rm i+1}}$*.  The y$_{{\rm i}}$*'s are
determined by a least squares technique that minimizes the sum of the error
terms given in Eq.(1), so that

\begin{equation}
\epsilon =\sum _{i=1}^{n-1} \epsilon_{i} ;\hspace{0.32cm}{\frac{\partial \epsilon }{\partial y_{i} *}} =0,\hspace{0.05cm} i=1...n.
\end{equation}
This gives n equations for the n unknown y$_{{\rm i}}$*'s.

\hspace*{1.28cm}The work between here and Eq.(8) contains the algebra necessary to
formulate these n equations so that we can program their solution.  The reader
who doesn't care to follow this derivation can skip to Eq.(8) which contains the
required result.  If you do skip there, note that y$_{{\rm i}{\rm \pm}{}{\rm \frac{1}{2}}{}}$=f(x$\pm$$\frac{1}{2}$$\Delta$x) denotes the
middle of each $\Delta$x interval.

\hspace*{1.28cm}To keep the analysis tractable we assume that x$_{{\rm i+1}}$$-$x$_{{\rm i}}$=$\Delta$x is constant.  Then
the line segments can be written

\begin{displaymath}
u_{i} (x) ={\frac{1}{\Delta x}} \left [ (x_{i+1} -x)y_{i} *\hspace{0.05cm} +\hspace{0.05cm} (x-x_{i} )y_{i+1} * \right ],\hspace{0.11cm} x_{i} \leq x\leq x_{i+1} ,
\end{displaymath}
and the sum of the error terms given by Eqs.(1) and (2) is

\begin{equation}
\epsilon = \sum_{i=1}^{n-1} {\frac{1}{\Delta x}} \int_{x_{i}}^{x_{i+1} } \left\{{\frac{1}{\Delta x}} \left [ (x_{i+1} -x)y_{i} * + (x-x_{i} )y_{i+1} * \right ] - f(x) \right\}^{2} dx
\end{equation}

 \bigskip
\hspace*{1.28cm}At the interior points 1$<$i$<$n two terms from the sum in Eq.(3) contain
subscript i, so that

\begin{equation}
{\frac{\partial \epsilon }{\partial y_{i} *}} \hspace{0.11cm} =\hspace{0.11cm}{\frac{\partial \epsilon_{i-1} }{\partial y_{i} *}} \hspace{0.11cm} +\hspace{0.11cm}{\frac{\partial \epsilon _{i}} {\partial y_{i} *}} ,\hspace{0.11cm} 1<i<n,
\end{equation}
while at the boundaries i=1 and i=n only one term contributes, so that

\begin{equation}
{\frac{\partial \epsilon }{\partial y_{1} *}} ={\frac{\partial \epsilon _{1} }{\partial y_{1} *}} ,\hspace{0.11cm}{\frac{\partial \epsilon }{\partial y_{n} *}} ={\frac{\partial \epsilon _{n-1} }{\partial y_{n} *}} .
\end{equation}

\newpage
\bigskip
\hspace*{1.28cm}These derivatives can be evaluated directly as follows:

\begin{displaymath}
{\frac{\partial \epsilon _{i-1} }{\partial y_{i} *}} ={\frac{2}{\Delta x^{2} }} \int_{x_{i-1} }^{x_{i}} \left\{{\frac{1}{\Delta x}} \left [ (x-x_{i-1} )y_{i} * + (x_{i} -x)y_{i-1} * \right ] - f(x) \right\} (x-x_{i-1} )dx,
\end{displaymath}
\begin{displaymath}
{\frac{\partial \epsilon _{i}} {\partial y_{i} *}} ={\frac{2}{\Delta x^{2} }} \int_{x_{i}}^{x_{i+1} } \left\{{\frac{1}{\Delta x}} \left [ (x_{i+1} -x)y_{i} * + (x-x_{i} )y_{i+1} * \right ] - f(x) \right\} (x_{i+1} -x)dx.
\end{displaymath}
Substitute x$_{{\rm i}{\rm -}{1}}$=x$_{{\rm i}}$$-$$\Delta$x and x$_{{\rm i+1}}$=x$_{{\rm i}}$+$\Delta$x to obtain:

\begin{displaymath}
{\frac{\partial \epsilon _{i-1} }{\partial y_{i} *}} ={\frac{2}{\Delta x^{3} }} \int_{x_{i} -\Delta x}^{x_{i}} \left\{ (y_{i} *-y_{i-1} *)(x-x_{i} +\Delta x) + y_{i-1} *\Delta x - f(x)\Delta x \right\} (x-x_{i} +\Delta x)dx,
\end{displaymath}
\begin{displaymath}
{\frac{\partial \epsilon _{i}} {\partial y_{i} *}} ={\frac{2}{\Delta x^{3} }} \int_{x_{i}}^{x_{i} +\Delta x} \left\{ (y_{i+1} *-y_{i} *)(x-x_{i} -\Delta x) + y_{i+1} *\Delta x - f(x)\Delta x \right\} (x-x_{i} -\Delta x)dx.
\end{displaymath}
Now separate the f(x) term from each integral and substitute x*=x$-$x$_{{\rm i}}$+$\Delta$x on the
top line, x*=x$-$x$_{{\rm i}}$$-$$\Delta$x on the bottom line:

\begin{displaymath}
{\frac{\partial \epsilon _{i-1} }{\partial y_{i} *}} ={\frac{2}{\Delta x^{3} }} \int_{0}^{+\Delta x} \left\{ (y_{i} *-y_{i-1} *)x* + y_{i-1} *\Delta x \right\} x*dx*\hspace{0.11cm} -\hspace{0.11cm}{\frac{2}{\Delta x^{2} }} \int_{x_{i} -\Delta x}^{x_{i}} f(x)(x-x_{i} +\Delta x)dx,
\end{displaymath}
\begin{displaymath}
{\frac{\partial \epsilon _{i}} {\partial y_{i} *}} ={\frac{2}{\Delta x^{3} }} \int_{0}^{-\Delta x} \left\{ (y_{i+1} *-y_{i} *)x* + y_{i+1} *\Delta x \right\} x*dx*\hspace{0.11cm} -\hspace{0.11cm}{\frac{2}{\Delta x^{2} }} \int_{x_{i}}^{x_{i} +\Delta x} f(x)(x_{i} -x+\Delta x)dx.
\end{displaymath}
The left$-$hand integrals can be evaluated directly:

\begin{equation}
{\frac{\partial \epsilon _{i-1} }{\partial y_{i} *}} \hspace{0.11cm} =\hspace{0.11cm}{\frac{2}{3}} y_{i} * +{\frac{1}{3}} y_{i-1} *\hspace{0.11cm} -\hspace{0.11cm}{\frac{2}{\Delta x^{2} }} \int_{x_{i} -\Delta x}^{x_{i}} f(x)(x-x_{i} +\Delta x)dx,
\end{equation}
\begin{equation}
{\frac{\partial \epsilon _{i}} {\partial y_{i} *}} \hspace{0.11cm} =\hspace{0.11cm}{\frac{2}{3}} y_{i} * +{\frac{1}{3}} y_{i+1} *\hspace{0.11cm} -\hspace{0.11cm}{\frac{2}{\Delta x^{2} }} \int_{x_{i}}^{x_{i} +\Delta x} f(x)(x_{i} -x+\Delta x)dx.
\end{equation}

\newpage
\bigskip
\hspace*{1.28cm}So far the analysis is exact.  To finish the problem we approximate the
right$-$hand integrals using Simpson's rule applied at the midpoint of each
interval.  Let y$_{{\rm i}{\rm \pm}{}{\rm \frac{1}{2}}{}}$=f(x$\pm$$\frac{1}{2}$$\Delta$x) so that

\begin{displaymath}
\int_{x_{i} -\Delta x}^{x_{i}} f(x)dx\hspace{0.11cm} \approx\hspace{0.11cm}{\frac{\Delta x}{6}} (y_{i-1} + 4y_{i-\frac{1}{2} } + y_{i} ),
\end{displaymath}
from which we can evaluate

\begin{equation}
{\frac{2}{\Delta x^{2} }} \int_{x_{i} -\Delta x}^{x_{i}} f(x)(x-x_{i} +\Delta x)dx\hspace{0.11cm} \approx\hspace{0.11cm}{\frac{1}{3}} (y_{i} + 2y_{i-\frac{1}{2} } ),
\end{equation}
\begin{equation}
{\frac{2}{\Delta x^{2} }} \int_{x_{i}}^{x_{i} +\Delta x} f(x)(x_{i} -x+\Delta x)dx\hspace{0.11cm} \approx\hspace{0.11cm}{\frac{1}{3}} (y_{i} + 2y_{i+\frac{1}{2} } ).
\end{equation}

 \bigskip
\hspace*{1.28cm}Any desired degree of accuracy can be obtained by using better
approximations in Eqs.(6).  The results are then substituted back into Eqs.(5),
and in this case

\begin{equation}
{\frac{\partial \epsilon _{i-1} }{\partial y_{i} *}} \hspace{0.11cm} \approx\hspace{0.11cm}{\frac{2}{3}} y_{i} * +{\frac{1}{3}} y_{i-1} *\hspace{0.11cm} -\hspace{0.11cm}{\frac{1}{3}} ( y_{i} + 2 y_{i-\frac{1}{2} } ),
\end{equation}
\begin{equation}
{\frac{\partial \epsilon _{i}} {\partial y_{i} *}} \hspace{0.11cm} \approx\hspace{0.11cm}{\frac{2}{3}} y_{i} * +{\frac{1}{3}} y_{i+1} *\hspace{0.11cm} -\hspace{0.11cm}{\frac{1}{3}} ( y_{i} + 2 y_{i+\frac{1}{2} } ).
\end{equation}

 \bigskip
\hspace*{1.28cm}Eqs.(2), (4) and (7) now let us write down n equations for the n unknown
y$_{{\rm i}}$*'s in terms of the known y$_{{\rm i}}$'s:

\begin{eqnarray*}
\hspace{0.85cm} 2y_{1} * + y_{2} *\hspace{0.11cm} =\hspace{0.11cm} y_{1} + 2y_{1+\frac{1}{2} } ,\hspace{1.06cm} for\hspace{0.05cm} i=1,  \nonumber  \\
y_{i-1} * + 4y_{i} * + y_{i+1} *\hspace{0.11cm} =\hspace{0.11cm} 2y_{i-\frac{1}{2} } + 2y_{i} + 2y_{i+\frac{1}{2} } ,\hspace{0.43cm} 1<i<n,  \nonumber  \\
\hspace{0.64cm} 2y_{n} * + y_{n-1} *\hspace{0.11cm} =\hspace{0.11cm} y_{n} + 2y_{n-\frac{1}{2} } ,\hspace{1.49cm} i=n.
\end{eqnarray*}

\newpage
These results are most easily written in matrix form:

\begin{equation}
\left (\begin{array}{ccccccc}2&1&.&.&.&.&.  \\
1&4&1&.&.&0&.  \\
.&1&4&1&.&.&.  \\
.&.&1&4&1&.&.  \\
.&.&.&1&4&1&.  \\
.&0&.&.&1&4&1  \\
.&.&.&.&.&1&2\end{array} \right )\hspace{0.11cm} \left (\begin{array}{c}y_{1} * \\
. \\
. \\
y_{i} * \\
. \\
. \\
y_{n} *\end{array} \right )\hspace{0.11cm} =\hspace{0.11cm} \left (\begin{array}{c}y_{1} +2y_{1+\frac{1}{2} }  \\
. \\
. \\
2y_{i-\frac{1}{2} } +2y_{i} +2y_{i+\frac{1}{2} }  \\
. \\
. \\
y_{n} +2y_{n-\frac{1}{2} } \end{array} \right )
\end{equation}

 \bigskip
\hspace*{1.28cm}Eq.(8) specifies how to construct a table of y*'s, for linearly
interpolating any function f(x), that minimizes (modulo Simpson's rule) the mean
square error due to randomly sampling points in the function's domain.  However,
f(x) must be a function of a single variable.  The application for which this
alogorithm was originally designed requires interpolation of a function z=f(x,y)
of two variables.  This more general situation is conceptually quite similar to
the one$-$dimensional case we discussed; however, the required algebra becomes
noticeably more complicated.  I therefore decided to stop the formal discussion
of the algorithm at this point.  However, if some unanticipated level of
interest arises, I'll be happy to supply CUJ with a short discussion and the
necessary code.

\nwln
\begin{center}
III. Discussion of Code\\
$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$$-$
\end{center}
\hspace*{1.28cm}The implementation of Eq.(8) is demonstrated in Listings~1 and 2.  Listing
1 contains a "C" function lint1d(x0,dx,n,f,y) that returns a table of y*'s,
suitable for one$-$dimensional linear interpolation, as prescribed by Eq.(8).  It
simply malloc's a temporary nxn array in which the coefficient matrix is set up,
then calls the user$-$supplied function f() to initialize array y[] with the
vector of constants, and finally solves the simultaneous equations, returning
the y* solutions to the caller in the same array y[].  Listing 1 also contains
a test driver for lint1d(), and some sample output for f(x)=x$^{{\rm ^{2}}{}}$ as discussed in
the introduction.  Listing 2 contains the simultaneous equation solver used by
lint1d().  It's a Fortran$-$to$-$C port of an old IBM Scientific Subroutine Package
routine, simq(), as noted in the code.  As you can see, Eq.(8) specifies a very
straightforward procedure, and the resulting code for lint1d() is
correspondingly simple.  The code for simq() is a bit dense, but can be treated
as a black box for the purposes of this application.

\hspace*{1.28cm}By supplying the address of your own function to lint1d()'s f$-$argument, an
interpolation table can very easily be constructed to suit any purpose
whatsoever.  Thus, almost any application that employs one$-$dimensional linear
interpolation at fixed intervals should be able to realize noticeable
improvements in accuracy.  All it takes is a simple one$-$line call to lint1d(). 
Afterwards, interpolation on the improved table is performed totally
transparently, so absolutely no application$-$level recoding is required.

\hspace*{1.28cm}Both listings contain embedded drivers (that are compiled when the \#defined
flag TESTDRIVE is true) to facilitate module$-$level testing and to illustrate
usage.  I've also tried to carefully document the code (except for some of the
more elaborate indexing arithmetic required near the bottom of simq()).  My
personal experience with scientific code is that it frequently tends to be
under$-$documented and also under$-$robust.  The latter problem can be particularly
pernicious since the code usually won't blow up, but will simply supply the
wrong result due to an inappropriate numerical method or due to some more
trivial problem.  (For instance, I once worked on a model where the value of $\pi$
had been inadvertantly entered incorrectly.)  Since numerical results are
usually not analytically verifiable, and since feedback mechanisms frequently
prevent small errors from growing uncontrollably, it may take quite some time
(if ever) before a researcher realizes that his results are little more than
massaged noise.  The only effective procedure to guard against such difficulties
is a design that carefully decomposes a complicated calculation into easily
describable functional components that can be individually tested for known
correct results.
\end{document}