c4-4 (779481), страница 2
Текст из файла (страница 2)
aa andbb must have the same sign.{float x,tnm,sum,del,ddel,b,a;static float s;int it,j;b=1.0/aa;These two statements change the limits of integration.a=1.0/bb;if (n == 1) {From this point on, the routine is identical to midpnt.return (s=(b-a)*FUNC(0.5*(a+b)));} else {for(it=1,j=1;j<n-1;j++) it *= 3;tnm=it;del=(b-a)/(3.0*tnm);ddel=del+del;x=a+0.5*del;sum=0.0;for (j=1;j<=it;j++) {sum += FUNC(x);x += ddel;sum += FUNC(x);x += del;}return (s=(s+(b-a)*sum/tnm)/3.0);}}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).The differences between qromo and qromb (§4.3) are so slight that it is perhapsgratuitous to list qromo in full. It, however, is an excellent driver routine for solvingall the other problems of improper integrals in our first list (except the intractablefifth), as we shall now see.The basic trick for improper integrals is to make a change of variables toeliminate the singularity, or to map an infinite range of integration to a finite one.For example, the identity1454.4 Improper IntegralsIf you need to integrate from a negative lower limit to positive infinity, you dothis by breaking the integral into two pieces at some positive value, for example,answer=qromo(funk,-5.0,2.0,midpnt)+qromo(funk,2.0,1.0e30,midinf);If the singularity is at the upper limit, use the identityZ (b−a)1−γZ bγ11f(x)dx =t 1−γ f(b − t 1−γ )dt1−γ 0a(b > a)(4.4.4)If there is a singularity at both limits, divide the integral at an interior breakpointas in the example above.Equations (4.4.3) and (4.4.4) are particularly simple in the case of inversesquare-root singularities, a case that occurs frequently in practice:Z √b−aZ bf(x)dx =2tf(a + t2 )dt(b > a)(4.4.5)a0for a singularity at a, andZ bZf(x)dx =a√b−a2tf(b − t2 )dt(b > a)(4.4.6)0for a singularity at b.
Once again, we can implement these changes of variabletransparently to the user by defining substitute routines for midpnt which make thechange of variable automatically:#include <math.h>#define FUNC(x) (2.0*(x)*(*funk)(aa+(x)*(x)))float midsql(float (*funk)(float), float aa, float bb, int n)This routine is an exact replacement for midpnt, except that it allows for an inverse square-rootsingularity in the integrand at the lower limit aa.{float x,tnm,sum,del,ddel,a,b;static float s;int it,j;b=sqrt(bb-aa);a=0.0;if (n == 1) {The rest of the routine is exactly like midpnt and is omitted.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).Where should you choose the breakpoint? At a sufficiently large positive value sothat the function funk is at least beginning to approach its asymptotic decrease tozero value at infinity.
The polynomial extrapolation implicit in the second call toqromo deals with a polynomial in 1/x, not in x.To deal with an integral that has an integrable power-law singularity at its lowerlimit, one also makes a change of variable. If the integrand diverges as (x − a)γ ,0 ≤ γ < 1, near x = a, use the identityZ (b−a)1−γZ bγ11f(x)dx =t 1−γ f(t 1−γ + a)dt(b > a)(4.4.3)1−γ 0a146Chapter 4.Integration of FunctionsSimilarly,#include <math.h>#define FUNC(x) (2.0*(x)*(*funk)(bb-(x)*(x)))b=sqrt(bb-aa);a=0.0;if (n == 1) {The rest of the routine is exactly like midpnt and is omitted.One last example should suffice to show how these formulas are derived ingeneral. Suppose the upper limit of integration is infinite, and the integrand falls offexponentially.
Then we want a change of variable that maps e−x dx into (±)dt (withthe sign chosen to keep the upper limit of the new variable larger than the lowerlimit). Doing the integration gives by inspectiont = e−xx = − log tor(4.4.7)so thatZZx=∞t=e−af(x)dx =x=af(− log t)t=0dtt(4.4.8)The user-transparent implementation would be#include <math.h>#define FUNC(x) ((*funk)(-log(x))/(x))float midexp(float (*funk)(float), float aa, float bb, int n)This routine is an exact replacement for midpnt, except that bb is assumed to be infinite(value passed not actually used). It is assumed that the function funk decreases exponentiallyrapidly at infinity.{float x,tnm,sum,del,ddel,a,b;static float s;int it,j;b=exp(-aa);a=0.0;if (n == 1) {The rest of the routine is exactly like midpnt and is omitted.CITED REFERENCES AND FURTHER READING:Acton, F.S. 1970, Numerical Methods That Work; 1990, corrected edition (Washington: Mathematical Association of America), Chapter 4.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).float midsqu(float (*funk)(float), float aa, float bb, int n)This routine is an exact replacement for midpnt, except that it allows for an inverse square-rootsingularity in the integrand at the upper limit bb.{float x,tnm,sum,del,ddel,a,b;static float s;int it,j;4.5 Gaussian Quadratures and Orthogonal Polynomials147Dahlquist, G., and Bjorck, A.
1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall),§7.4.3, p. 294.Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag),§3.7, p. 152.In the formulas of §4.1, the integral of a function was approximated by the sumof its functional values at a set of equally spaced points, multiplied by certain aptlychosen weighting coefficients. We saw that as we allowed ourselves more freedomin choosing the coefficients, we could achieve integration formulas of higher andhigher order. The idea of Gaussian quadratures is to give ourselves the freedom tochoose not only the weighting coefficients, but also the location of the abscissas atwhich the function is to be evaluated: They will no longer be equally spaced.
Thus,we will have twice the number of degrees of freedom at our disposal; it will turn outthat we can achieve Gaussian quadrature formulas whose order is, essentially, twicethat of the Newton-Cotes formula with the same number of function evaluations.Does this sound too good to be true? Well, in a sense it is. The catch is afamiliar one, which cannot be overemphasized: High order is not the same as highaccuracy. High order translates to high accuracy only when the integrand is verysmooth, in the sense of being “well-approximated by a polynomial.”There is, however, one additional feature of Gaussian quadrature formulas thatadds to their usefulness: We can arrange the choice of weights and abscissas to makethe integral exact for a class of integrands “polynomials times some known functionW (x)” rather than for the usual class of integrands “polynomials.” The functionW (x) can then be chosen to remove integrable singularities from the desired integral.Given W (x), in other words, and given an integer N , we can find a set of weightswj and abscissas xj such that the approximationZbW (x)f(x)dx ≈aNXwj f(xj )(4.5.1)j=1is exact if f(x) is a polynomial.
For example, to do the integralZ1−1exp(− cos2 x)√dx1 − x2(4.5.2)(not a very natural looking integral, it must be admitted), we might well be interestedin a Gaussian quadrature formula based on the choiceW (x) = √11 − x2(4.5.3)in the interval (−1, 1). (This particular choice is called Gauss-Chebyshev integration,for reasons that will become clear shortly.)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).4.5 Gaussian Quadratures and OrthogonalPolynomials.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
