Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 37
Текст из файла (страница 37)
One example should give the idea, the formula with error termdecreasing as 1/N 3 which is closed on the right and open on the left:& xN723f2 + f3 + f4 + f5 +f(x)dx = h1212x1 (4.1.20)1351+O· · · + fN−2 + fN−1 + fN1212N3CITED 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.4. [1]Isaacson, E., and Keller, H.B. 1966, Analysis of Numerical Methods (New York: Wiley), §7.1.4.2 Elementary AlgorithmsOur starting point is equation (4.1.11), the extended trapezoidal rule. There aretwo facts about the trapezoidal rule which make it the starting point for a variety ofalgorithms. One fact is rather obvious, while the second is rather “deep.”The obvious fact is that, for a fixed function f(x) to be integrated between fixedlimits a and b, one can double the number of intervals in the extended trapezoidalrule without losing the benefit of previous work.
The coarsest implementation ofthe trapezoidal rule is to average the function at its endpoints a and b. The firststage of refinement is to add to this average the value of the function at the halfwaypoint. The second stage of refinement is to add the values at the 1/4 and 3/4 points.And so on (see Figure 4.2.1).Without further ado we can write a routine with this kind of logic to it:4.2 Elementary Algorithms137#define FUNC(x) ((*func)(x))float trapzd(float (*func)(float), float a, float b, int n)This routine computes the nth stage of refinement of an extended trapezoidal rule. func is inputas a pointer to the function to be integrated between limits a and b, also input.
When called with'n=1, the routine returns the crudest estimate of ab f (x)dx. Subsequent calls with n=2,3,...(in that sequential order) will improve the accuracy by adding 2n-2 additional interior points.{float x,tnm,sum,del;static float s;int it,j;if (n == 1) {return (s=0.5*(b-a)*(FUNC(a)+FUNC(b)));} else {for (it=1,j=1;j<n-1;j++) it <<= 1;tnm=it;del=(b-a)/tnm;This is the spacing of the points to be added.x=a+0.5*del;for (sum=0.0,j=1;j<=it;j++,x+=del) sum += FUNC(x);s=0.5*(s+(b-a)*sum/tnm);This replaces s by its refined value.return s;}}The above routine (trapzd) is a workhorse that can be harnessed in severalways.
The simplest and crudest is to integrate a function by the extended trapezoidalrule where you know in advance (we can’t imagine how!) the number of steps youwant. If you want 2M + 1, you can accomplish this by the fragmentfor(j=1;j<=m+1;j++) s=trapzd(func,a,b,j);with the answer returned as s.Much better, of course, is to refine the trapezoidal rule until some specifieddegree of accuracy has been achieved:#include <math.h>#define EPS 1.0e-5#define JMAX 20float qtrap(float (*func)(float), float a, float b)Returns the integral of the function func from a to b.
The parameters EPS can be set to thedesired fractional accuracy and JMAX so that 2 to the power JMAX-1 is the maximum allowednumber of steps. Integration is performed by the trapezoidal rule.{float trapzd(float (*func)(float), float a, float b, int n);void nrerror(char error_text[]);int j;float s,olds;olds = -1.0e30;Any number that is unlikely to be the average of thefor (j=1;j<=JMAX;j++) {function at its endpoints will do here.s=trapzd(func,a,b,j);if (j > 5)Avoid spurious early convergence.if (fabs(s-olds) < EPS*fabs(olds) ||(s == 0.0 && olds == 0.0)) return s;olds=s;}nrerror("Too many steps in routine qtrap");return 0.0;Never get here.}138Chapter 4.Integration of FunctionsUnsophisticated as it is, routine qtrap is in fact a fairly robust way of doingintegrals of functions that are not very smooth. Increased sophistication will usuallytranslate into a higher-order method whose efficiency will be greater only forsufficiently smooth integrands. qtrap is the method of choice, e.g., for an integrandwhich is a function of a variable that is linearly interpolated between measured datapoints.
Be sure that you do not require too stringent an EPS, however: If qtrap takestoo many steps in trying to achieve your required accuracy, accumulated roundofferrors may start increasing, and the routine may never converge. A value 10−6is just on the edge of trouble for most 32-bit machines; it is achievable when theconvergence is moderately rapid, but not otherwise.We come now to the “deep” fact about the extended trapezoidal rule, equation(4.1.11). It is this: The error of the approximation, which begins with a term oforder 1/N 2 , is in fact entirely even when expressed in powers of 1/N . This followsdirectly from the Euler-Maclaurin Summation Formula,&xNx111f(x)dx = h f1 + f2 + f3 + · · · + fN−1 + fN22B2 h2 B2k h2k (2k−1)(2k−1)(fN − f1 ) − · · · −(f−− f1) −···2!(2k)! N(4.2.1)Here B2k is a Bernoulli number, defined by the generating function∞tnt=Bnet − 1 n=0n!(4.2.2)with the first few even values (odd values vanish except for B1 = −1/2)B0 = 1B2 =1B8 = −3016B10B4 = −5=66130B12B6 =142691=−2730(4.2.3)Equation (4.2.1) is not a convergent expansion, but rather only an asymptoticexpansion whose error when truncated at any point is always less than twice themagnitude of the first neglected term.
The reason that it is not convergent is thatthe Bernoulli numbers become very large, e.g.,B50 =49505720524107964821247752566The key point is that only even powers of h occur in the error series of (4.2.1).This fact is not, in general, shared by the higher-order quadrature rules in §4.1.For example, equation (4.1.12) has an error series beginning with O(1/N 3 ), butcontinuing with all subsequent powers of N : 1/N 4 , 1/N 5 , etc.Suppose we evaluate (4.1.11) with N steps, getting a result SN , and then againwith 2N steps, getting a result S2N . (This is done by any two consecutive calls of4.2 Elementary Algorithms139trapzd.) The leading error term in the second evaluation will be 1/4 the size of theerror in the first evaluation.
Therefore the combinationS=14S2N − SN33(4.2.4)will cancel out the leading order error term. But there is no error term of order1/N 3 , by (4.2.1). The surviving error is of order 1/N 4 , the same as Simpson’s rule.In fact, it should not take long for you to see that (4.2.4) is exactly Simpson’s rule(4.1.13), alternating 2/3’s, 4/3’s, and all. This is the preferred method for evaluatingthat rule, and we can write it as a routine exactly analogous to qtrap above:#include <math.h>#define EPS 1.0e-6#define JMAX 20float qsimp(float (*func)(float), float a, float b)Returns the integral of the function func from a to b.
The parameters EPS can be set to thedesired fractional accuracy and JMAX so that 2 to the power JMAX-1 is the maximum allowednumber of steps. Integration is performed by Simpson’s rule.{float trapzd(float (*func)(float), float a, float b, int n);void nrerror(char error_text[]);int j;float s,st,ost,os;ost = os = -1.0e30;for (j=1;j<=JMAX;j++) {st=trapzd(func,a,b,j);s=(4.0*st-ost)/3.0;Compare equation (4.2.4), above.if (j > 5)Avoid spurious early convergence.if (fabs(s-os) < EPS*fabs(os) ||(s == 0.0 && os == 0.0)) return s;os=s;ost=st;}nrerror("Too many steps in routine qsimp");return 0.0;Never get here.}The routine qsimp will in general be more efficient than qtrap (i.e., requirefewer function evaluations) when the function to be integrated has a finite 4thderivative (i.e., a continuous 3rd derivative).
The combination of qsimp and itsnecessary workhorse trapzd is a good one for light-duty work.CITED REFERENCES AND FURTHER READING:Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag),§3.3.Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall),§§7.4.1–7.4.2.Forsythe, G.E., Malcolm, M.A., and Moler, C.B. 1977, Computer Methods for MathematicalComputations (Englewood Cliffs, NJ: Prentice-Hall), §5.3.140Chapter 4.Integration of Functions4.3 Romberg IntegrationWe can view Romberg’s method as the natural generalization of the routineqsimp in the last section to integration schemes that are of higher order thanSimpson’s rule. The basic idea is to use the results from k successive refinementsof the extended trapezoidal rule (implemented in trapzd) to remove all terms inthe error series up to but not including O(1/N 2k ).
The routine qsimp is the caseof k = 2. This is one example of a very general idea that goes by the name ofRichardson’s deferred approach to the limit: Perform some numerical algorithm forvarious values of a parameter h, and then extrapolate the result to the continuumlimit h = 0.Equation (4.2.4), which subtracts off the leading error term, is a special case ofpolynomial extrapolation. In the more general Romberg case, we can use Neville’salgorithm (see §3.1) to extrapolate the successive refinements to zero stepsize.Neville’s algorithm can in fact be coded very concisely within a Romberg integrationroutine.
For clarity of the program, however, it seems better to do the extrapolationby function call to polint, already given in §3.1.#include <math.h>#define EPS 1.0e-6#define JMAX 20#define JMAXP (JMAX+1)#define K 5Here EPS is the fractional accuracy desired, as determined by the extrapolation error estimate;JMAX limits the total number of steps; K is the number of points used in the extrapolation.float qromb(float (*func)(float), float a, float b)Returns the integral of the function func from a to b.
Integration is performed by Romberg’smethod of order 2K, where, e.g., K=2 is Simpson’s rule.{void polint(float xa[], float ya[], int n, float x, float *y, float *dy);float trapzd(float (*func)(float), float a, float b, int n);void nrerror(char error_text[]);float ss,dss;float s[JMAXP],h[JMAXP+1];These store the successive trapezoidal approxiint j;mations and their relative stepsizes.h[1]=1.0;for (j=1;j<=JMAX;j++) {s[j]=trapzd(func,a,b,j);if (j >= K) {polint(&h[j-K],&s[j-K],K,0.0,&ss,&dss);if (fabs(dss) <= EPS*fabs(ss)) return ss;}h[j+1]=0.25*h[j];This is a key step: The factor is 0.25 even though the stepsize is decreased by only0.5.
This makes the extrapolation a polynomial in h2 as allowed by equation (4.2.1),not just a polynomial in h.}nrerror("Too many steps in routine qromb");return 0.0;Never get here.}The routine qromb, along with its required trapzd and polint, is quitepowerful for sufficiently smooth (e.g., analytic) integrands, integrated over intervals4.4 Improper Integrals141which contain no singularities, and where the endpoints are also nonsingular. qromb,in such circumstances, takes many, many fewer function evaluations than either ofthe routines in §4.2.
For example, the integral&2x4 log(x +(x2 + 1)dx0converges (with parameters as shown above) on the very first extrapolation, afterjust 5 calls to trapzd, while qsimp requires 8 calls (8 times as many evaluations ofthe integrand) and qtrap requires 13 calls (making 256 times as many evaluationsof the integrand).CITED REFERENCES AND FURTHER READING:Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag),§§3.4–3.5.Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall),§§7.4.1–7.4.2.Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis, 2nd ed.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
