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@   d More 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 ddddd ddwxAK` 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#$ ddddd ddwxAK`Z Eq.(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 AA A A>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. x is y idxixxi()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< c c _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 #y i #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 #y i #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 #x i(#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> #x i#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 #y4 i5#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 Lv L,s,#4##\#8 #z GP## O  G #Jsiyloig#yi#yli xl ;xdd i Gxdd> %i0 Gx #fp#x#x#xEi#x#dxs123_1_3 1 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& Lj L,,#(##P#, #n ;##3 C  ;#siyfoi[#yi#y`i x ;xdd2 i$ ;x` Gxdd %it #fd#x#xIi#x#x#dx23S1S31 2p 2 #(#)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 #y iG# q #  #G0#( #)O61#(1 #4 #)I #, k$%%%%!%$from which we can evaluate  #$sz dddddWddwxJ`7 Eq.(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#(#)|#(: #) 1 3 #(o#2r#)#,kG)#x;xddiGxddZ%iLGx#f#x#x#xai#x #dx* #y i#y_iGl## #  # $%%%%!%$ !#$sbdddddWddwxJ`7 Eq.(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#(#)|#(: #) 1 3 #(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 )  c J/i_yl+igyiyli\ y i y i/1<2_3_<1__3 1O <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`iP y i 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 2 22( o.( g. _25y3yji3_y+n y{ y y7 ji: y jiyoji _y= +n} _y +n_+  C jJ j _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*%t0 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. 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.