Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 34
Текст из файла (страница 34)
, xn )from an n-dimensional grid of tabulated values y and n one-dimensional vectors 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 notCartesian, i.e., has tabulated function values at “random” points in n-dimensionalspace rather than at the vertices of a rectangular array. For clarity, we will considerexplicitly only the case of two dimensions, the cases of three or more dimensionsbeing analogous in every way.In two dimensions, we imagine that we are given a matrix of functional valuesya[1..m][1..n]. We are also given an array x1a[1..m], and an array x2a[1..n].The relation of these input quantities to an underlying function y(x1 , x2 ) isya[j][k] = y(x1a[j], x2a[k])(3.6.1)We want to estimate, by interpolation, the function y at 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.
Forconvenience, we will number these points from 1 to 4, counterclockwise startingfrom the lower left (see Figure 3.6.1). More precisely, ifx1a[j] ≤ x1 ≤ x1a[j+1]x2a[k] ≤ x2 ≤ x2a[k+1](3.6.2)defines j and k, theny1 ≡ ya[j][k]y2 ≡ ya[j+1][k]y3 ≡ ya[j+1][k+1](3.6.3)y4 ≡ ya[j][k+1]The simplest interpolation in two dimensions is bilinear interpolation on thegrid 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 t and u each lie between 0 and 1), andy(x1 , x2) = (1 − t)(1 − u)y1 + t(1 − u)y2 + tuy3 + (1 − t)uy4(3.6.5)Bilinear interpolation is frequently “close enough for government work.” Asthe interpolating point wanders from grid square to grid square, the interpolated124Chapter 3.Interpolation and Extrapolationpt.
number⊗pt. 3desired pt.(x1, x2)x2 = x2upt. 1pt. 234∂y / ∂x2∂2y / ∂x1∂x2x1 = x1ux1 = x1lx2 = x2l2y∂y / ∂x1d2d11usth er ses upev pal liesuespt. 4(a)(b)Figure 3.6.1. (a) Labeling of points used in the two-dimensional interpolation routines bcuint andbcucof. (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 interpolatedfunction changes discontinuously at the boundaries of each grid square.There are two distinctly different directions that one can take in going beyondbilinear interpolation to higher-order methods: One can use higher order to obtainincreased accuracy for the interpolated function (for sufficiently smooth functions!),without necessarily trying to fix up the continuity of the gradient and higherderivatives.
Or, one can make use of higher order to enforce smoothness of some ofthese derivatives as the interpolating point crosses grid-square boundaries. We willnow consider each of these two directions in turn.Higher Order for AccuracyThe basic idea is to break up the problem into a succession of one-dimensionalinterpolations.
If we want to do m-1 order interpolation in the x1 direction, and n-1order in the x2 direction, we first locate an m × n sub-block of the tabulated functionmatrix that contains our desired point (x1 , x2 ). We then do m one-dimensionalinterpolations in the x2 direction, i.e., on the rows of the sub-block, to get functionvalues at the points (x1a[j], x2 ), j = 1, . . . , m. Finally, we do a last interpolationin the x1 direction to get the answer.
If we use the polynomial interpolation routinepolint of §3.1, and a sub-block which is presumed to be already located (andaddressed 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 functionvalues ya[1..m][1..n], tabulated at the grid points defined by x1a and x2a; and given valuesx1 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);3.6 Interpolation in Two or More Dimensions125int j;float *ymtmp;ymtmp=vector(1,m);for (j=1;j<=m;j++) {polint(x2a,ya[j],n,x2,&ymtmp[j],dy);}polint(x1a,ymtmp,m,x1,y,dy);free_vector(ymtmp,1,m);Loop over rows.Interpolate answer into temporary storage.Do the final interpolation.}Higher Order for Smoothness: Bicubic InterpolationWe will give two methods that are in common use, and which are themselvesnot unrelated.
The first is usually called bicubic interpolation.Bicubic interpolation requires the user to specify at each grid point not justthe function y(x1 , x2), but also the gradients ∂y/∂x1 ≡ y,1 , ∂y/∂x2 ≡ y,2 andthe cross derivative ∂ 2 y/∂x1 ∂x2 ≡ y,12 .
Then an interpolating function that iscubic in the scaled coordinates t and u (equation 3.6.4) can be found, with thefollowing properties: (i) The values of the function and the specified derivativesare reproduced exactly on the grid points, and (ii) the values of the function andthe specified derivatives change continuously as the interpolating point crosses fromone grid square to another.It is important to understand that nothing in the equations of bicubic interpolationrequires you to specify the extra derivatives correctly! The smoothness properties aretautologically “forced,” and have nothing to do with the “accuracy” of the specifiedderivatives.
It is a separate problem for you to decide how to obtain the values thatare specified. The better you do, the more accurate the interpolation will be. Butit will be smooth no matter what you do.Best of all is to know the derivatives analytically, or to be able to compute themaccurately by numerical means, at the grid points. Next best is to determine them bynumerical differencing from the functional values already tabulated on the grid. Therelevant 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 thederivatives 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 bcucofbelow. (The formulas that obtain the c’s from the function and derivative valuesare just a complicated linear transformation, with coefficients which, having beendetermined 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 functionand derivatives, as desired:126Chapter 3.y(x1 , x2 ) =44 Interpolation and Extrapolationcij ti−1 uj−1i=1 j=1y,1 (x1 , x2 ) =44 (i − 1)cij ti−2 uj−1 (dt/dx1)i=1 j=144 y,2 (x1 , x2 ) =(j − 1)cij ti−1 uj−2 (du/dx2)(3.6.6)i=1 j=1y,12 (x1 , x2 ) =4 4(i − 1)(j − 1)cij ti−2 uj−2 (dt/dx1)(du/dx2 )i=1 j=1where t and u are again given by equation (3.6.4).void bcucof(float y[], float y1[], float y2[], float y12[], float d1, float d2,float **c)Given arrays y[1..4], y1[1..4], y2[1..4], and y12[1..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 d1 and d2, the length of the grid cell in the 1- and2-directions, this routine returns the table c[1..4][1..4] that is used by routine bcuintfor bicubic interpolation.{static int wt[16][16]={ 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,-3,0,0,3,0,0,0,0,-2,0,0,-1,0,0,0,0,2,0,0,-2,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,-3,0,0,3,0,0,0,0,-2,0,0,-1,0,0,0,0,2,0,0,-2,0,0,0,0,1,0,0,1,-3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0,9,-9,9,-9,6,3,-3,-6,6,-6,-3,3,4,2,1,2,-6,6,-6,6,-4,-2,2,4,-3,3,3,-3,-2,-1,-1,-2,2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0,-6,6,-6,6,-3,-3,3,3,-4,4,2,-2,-2,-2,-1,-1,4,-4,4,-4,2,2,-2,-2,2,-2,-2,2,1,1,1,1};int l,k,j,i;float xx,d1d2,cl[16],x[16];d1d2=d1*d2;for (i=1;i<=4;i++) {Pack a temporary vector x.x[i-1]=y[i];x[i+3]=y1[i]*d1;x[i+7]=y2[i]*d2;x[i+11]=y12[i]*d1d2;}for (i=0;i<=15;i++) {Matrix multiply by the stored table.xx=0.0;for (k=0;k<=15;k++) xx += wt[i][k]*x[k];cl[i]=xx;}l=0;for (i=1;i<=4;i++)Unpack the result into the output table.for (j=1;j<=4;j++) c[i][j]=cl[l++];}3.6 Interpolation in Two or More Dimensions127The implementation of equation (3.6.6), which performs a bicubic interpolation,gives back the interpolated function value and the two gradient values, and uses theabove routine bcucof, is simply:#include "nrutil.h"void bcuint(float y[], float y1[], float y2[], float y12[], float x1l,float x1u, float x2l, float x2u, float x1, float x2, float *ansy,float *ansy1, float *ansy2)Bicubic interpolation within a grid square.
Input quantities are y,y1,y2,y12 (as described inbcucof); x1l and x1u, the lower and upper coordinates of the grid square in the 1-direction;x2l and x2u likewise for the 2-direction; and x1,x2, the coordinates of the desired point forthe interpolation. The interpolated function value is returned as ansy, and the interpolatedgradient values as ansy1 and ansy2.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
