3.6 Interpolation in Two or More Dimensions 125 int ji float *ymtmp; ymtmp=vector(1,m); for(j=1;j<=m;j++)[ Loop over rows. polint(x2a,ya[j],n,x2,&ymtmp[j],dy); Interpolate answer into temporary stor- age. polint(x1a,ymtmp,m,x1,y,dy); Do the final interpolation free_vector(ymtmp,1,m); Higher Order for Smoothness:Bicubic Interpolation nted for 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(r1,2),but also the gradients oy/ox1 =y.1,oy/ox2 =y.2 and 令 the cross derivative 02y/ox1ox2=y.12.Then an interpolating function that is cubic in the scaled coordinates t and u (equation 3.6.4)can be found,with the following properties:(i)The values of the function and the specified derivatives are reproduced exactly on the grid points,and(ii)the values of the function and the specified derivatives change continuously as the interpolating point crosses from one grid square to another. N岂g合 It is important to understand that nothing in the equations of bicubic interpolation requires you to specify the extra derivatives correctly!The smoothness properties are tautologically "forced,"and have nothing to do with the"accuracy"of the specified 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 smooth no matter what you do. Best of all is to know the derivatives analytically,or to be able to compute them accurately by numerical means,at the grid points.Next best is to determine them by numerical differencing from the functional values already tabulated on the grid.The 10621 relevant code would be something like this (using centered differencing): Numerical Recipes 43106 y1a[j][k]=(ya[j+1][]-ya[j-1]k])/(x1a[j+1]-x1a[j-1]); (outside y2a[j][k]=(ya[j][k+1]-ya[j][k-1])/(x2ak+1]-x2a[k-1]); y12a[j][]=(ya[j+1][k+1]-ya[j+1][k-1]-ya[j-1][k+1]+ya[j-1][k-1]) Software. /((x1a[j+1]-x1a[j-1])*(x2ak+1]-x2a[k-1])): To do a bicubic interpolation within a grid square,given the function y and the derivatives y1,y2,y12 at each of the four corners of the square,there are two steps: First obtain the sixteen quantities ci,i,j=1,...,4 using the routine bcucof below.(The formulas that obtain the c's from the function and derivative values are just a complicated linear transformation,with coefficients which,having been determined once in the mists of numerical history,can be tabulated and forgotten.) Next,substitute the c's into any or all of the following bicubic formulas for function and derivatives,as desired:3.6 Interpolation in Two or More Dimensions 125 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). int j; float *ymtmp; ymtmp=vector(1,m); for (j=1;j<=m;j++) { Loop over rows. polint(x2a,ya[j],n,x2,&ymtmp[j],dy); Interpolate answer into temporary stor- } age. polint(x1a,ymtmp,m,x1,y,dy); Do the final interpolation. free_vector(ymtmp,1,m); } Higher Order for Smoothness: 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/∂x1 ≡ y,1, ∂y/∂x2 ≡ y,2 and the cross derivative ∂2y/∂x1∂x2 ≡ y,12. Then an interpolating function that is cubic in the scaled coordinates t and u (equation 3.6.4) can be found, with the following properties: (i) The values of the function and the specified derivatives are reproduced exactly on the grid points, and (ii) the values of the function and the specified derivatives change continuously as the interpolating point crosses from one grid square to another. It is important to understand that nothing in the equations of bicubic interpolation requires you to specify the extra derivatives correctly! The smoothness properties are tautologically “forced,” and have nothing to do with the “accuracy” of the specified 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 smooth no matter what you do. Best of all is to know the derivatives analytically, or to be able to compute them accurately by numerical means, at the grid points. Next best is to determine them by numerical differencing from the functional values already tabulated on the grid. 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 bicubic interpolation within a grid square, given the function y and the derivatives y1, y2, y12 at each of the four corners of the square, there are two steps: First obtain the sixteen quantities cij , i, j = 1,..., 4 using the routine bcucof below. (The formulas that obtain the c’s from the function and derivative values are just a complicated linear transformation, with coefficients which, having been determined once in the mists of numerical history, can be tabulated and forgotten.) Next, substitute the c’s into any or all of the following bicubic formulas for function and derivatives, as desired: