Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 35
Текст из файла (страница 35)
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;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 an128Chapter 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: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;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.Chapter 4.Integration of Functions4.0 IntroductionNumerical integration, which is also called quadrature, has a history extendingback to the invention of calculus and before. The fact that integrals of elementaryfunctions could not, in general, be computed analytically, while derivatives couldbe, served to give the field a certain panache, and to set it a cut above the arithmeticdrudgery of numerical analysis during the whole of the 18th and 19th centuries.With the invention of automatic computing, quadrature became just one numerical task among many, and not a very interesting one at that.
Automatic computing,even the most primitive sort involving desk calculators and rooms full of “computers”(that were, until the 1950s, people rather than machines), opened to feasibility themuch richer field of numerical integration of differential equations. Quadrature ismerely the simplest special case: The evaluation of the integral&bf(x)dxI=(4.0.1)ais precisely equivalent to solving for the value I ≡ y(b) the differential equationdy= f(x)dx(4.0.2)y(a) = 0(4.0.3)with the boundary conditionChapter 16 of this book deals with the numerical integration of differentialequations.
In that chapter, much emphasis is given to the concept of “variable” or“adaptive” choices of stepsize. We will not, therefore, develop that material here.If the function that you propose to integrate is sharply concentrated in one or morepeaks, or if its shape is not readily characterized by a single length-scale, then itis likely that you should cast the problem in the form of (4.0.2)–(4.0.3) and usethe methods of Chapter 16.The quadrature methods in this chapter are based, in one way or another, on theobvious device of adding up the value of the integrand at a sequence of abscissaswithin the range of integration.
The game is to obtain the integral as accuratelyas possible with the smallest number of function evaluations of the integrand. Justas in the case of interpolation (Chapter 3), one has the freedom to choose methods129130Chapter 4.Integration of Functionsof various orders, with higher order sometimes, but not always, giving higheraccuracy.
“Romberg integration,” which is discussed in §4.3, is a general formalismfor making use of integration methods of a variety of different orders, and werecommend it highly.Apart from the methods of this chapter and of Chapter 16, there are yetother methods for obtaining integrals. One important class is based on functionapproximation. We discuss explicitly the integration of functions by Chebyshevapproximation (“Clenshaw-Curtis” quadrature) in §5.9. Although not explicitlydiscussed here, you ought to be able to figure out how to do cubic spline quadratureusing the output of the routine spline in §3.3. (Hint: Integrate equation 3.3.3over x analytically. See [1].)Some integrals related to Fourier transforms can be calculated using the fastFourier transform (FFT) algorithm.
This is discussed in §13.9.Multidimensional integrals are another whole multidimensional bag of worms.Section 4.6 is an introductory discussion in this chapter; the important technique ofMonte-Carlo integration is treated in Chapter 7.CITED REFERENCES AND FURTHER READING:Carnahan, B., Luther, H.A., and Wilkes, J.O. 1969, Applied Numerical Methods (New York:Wiley), Chapter 2.Isaacson, E., and Keller, H.B. 1966, Analysis of Numerical Methods (New York: Wiley), Chapter 7.Acton, F.S. 1970, Numerical Methods That Work; 1990, corrected edition (Washington: Mathematical Association of America), Chapter 4.Stoer, J., and Bulirsch, R.
1980, Introduction to Numerical Analysis (New York: Springer-Verlag),Chapter 3.Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis, 2nd ed. (New York:McGraw-Hill), Chapter 4.Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall),§7.4.Kahaner, D., Moler, C., and Nash, S. 1989, Numerical Methods and Software (Englewood Cliffs,NJ: Prentice Hall), Chapter 5.Forsythe, G.E., Malcolm, M.A., and Moler, C.B. 1977, Computer Methods for MathematicalComputations (Englewood Cliffs, NJ: Prentice-Hall), §5.2, p. 89.
[1]Davis, P., and Rabinowitz, P. 1984, Methods of Numerical Integration, 2nd ed. (Orlando, FL:Academic Press).4.1 Classical Formulas for Equally SpacedAbscissasWhere would any book on numerical analysis be without Mr. Simpson and his“rule”? The classical formulas for integrating a function whose value is known atequally spaced steps have a certain elegance about them, and they are redolent withhistorical association. Through them, the modern numerical analyst communes withthe spirits of his or her predecessors back across the centuries, as far as the timeof Newton, if not farther.
Alas, times do change; with the exception of two of themost modest formulas (“extended trapezoidal rule,” equation 4.1.11, and “extended1314.1 Classical Formulas for Equally Spaced Abscissashx0x1x2xNxN + 1open formulas use these pointsclosed formulas use these pointsFigure 4.1.1. Quadrature formulas with equally spaced abscissas compute the integral of a functionbetween x0 and xN +1 . Closed formulas evaluate the function on the boundary points, while openformulas refrain from doing so (useful if the evaluation algorithm breaks down on the boundary points).midpoint rule,” equation 4.1.19, see §4.2), the classical formulas are almost entirelyuseless.
They are museum pieces, but beautiful ones.Some notation: We have a sequence of abscissas, denoted x0 , x1 , . . . , xN ,xN+1 which are spaced apart by a constant step h,xi = x0 + ihi = 0, 1, . . . , N + 1(4.1.1)A function f(x) has known values at the xi ’s,f(xi ) ≡ fi(4.1.2)We want to integrate the function f(x) between a lower limit a and an upper limitb, where a and b are each equal to one or the other of the xi ’s.
An integrationformula that uses the value of the function at the endpoints, f(a) or f(b), is calleda closed formula. Occasionally, we want to integrate a function whose value at oneor both endpoints is difficult to compute (e.g., the computation of f goes to a limitof zero over zero there, or worse yet has an integrable singularity there). In thiscase we want an open formula, which estimates the integral using only xi ’s strictlybetween a and b (see Figure 4.1.1).The basic building blocks of the classical formulas are rules for integrating afunction over a small number of intervals.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
