c3-6 (779476), страница 2
Текст из файла (страница 2)
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).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);126Chapter 3.y(x1 , x2 ) =44 XXInterpolation and Extrapolationcij ti−1 uj−1i=1 j=1y,1 (x1 , x2 ) =44 XX(i − 1)cij ti−2 uj−1 (dt/dx1)i=1 j=1(3.6.6)i=1 j=1y,12 (x1 , x2 ) =4 X4X(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++];}Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).44 XXy,2 (x1 , x2 ) =(j − 1)cij ti−1 uj−2 (du/dx2)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"c=matrix(1,4,1,4);d1=x1u-x1l;d2=x2u-x2l;bcucof(y,y1,y2,y12,d1,d2,c);Get the c’s.if (x1u == x1l || x2u == x2l) nrerror("Bad input in routine bcuint");t=(x1-x1l)/d1;Equation (3.6.4).u=(x2-x2l)/d2;*ansy=(*ansy2)=(*ansy1)=0.0;for (i=4;i>=1;i--) {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.0*c[i][4]*u+2.0*c[i][3])*u+c[i][2];*ansy1=u*(*ansy1)+(3.0*c[4][i]*t+2.0*c[3][i])*t+c[2][i];}*ansy1 /= d1;*ansy2 /= d2;free_matrix(c,1,4,1,4);}Higher Order for Smoothness: Bicubic SplineThe other common technique for obtaining smoothness in two-dimensionalinterpolation is the bicubic spline.
Actually, this is equivalent to a special caseof bicubic interpolation: The interpolating function is of the same functional formas equation (3.6.6); the values of the derivatives at the grid points are, however,determined “globally” by one-dimensional splines. However, bicubic splines areusually implemented in a form that looks rather different from the above bicubicinterpolation routines, instead looking much closer in form to the routine polin2above: To interpolate one functional value, one performs m one-dimensional splinesacross the rows of the table, followed by one additional one-dimensional splinedown the newly created column. It is a matter of taste (and trade-off between timeand memory) as to how much of this process one wants to precompute and store.Instead of precomputing and storing all the derivative information (as in bicubicinterpolation), spline users typically precompute 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 and anSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).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. This routine calls bcucof.{void bcucof(float y[], float y1[], float y2[], float y12[], float d1,float d2, float **c);int i;float t,u,d1,d2,**c;128Chapter 3.Interpolation and Extrapolationevaluation for the final column spline. (Recall that a spline construction is a processof order N , while a spline evaluation is only of order log N — and that is just tofind the place in the table!)Here is a routine to precompute the auxiliary second-derivative table:for (j=1;j<=m;j++)spline(x2a,ya[j],n,1.0e30,1.0e30,y2a[j]);}Values 1×1030 signal a natural spline.(If you want to interpolate on a sub-block of a bigger matrix, see §1.2.)After the above routine has been executed once, any number of bicubic splineinterpolations can be performed by successive calls of the following routine:#include "nrutil.h"void splin2(float x1a[], float x2a[], float **ya, float **y2a, int m, int n,float x1, float x2, float *y)Given x1a, x2a, ya, m, n as described in splie2 and y2a as produced by that routine; andgiven a desired interpolating point x1,x2; this routine returns an interpolated function value yby bicubic spline interpolation.{void spline(float x[], float y[], int n, float yp1, float ypn, float y2[]);void splint(float xa[], float ya[], float y2a[], int n, float x, float *y);int j;float *ytmp,*yytmp;ytmp=vector(1,m);yytmp=vector(1,m);Perform m evaluations of the row splines constructed byfor (j=1;j<=m;j++)splie2, using the one-dimensional spline evaluatorsplint(x2a,ya[j],y2a[j],n,x2,&yytmp[j]);splint.spline(x1a,yytmp,m,1.0e30,1.0e30,ytmp);Construct the one-dimensional colsplint(x1a,yytmp,ytmp,m,x1,y);umn spline and evaluate it.free_vector(yytmp,1,m);free_vector(ytmp,1,m);}CITED REFERENCES AND FURTHER READING:Abramowitz, M., and Stegun, I.A.
1964, Handbook of Mathematical Functions, Applied Mathematics 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: AddisonWesley), §5.2.7.Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall),§7.7.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).void splie2(float x1a[], float x2a[], float **ya, int m, int n, float **y2a)Given an m by n tabulated function ya[1..m][1..n], and tabulated independent variablesx2a[1..n], this routine constructs one-dimensional natural cubic splines of the rows of yaand returns the second-derivatives in the array y2a[1..m][1..n]. (The array x1a[1..m] isincluded in the argument list merely for consistency with routine splin2.){void spline(float x[], float y[], int n, float yp1, float ypn, float y2[]);int j;.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
