124 Chapter 3. Interpolation and Extrapolation pt.number pt④ pt③ ①②③④ desired pt x=x2u y (x1,x2） dy/ax d2 dy/dx2 user pt. ① pt② these 2=x21 dy/dx1dx2 三 dl (a) (b) %时 83g Figure 3.6.1. (a)Labeling of points used in the two-dimensional interpolation routines bcuint and 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. (Nort server function value changes continuously.However,the gradient of the interpolated function changes discontinuously at the boundaries of each grid square. America Press. 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 increased accuracy for the interpolated function (for sufficiently smooth functions!) 9 without necessarily trying to fix up the continuity of the gradient and higher derivatives.Or,one can make use of higher order to enforce smoothness of some of these derivatives as the interpolating point crosses grid-square boundaries.We will 61 now consider each of these two directions in turn. Higher Order for Accuracy The basic idea is to break up the problem into a succession of one-dimensional interpolations.If we want to do m-1 order interpolation in the z direction,and n-1 Numerical 10521 order in the x2 direction,we first locate an m x n sub-block of the tabulated function matrix that contains our desired point(1,z2).We then do m one-dimensional uctio 4310 interpolations in the z2 direction,i.e.,on the rows of the sub-block,to get function values at the points(x1a[j],2),j=1,...,m.Finally,we do a last interpolation (outside Recipes in the direction to get the answer.If we use the polynomial interpolation routine polint of $3.1,and a sub-block which is presumed to be already located (and North Software. addressed through the pointer float **ya,see $1.2),the procedure looks like this: #include "nrutil.h" void polin2(float xia[],float x2a[],float **ya,int m,int n,float x1, float x2,float y,float *dy) Given arrays x1a[1..m]and x2a[1..n]of independent variables,and a submatrix of function values ya[1..m][1..n],tabulated at the grid points defined by x1a and x2a:and given values x1 and x2 of 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). void polint(float xa[],float ya[],int n,float x,float *y,float *dy);124 Chapter 3. Interpolation and Extrapolation Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). y ∂y/∂x1 ∂y/∂x2 ∂2y/∂x1∂x2 x2 = x2u x2 = x2l x1 = x1u x1 = x1l 1234 pt. 1 user supplies these values pt. 4 pt. 2 pt. 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 bcuint and 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. 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 increased accuracy for the interpolated function (for sufficiently smooth functions!), without necessarily trying to fix up the continuity of the gradient and higher derivatives. Or, one can make use of higher order to enforce smoothness of some of these derivatives as the interpolating point crosses grid-square boundaries. We will now consider each of these two directions in turn. Higher Order for Accuracy The basic idea is to break up the problem into a succession of one-dimensional interpolations. If we want to do m-1 order interpolation in the x1 direction, and n-1 order in the x2 direction, we first locate an m × n sub-block of the tabulated function matrix that contains our desired point (x1, x2). We then do m one-dimensional interpolations in the x2 direction, 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 x1 direction to get the answer. If we use the polynomial interpolation routine polint of §3.1, and a sub-block which is presumed to be already located (and addressed through the pointer float **ya, see §1.2), the procedure looks like this: #include "nrutil.h" void polin2(float x1a[], float x2a[], float **ya, int m, int n, float x1, float x2, float *y, float *dy) Given arrays x1a[1..m] and x2a[1..n] of independent variables, and a submatrix of function values ya[1..m][1..n], tabulated at the grid points defined by x1a and x2a; and given values x1 and x2 of 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). { void polint(float xa[], float ya[], int n, float x, float *y, float *dy);