WPC4
2BPZCourier 10cpi3|xx6X@8;X@HP LaserJet IIIHPLASIII.PRSx
@,\,~qX@22V@ZCourier 10cpi?xxx,wx6X@8;X@HP LaserJet IIIHPLASIII.PRSx
@,\,~qX@3|x22F`2@
dMore Accurate Linear Interpolation For Table Lookup Algorithms
$John Forkosh
L285 Stegman Parkway #309
Jersey City, NJ 07305
#(201) 860-9677
'Resume
&--------
88 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).
h#I. Introduction
"-----------------
88 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.
88 A lookup table for a function f(x) is a sequence of numbers yi=f(xi),
i=1...n. Using linear interpolation, f(x) is approximated by straight line
segments connecting the yi's. Consider the very simple function f(x)=x.
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.: 0*0*0*
88 To avoid such systematic errors, we suggest line segments connecting a
sequence of "pseudo-values" yi* 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 yi*'s. The result is that the
table does not contain exact values for the function f(x) at the points xi.
However, the mean square error due to sampling random points in the function's
domain is minimized. For instance, for f(x)=x between -10x10 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).
88 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 yi*'s as n simultaneous
equations in terms of the yi's. Complete code is provided so that the reader
can apply the algorithm to any function f(x) whatsoever.
II. Functions of a Single Variable
|------------------------------------
88 Interpolation techniques represent a function y=f(x) by a table of values
yi=f(xi), i=1...n. Between yi and yi+1, standard linear interpolation
approximates f(x) by a line segment, which we'll denote by ui(x)=aix+bi,
xixxi+1, with ai and bi chosen so that ui(xi)=yi and ui(xi+1)=yi+1. This
provides exact values of f(x) at the xi's, but does not constrain the mean
square errors
c!#$s dddddddwxAK`
Eq.(1)Oepsilon_i = 1 over {x_{i+1}-x_i} INT from{x_i}to{x_{i+1}}
{(u_i(x)-f(x))^2dx}.
x6X@8;X@x6X@8;X@x6X@8;X@!#ixoiexoiO;xddiGxdd%i#uHi#xx#fh #x #dx#Modd#hLL1o1dd1S#(#(#)#( #)X
#)
e2#.c$%%X%%!!%$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.0."0*0*0* %#!0
88 Figure 1b illustrates the procedure to be developed. In this case the line
segments ui(x) connect points ui(xi)=yi* and ui(xi+1)=yi+1*. The yi*'s are
determined by a least squares technique that minimizes the sum of the error
terms given in Eq.(1), so that
GA#$ dddddddwxAK`ZEq.(2)oepsilon = smallsum from{i=1}to{n-1} epsilon_i; horz 150
{partial epsilon}over{partial y_i*}=0, horz 25 i=1...n.x6X@8;X@x6X@8;X@x6X@8;X@!L<JJR<,z_,_<N()]nRiKi_yW+ii
n1R1;0,,1> . ..
..G$%%%%!A%$This gives n equations for the n unknown yi*'s.
88 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 yi=f(xx) denotes the
middle of each x interval.
88 To keep the analysis tractable we assume that xi+1-xi=x is constant. Then
the line segments can be written
AAAA>a#$s>ddddd}ddwxtu_i(x) = 1over{DELTA x}
left[ (x_{i+1}-x)y_i* horz 25 +
horz 25 (x-x_i)y_{i+1}* right], horz 50
x_ixx_{i+1},x6X@8;X@x6X@8;X@x6X@8;X@!uifx8xxMixy"i>
x.xisyidxixxi()U1](1-) ()
1,1,V=r
;
1!a7ksSkz 8>$%%%% !a%$and the sum of the error terms given by Eqs.(1) and (2) is
K#$sXdddddddwxJ`P Eq.(3)epsilon = sum from{i=1}to{n-1} 1 over{x}
int from x_i to x_{i+1}
left lbrace 1over{x}
left [ (x_{i+1}-x)y_i* + (x-x_i)y_{i+1}*
right ] - f(x)
right rbrace ^2 dxx6X@8;X@x6X@8;X@x6X@8;X@!#\#JJodd<#q
#
#Q#v###I?LL}{} L6sznoix;xddiGxdd%ix#xLi#x #y!
i#x#xF
i#yi#f#x##dx1o11ddR1T1\#(1, #)a#(
#)&1#(#)2K$%%%%!%$A
88 At the interior points 1*)<,z_,_c,/_,_Jc,_,-_L<cc_yW+ix/i_y+i/i`_y+i
i~n/1, 1
<<,:$%%~%%!%$while at the boundaries i=1 and i=n only one term contributes, so that
#$'ddddd_
ddwxJ`Eq.(4b)M{partial }over{partial y_1*}={partial _1}over{partial y_1*},
horz 50
{partial }over{partial y_n*}=
{partial _{n-1}}over{partial y_n*}.x6X@8;X@x6X@8;X@x6X@8;X@ l )>%<,_,_F@c,_,!_?<,_,_Bc,S /H_,u _X<c<c_yY_y-_y+n /n_y% +na+1/1+1, /1
.M$%%`"%%!%$pH&0*0*0*Q %A>%aX%.."%R%'%4*p
88 These derivatives can be evaluated directly as follows:
#$sdddddtddwx{partial _{i-1}}over{partial y_i*}=2over{x^2} int from
x_{i-1}to x_i left lbrace 1over{x} left[ (x-x_{i-1})y_i* +
(x_i-x)y_{i-1}* right] - f(x) right rbrace (x-x_{i-1})dx,x6X@8;X@x6X@8;X@x6X@8;X@>NLLL}w} L2szk,|sq,9#dd%H# m##
#r###,siyNoix;xddiGxdd%ix#x#x8 i
#yi#xR
i#x
#yi#f#xG#x7#xi#dxs1n22ddN%1P1X#( 1(
#)]#(#)"1#(#)#(O1#)#,$%%%%!%$#$s
dddddhddwx]{partial _i}over{partial y_i*}=2over{x^2} int from x_i to
x_{i+1} left lbrace 1over{x} left[ (x_{i+1}-x)y_i* +
(x-x_i)y_{i+1}* right] - f(x) right rbrace (x_{i+1}-x)dx.x6X@8;X@x6X@8;X@x6X@8;X@)NLLL~}k} L&sz,k,-#dd, #a##A
#f###vsiyHoiyxz;xddiGxdd%ix#x<i #x
#yi#x
#x6i#yvi#f{#x;#xi#x#dxb2
2ddB1D1L#(1
#)Q#(#)1#(#)#(S1#)#.]$%%%%!%$Substitute xi-1=xi-x and xi+1=xi+x to obtain:
-!#$sndddddddwx{partial _{i-1}}over{partial y_i*}=2over{x^3} int from
{x_i-x}to x_i left lbrace (y_i*-y_{i-1}*)(x-x_i+x) +
y_{i-1}*x - f(x)x right rbrace (x-x_i+x)dx,x6X@8;X@x6X@8;X@x6X@8;X@>NLLLk,|sq,9#G##R ##
##X##!#f#gG
#p###,siyNoix ;xdd^iGxdd%iGx#yJi#y iZ#xJ#xi(#x#yi #x#f#xr#x#x#xiw#xg#dxs1n23U#( 1j
#)
#(#)1q#(a#)1#(#)W#,-$%%
%%!!%$ A#$sVdddddddwx{partial _i}over{partial y_i*}=2over{x^3} int from x_i to
{x_i+x} left lbrace (y_{i+1}*-y_i*)(x-x_i-x) + y_{i+1}*x -
f(x)x right rbrace (x-x_i-x)dx.x6X@8;X@x6X@8;X@x6X@8;X@)NLLL,k,-#;)## ##
##L#u##Z#[;
#d###vsiyHoiyxz;xddi;xGxddR%i#y9i #y iN#x>#xi#x#yi#x#f#xf#x#x#x
ik#x[#dxb2
3I#(1^
#)
#(#)1e#(U#)%#(#)K#. $%%%%!A%$Now separate the f(x) term from each integral and substitute x*=x-xi+x on the
top line, x*=x-xi-x on the bottom line:
a#$s2dddddhddwx0{partial _{i-1}}over{partial y_i*}=2over{x^3} int from 0 to
{+x} left lbrace (y_i*-y_{i-1}*)x* + y_{i-1}*x right rbrace
x*dx* horz 50 - horz 50 2over{x^2} int from{x_i-x} to x_i
f(x)(x-x_i+x)dx,x6X@8;X@x6X@8;X@x6X@8;X@>]LLLLLz,s,H###p #
#P#0
#x###G##
#NG#;siy]oixKxK#yi#yi`
#x#y@iA#x#x#dxlx;xddEimGxdd%iGx#f#x\#xL#xi*#x#dxs1}2(3G0#( 1 #)1U22|#(l#)#(#)
#,0߭$%%%%!a%$#$sddddd\ddwx{partial _i}over{partial y_i*}=2over{x^3} int from 0 to
{-x} left lbrace (y_{i+1}*-y_i*)x* + y_{i+1}*x right rbrace
x*dx* horz 50 - horz 50 2over{x^2} int from x_i to {x_i+x}
f(x)(x_i-x+x)dx.x6X@8;X@x6X@8;X@x6X@8;X@)]LLLLL,z,<###d #
#D#$
#l###;#
#$
#B;#siyWoix?x?#yi#y iT
#x#y4i5#x#x#dx`xa;xddi;xGxdd9%i#f#xP#xi#x#x#dxq23G0#(W1 #)1I22p#(`#)#(#)#.߉$%%j%%!%$The left-hand integrals can be evaluated directly:
#$s!dddddddwxJ`bEq.(5a){partial _{i-1}}over{partial y_i*} horz 50 = horz 50 2 over 3
y_i* + 1 over 3 y_{i-1}* horz 50 - horz 50 2over{x^2} int
from{x_i-x} to x_i f(x)(x-x_i+x)dx,x6X@8;X@x6X@8;X@x6X@8;X@>lL
L
BL2
LvL,s,#4##\#8 #zGP##O
G
#Jsiyloig#yi#yli
xl;xddiGxdd>%i0
Gx
#fp#x#x#xEi#x#dxs123_1_31
2|2
#(#)`#(#)#,$%%L%%!%$#$s%dddddddwxJ`gEq.(5b)l{partial _i}over{partial y_i*} horz 50 = horz 50 2 over 3
y_i* + 1 over 3 y_{i+1}* horz 50 - horz 50 2over{x^2} int
from x_i to {x_i+x} f(x)(x_i-x+x)dx.x6X@8;X@x6X@8;X@x6X@8;X@)lL
L
6L&
LjL,,#(##P#, #n;##3C
;#siyfoi[#yi#y`i
x;xdd2i$
;x`Gxdd%it
#fd#x#xIi#x#x#dx23S1S31
2p2
#(#)T#(#)z#.l$%%4!%%!%$%0*0*0*%z
%bn%D!V%,A2%a% !%%%%)
88 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 yi=f(xx) so that
k#$sddddd
ddwxdint from{x_i-x} to x_i f(x)dx horz 50
horz 50
{x} over 6 (y_{i-1} + 4y_{i-} + y_i),x6X@8;X@x6X@8;X@x6X@8;X@LL;xddi!Gxddv%ihGx#f#x#dxx#y!i
#yy
i#yiG#
q #
#G0#( #)O61#(1 #4#)I
#,k$%%%%!%$from which we can evaluate
#$sz
dddddWddwxJ`7Eq.(6a)s2 over {x^2} int from {x_i-x} to x_i f(x)(x-x_i+x)dx
horz 50
horz 50
1 over 3 (y_i + 2y_{i-}),x6X@8;X@x6X@8;X@x6X@8;X@NLL
L22#(#)|#(: #)13#(o#2r#)#,kG)#x;xddiGxddZ%iLGx#f#x#x#xai#x #dx*
#y
i#y_iGl###
#
$%%%%!%$
!#$sbdddddWddwxJ`7Eq.(6b)s2 over {x^2} int from x_i to {x_i+x} f(x)(x_i-x+x)dx
horz 50
horz 50
1 over 3 (y_i + 2y_{i+}).x6X@8;X@x6X@8;X@x6X@8;X@NLL
L22#(#)|#(: #)13#(o#2r#)#.k;)#x;xddZiL;xGxdd%i#f#x#xqi9#x#x #dx*
#y
i#y_i;###
#
$%%%%!!%$
88 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
A#$2dddddkddwxJ``
Eq.(7a){partial _{i-1}}over{partial y_i*} horz 50
horz 50
2 over 3 y_i* + 1 over 3 y_{i-1}* horz 50 - horz 50
1 over 3 ( y_i + 2 y_{i-}),x6X@8;X@x6X@8;X@x6X@8;X@>l
B
2
c,/_,_
4\8 )
cJ/i_yl+igyiyli\yi
y
i/1<2_3_<1__31O
<1O
_3
(2),1$%%%%!A%$a#$ddddd_ddwxJ`d
Eq.(7b)j{partial _i}over{partial y_i*} horz 50
horz 50
2 over 3 y_i* + 1 over 3 y_{i+1}* horz 50 - horz 50
1 over 3 ( y_i + 2 y_{i+}).x6X@8;X@x6X@8;X@x6X@8;X@)l
6
&
c,_,_
(P,
3c/i_yf+i[yiy`iPyi
y
i<2_3S<1S_31C
<1C
_3
(2).%j$%%j%%!a%$
88 Eqs.(2), (4) and (7) now let us write down n equations for the n unknown
yi*'s in terms of the known yi's:
b#$q#dddddddwxhorz 400 2y_1* + y_2* horz 50 = horz 50 y_1 + 2y_{1+}, horz
500 for horz 25 i=1, #
y_{i-1}* + 4y_i* + y_{i+1}* horz 50 = horz 50 2y_{i-} + 2y_i
+ 2y_{i+}, horz 200 1** a#%'These results are most easily written in matrix form:
#$gdddddddwxJ`Eq.(8)h&left( matrix{2&1&.&.&.&.&. # 1&4&1&.&.&0&. # .&1&4&1&.&.&. #
.&.&1&4&1&.&. # .&.&.&1&4&1&. # .&0&.&.&1&4&1 # .&.&.&.&.&1&2}
right) horz 50 left( matrix{y_1*#.#.#y_i*#.#.#y_n*} right)
horz 50 = horz 50 left( matrix{y_1+2y_{1+} #.#.#
2y_{i-}+2y_i+2y_{i+} #.#.# y_n+2y_{n-}} right)x6X@8;X@x6X@8;X@x6X@8;X@*hjj=jjjPjjjcjjjvj'jjj:jMi@*o@q@q@=q@q@q@Pq@q@q@cq@q@q@vq@'q@q@q@:q@Mph>jjjQjjjdjjjwj(jjj;jjjNjjixox>qxqxqxQqxqxqxdqxqxqxwqx(qxqxqx;qxqxqxNqxqxp
h
>j
j
j
Qj
j
j
dj
j
j
wj
(j
j
j
;j
j
j
Nj
j
io>qqqQqqqdqqqwq(qqq;qqqNqqpR2;1$.
....R1;4$1
..0.R.;1$4
1...R.;.$1
41..R.;.$.
141.Rz.;z0$z.
z.z1z4z1Rr.;r.$r.
r.r.r1r21..o.g.;1
2
1(
.(
.G
222(
o.(
g.
_25y3yji3_y+ny{
y
y7ji:
y
jiyoji_y=+n}
_y
+n_+ CjJj_E+jj+h$%%%%` !%$
88 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.
III. Discussion of Code
-------------------------
88 The implementation of Eq.(8) is demonstrated in Listings1 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 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.
88 By supplying the address of your own function to lint1d()'s f-argument, an0(0*0*0*%t0interpolation 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.
88 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 !
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.
*