f3-6

116 Chapter3. InterpolationandExtrapolationSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).If the point x=0is not in (or at least close to) the range of the tabulated xi’s, thenthecoefficientsoftheinterpolatingpolynomialwillingeneralbecomeverylarge.However, the real “information content” of the coefficients is in small differencesfrom the “translation-induced” large values. This is one cause of ill-conditioning, resulting in loss of significance and poorly determined coefficients. You should consider redefiningthe origin of the problem,to put x=0in a sensible place. Another pathologyis that, if too high a degree of interpolationis attempted on a smooth function, the interpolating polynomial will attempt to use its high-degree coefficients,incombinationswithlargeandalmostpreciselycancelingcombinations, to match the tabulated values down to the last possible epsilon of accuracy. This effect is the same as the intrinsic tendencyof the interpolatingpolynomialvalues tooscillate (wildly) between its constrained points, and would be present even if the machine’s floating precision were infinitely good. The above routines polcoeand polcofhave slightly different sensitivities to the pathologies that can occur. Are you still quite certain that using the coefficients is a good idea? CITED REFERENCES AND FURTHER READING: Isaacson, E., and Keller, H.B. 1966, Analysis of Numerical Methods (New York: Wiley), §5.2. 3.6 Interpolation in Two or More Dimensions In multidimensional interpolation, we seek an estimate of y(x1,x2,...,x n) from an n-dimensional grid of tabulated values yand none-dimensional vec- tors giving the tabulated values of each of the independent variables x1,x2,..., xn. We will not here consider the problem of interpolating on a mesh that is not Cartesian, i.e., has tabulated function values at “random” points in n-dimensional space rather than at the vertices of a rectangular array. For clarity, we will consider explicitly only the case of two dimensions, the cases of three or more dimensions being analogous in every way. In two dimensions, we imagine that we are given a matrix of functionalvalues ya(j,k), where jvaries from 1 to m, and kvaries from 1 to n. We are also given an array x1aof length m, and an array x2aof length n. The relation of these input quantities to an underlying function y(x1,x2)is ya(j,k) =y(x1a(j) ,x2a(k) )( 3.6.1 ) We want to estimate, by interpolation, the function yat some untabulated point (x1,x2). An important concept is that of the grid square in which the point (x1,x2) falls, that is, the four tabulated points that surround the desired interior point. For convenience, we will number these points from 1 to 4, counterclockwise startingfrom the lower left (see Figure 3.6.1). More precisely, if x1a(j) ≤x 1≤x1a(j+1) x2a(k) ≤x2≤x2a(k+1)(3.6.2 ) 3.6InterpolationinTwoor MoreDimensions 117Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).y ∂y/∂x1 ∂y/∂x2 ∂2y/∂x1∂x2x2 = x2u x2 = x2lx1 = x1ux1 = x1l1234 pt. 1user suppliesthese valuespt. 4 pt. 2pt. 3 d2 d1 (a) (b)⊗desired pt. (x1,x2)pt. number Figure 3.6.1. (a) Labeling of points used in the two-dimensional interpolation routines bcuintand bcucof. (b) For each of the four points in (a), the user supplies one function value, two first derivatives, and one cross-derivative, a total of 16 numbers. definesjandk, then y1≡ya(j,k) y2≡ya(j+1,k) y3≡ya(j+1,k+1) y4≡ya(j,k+1)(3.6.3 ) The simplest interpolation in two dimensions is bilinear interpolation on the grid square. Its formulas are: t≡(x1−x1a(j) )/(x1a(j+1) −x1a(j) ) u≡(x2−x2a(k) )/(x2a(k+1) −x2a(k) )(3.6.4 ) (so that tand ueach lie between 0 and 1), and y(x1,x2)=( 1 −t)(1−u)y1+t(1−u)y2+tuy 3+( 1−t)uy 4 (3.6.5 ) Bilinear interpolation is frequently “close enough for government work. ”As the interpolating point wanders from grid square to grid square, the interpolated function value changes continuously. However, the gradient of the interpolated function changes discontinuously at the boundaries of each grid square. There are two distinctly different directions that one can take in going beyond bilinear interpolation to higher-order methods: One can use higher order to obtain increasedaccuracyfortheinterpolatedfunction(forsuf ficientlysmoothfunctions!), without necessarily trying to fix up the continuity of the gradient and higher derivatives. Or,one canmake use of higherorderto enforcesmoothnessof someof these derivatives as the interpolating point crosses grid-square boundaries. We will now consider each of these two directions in turn. 118 Chapter3. InterpolationandExtrapolationSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).HigherOrderfor Accuracy Thebasic idea is to breakupthe problemintoa succession ofone-dimensional interpolations. If wewant to do m-1orderinterpolationin the x1direction,and n-1 orderinthe x2direction,we first locatean m×nsub-blockofthetabulatedfunction matrix that contains our desired point (x1,x2). We then do mone-dimensional interpolations in the x2direction, i.e., on the rows of the sub-block, to get function values at the points (x1a(j) ,x2),j=1,..., m. Finally, we do a last interpolation in the x1direction to get the answer. If we use the polynomialinterpolationroutine polintof§3.1,andasub-blockwhichispresumedtobealreadylocated(andcopied into an mbynarray ya), the procedure looks like this: SUBROUTINE polin2(x1a,x2a,ya,m,n,x1,x2,y,dy) INTEGER m,n,NMAX,MMAX REAL dy,x1,x2,y,x1a(m),x2a(n),ya(m,n) PARAMETER (NMAX=20,MMAX=20) Maximum expected values of nandm. C USES polint Given arrays x1a(1:m) andx2a(1:n) of independent variables, and an mbynarray of function values ya(1:m,1:n) , tabulated at the grid points defined by x1a andx2a;a n d given values x1andx2of the independent variables; this routine returns an interpolated function value y, and an accuracy indication dy(based only on the interpolation in the x1 direction, however). INTEGER j,kREAL ymtmp(MMAX),yntmp(NMAX)do 12j=1,m Loop over rows. do11k=1,n Copy the row into temporary storage. yntmp(k)=ya(j,k) enddo 11 call polint(x2a,yntmp,n,x2,ymtmp(j),dy) Interpolate answer into temporary stor- age. enddo 12 call polint(x1a,ymtmp,m,x1,y,dy) Do the final interpolation. return END HigherOrderforSmoothness: Bicubic Interpolation We will give two methods that are in common use, and which are themselves not unrelated. The first is usually called bicubic interpolation . Bicubic interpolation requires the user to specify at each grid point not just the function y(x1,x2), but also the gradients ∂y/∂x 1≡y,1,∂y/∂x 2≡y,2and the cross derivative ∂2y/∂x 1∂x 2≡y,12. Then an interpolating function that is cubicin the scaled coordinates tand u(equation 3.6.4) can be found, with the following properties: (i) The values of the function and the speci fied derivatives are reproduced exactly on the grid points, and (ii) the values of the function and thespecifiedderivativeschangecontinuouslyas theinterpolatingpointcrossesfrom one grid square to another. Itisimportanttounderstandthatnothingintheequationsofbicubicinterpolation requiresyoutospecifytheextraderivatives correctly! Thesmoothnesspropertiesare tautologically “forced,”and have nothing to do with the “accuracy”of the speci fied derivatives. It is a separate problem for you to decide how to obtain the values that are specified. The better you do, the more accurate the interpolation will be. But it will be smoothno matter what you do. 3.6InterpolationinTwoor MoreDimensions 119Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Bestofallistoknowthederivativesanalytically,ortobeabletocomputethem accuratelybynumericalmeans,atthegridpoints. Nextbestis todeterminethembynumericaldifferencingfromthefunctionalvaluesalreadytabulatedonthegrid. The relevant code would be something like this (using centered differencing): y1a(j,k)=(ya(j+1,k)-ya(j-1,k))/(x1a(j+1)-x1a(j-1)) y2a(j,k)=(ya(j,k+1)-ya(j,k-1))/(x2a(k+1)-x2a(k-1)) y12a(j,k)=(ya(j+1,k+1)-ya(j+1,k-1)-ya(j-1,k+1)+ya(j-1,k-1)) /((x1a(j+1)-x1a(j-1))*(x2a(k+1)-x2a(k-1))) To do a bicubicinterpolationwithin a gridsquare, giventhe function yand the derivatives y1,y2,y12ateachofthefourcornersofthesquare,therearetwosteps: First obtain the sixteen quantities cij,i , j =1,..., 4using the routine bcucof below. (The formulas that obtain the c’s from the function and derivative values are just a complicated linear transformation, with coef ficients which, having been determined once in the mists of numerical history, can be tabulated and forgotten.)Next,substitutethe c’s intoanyor allof thefollowingbicubicformulasforfunction and derivatives, as desired: y(x 1,x2)=4/summationdisplay i=14/summationdisplay j=1cijti−1uj−1 y,1(x1,x2)=4/summationdisplay i=14/summationdisplay j=1(i−1)cijti−2uj−1(dt/dx 1) y,2(x1,x2)=4/summationdisplay i=14/summationdisplay j=1(j−1)cijti−1uj−2(du/dx 2) y,12(x1,x2)=4/summationdisplay i=14/summationdisplay j=1(i−1)(j−1)cijti−2uj−2(dt/dx 1)(du/dx 2)(3.6.6 ) where tanduare again given by equation (3.6.4). SUBROUTINE bcucof(y,y1,y2,y12,d1,d2,c) REAL d1,d2,c(4,4),y(4),y1(4),y12(4),y2(4) Given arrays y,y1,y2 ,a n d y12, each of length 4, containing the function, gradients, and cross derivative at the four grid points of a rectangular grid cell (numbered counterclockwise from the lower left), and given d1andd2, the length of the grid cell in the 1- and 2- directions, this routine returns the table c(1:4,1:4) that is used by routine bcuint for bicubic interpolation. INTEGER i,j,k,lREAL d1d2,xx,cl(16),wt(16,16),x(16)SAVE wt DATA wt/1,0,-3,2,4*0,-3,0,9,-6,2,0,-6,4,8*0,3,0,-9,6,-2,0,6,-4 * ,10*0,9,-6,2*0,-6,4,2*0,3,-2,6*0,-9,6,2*0,6,-4* ,4*0,1,0,-3,2,-2,0,6,-4,1,0,-3,2,8*0,-1,0,3,-2,1,0,-3,2 * ,10*0,-3,2,2*0,3,-2,6*0,3,-2,2*0,-6,4,2*0,3,-2 * ,0,1,-2,1,5*0,-3,6,-3,0,2,-4,2,9*0,3,-6,3,0,-2,4,-2* ,10*0,-3,3,2*0,2,-2,2*0,-1,1,6*0,3,-3,2*0,-2,2* ,5*0,1,-2,1,0,-2,4,-2,0,1,-2,1,9*0,-1,2,-1,0,1,-2,1 * ,10*0,1,-1,2*0,-1,1,6*0,-1,1,2*0,2,-2,2*0,-1,1/ d1d2=d1*d2do 11i=1,4 Pack a temporary vector x. x(i)=y(i) 120 Chapter3. InterpolationandExtrapolationSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).x(i+4)=y1(i)*d1 x(i+8)=y2(i)*d2 x(i+12)=y12(i)*d1d2 enddo 11 do13i=1,16 Matrix multiply by the stored table. xx=0. do12k=1,16 xx=xx+wt(i,k)*x(k) enddo 12 cl(i)=xx enddo 13 l=0do 15i=1,4 Unpack the result into the output table. do14j=1,4 l=l+1c(i,j)=cl(l) enddo 14 enddo 15 return END Theimplementationofequation(3.6.6),whichperformsabicubicinterpolation, returns the interpolated function value and the two gradient values, and uses theabove routine bcucof, is simply: SUBROUTINE bcuint(y,y1,y2,y12,x1l,x1u,x2l,x2u,x1,x2,ansy, * ansy1,ansy2) REAL ansy,ansy1,ansy2,x1,x1l,x1u,x2,x2l,x2u,y(4),y1(4), * y12(4),y2(4) C USES bcucof Bicubic interpolation within a grid square. Input quantities are y,y1,y2,y12 (as described inbcucof );x1l andx1u, the lower and upper coordinates of the grid square in the 1- direction; x2l andx2u likewise for the 2-direction; and x1,x2 , the coordinates of the desired point for the interpolation. The interpolated function value is returned as ansy , and the interpolated gradient values as ansy1 andansy2 . This routine calls bcucof . INTEGER iREAL t,u,c(4,4) call bcucof(y,y1,y2,y12,x1u-x1l,x2u-x2l,c) Get the c’s. if(x1u.eq.x1l.or.x2u.eq.x2l)pause ’bad input in bcuint’t=(x1-x1l)/(x1u-x1l) Equation (3.6.4). u=(x2-x2l)/(x2u-x2l) ansy=0.ansy2=0.ansy1=0. do 11i=4,1,-1 Equation (3.6.6). ansy=t*ansy+((c(i,4)*u+c(i,3))*u+c(i,2))*u+c(i,1)ansy2=t*ansy2+(3.*c(i,4)*u+2.*c(i,3))*u+c(i,2) ansy1=u*ansy1+(3.*c(4,i)*t+2.*c(3,i))*t+c(2,i) enddo 11 ansy1=ansy1/(x1u-x1l)ansy2=ansy2/(x2u-x2l) return END HigherOrderfor Smoothness: Bicubic Spline The other common technique for obtaining smoothness in two-dimensional interpolation is the bicubic spline . Actually, this is equivalent to a special case 3.6InterpolationinTwoor MoreDimensions 121Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).of bicubic interpolation: The interpolating function is of the same functional form as equation (3.6.6); the values of the derivatives at the grid points are,however, determined “globally”by one-dimensional splines. However, bicubic splines are usually implemented in a form that looks rather different from the above bicubic interpolation routines, instead looking much closer in form to theroutine polin2above: To interpolate one functional value, one performs mone- dimensional splines across the rows of the table, followed by one additional one-dimensional spline down the newly created column. It is a matter of taste (and trade-off between time and memory) as to how much of this process one wants to precompute and store. Instead of precomputing and storing all thederivative information (as in bicubic interpolation), spline users typically precom- pute and store only one auxiliary table, of second derivatives in one direction only. Then one need only do spline evaluations (not constructions) for the m row splines; one must still do a construction andan evaluation for the final col- umn spline. (Recall that a spline construction is a process of order N, while a spline evaluation is only of order logN—and that is just to find the place in the table!) Here is a routine to precomputethe auxiliary second-derivativetable: SUBROUTINE splie2(x1a,x2a,ya,m,n,y2a) INTEGER m,n,NN REAL x1a(m),x2a(n),y2a(m,n),ya(m,n) PARAMETER (NN=100) Maximum expected value of nandm. C USES spline Given an mbyntabulated function ya(1:m,1:n) , and tabulated independent variables x2a(1:n) , this routine constructs one-dimensional natural cubic splines of the rows of ya and returns the second-derivatives in the array y2a(1:m,1:n) . (The array x1a is included in the argument list merely for consistency with routine splin2 .) INTEGER j,k REAL y2tmp(NN),ytmp(NN) do13j=1,m do11k=1,n ytmp(k)=ya(j,k) enddo 11 call spline(x2a,ytmp,n,1.e30,1.e30,y2tmp) Values 1×1030signal a natural spline. do12k=1,n y2a(j,k)=y2tmp(k) enddo 12 enddo 13 returnEND After the above routine has been executed once, any number of bicubic spline interpolationscan be performedby successive calls of the following routine: SUBROUTINE splin2(x1a,x2a,ya,y2a,m,n,x1,x2,y) INTEGER m,n,NN REAL x1,x2,y,x1a(m),x2a(n),y2a(m,n),ya(m,n) PARAMETER (NN=100) Maximum expected value of nand m. C USES spline,splint Given x1a,x2a,ya,m,nas described in splie2 andy2a as produced by that routine; and given a desired interpolating point x1,x2; this routine returns an interpolated function value yby bicubic spline interpolation. INTEGER j,k REAL y2tmp(NN),ytmp(NN),yytmp(NN) 122 Chapter3. InterpolationandExtrapolationSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).do12j=1,m Perform mevaluations of the row splines constructed by splie2 , using the one- dimensional spline evaluator splint .do11k=1,n ytmp(k)=ya(j,k) y2tmp(k)=y2a(j,k) enddo 11 call splint(x2a,ytmp,y2tmp,n,x2,yytmp(j)) enddo 12 call spline(x1a,yytmp,m,1.e30,1.e30,y2tmp) Construct the one-dimensional column spline and evaluate it. call splint(x1a,yytmp,y2tmp,m,x1,y) return END CITED REFERENCES AND FURTHER READING: Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions , Applied Mathe- matics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 byDover Publications, New York), §25.2. Kinahan, B.F., and Harm, R. 1975, Astrophysical Journal , vol. 200, pp. 330–335. Johnson, L.W., and Riess, R.D. 1982, Numerical Analysis , 2nd ed. (Reading, MA: Addison- Wesley), §5.2.7. Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall), §7.7.