《数字信号处理》教学参考资料（Numerical Recipes in C，The Art of Scientific Computing Second Edition）Chapter 04.2 Integration of Functions 4.2 Elementary Algorithms

团购合买资源类别：文库，文档格式：PDF，文档页数：5，文件大小：70.23KB

136 Chapter 4.Integration of Functions N=1 2 3 ● 4 ● -●(total after N=4) Figure 4.2.1.Sequential calls to the routine trapzd incorporate the information from previous calls and evaluate the integrand only at those new points necessary to refine the grid.The bottom line shows the totality of function evaluations after the fourth call.The routine gsimp,by weighting the intermediate results,transforms the trapezoid rule into Simpson's rule with essentially no additional overhead. There are also formulas of higher order for this situation,but we will refrain from giving them. The semi-open formulas are just the obvious combinations ofequations(4.1.11)- RECIPES (4.1.14)with(4.1.15)-(4.1.18),respectively.At the closed end of the integration, 令 use the weights from the former equations:at the open end use the weights from the latter equations.One example should give the idea,the formula with error term decreasing as 1/N3 which is closed on the right and open on the left: Press. fa=h++f+6+ 23 1 (4.1.20) 13 5 +fN-2+12fN-1+2fN +0(a) 61 CITED REFERENCES AND FURTHER READING: Abramowitz,M.,and Stegun,I.A.1964,Handbook of Mathematical Functions,Applied Mathe- matics Series,Volume 55 (Washington:National Bureau of Standards;reprinted 1968 by Dover Publications,New York),$25.4.[1] Isaacson,E.,and Keller,H.B.1966,Analysis of Numerical Methods(New York:Wiley).$7.1. 10.621 E喜 Numerical Recipes 43106 4.2 Elementary Algorithms (outside Our starting point is equation(4.1.11),the extended trapezoidal rule.There are North Software. two facts about the trapezoidal rule which make it the starting point for a variety of algorithms.One fact is rather obvious,while the second is rather"deep." The obvious fact is that,for a fixed function f(z)to be integrated between fixed limits a and b,one can double the number of intervals in the extended trapezoidal rule without losing the benefit of previous work.The coarsest implementation of the trapezoidal rule is to average the function at its endpoints a and b.The first stage of refinement is to add to this average the value of the function at the halfway point.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:

136 Chapter 4. Integration of Functions Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). N = 1 2 3 4 (total after N = 4) Figure 4.2.1. Sequential calls to the routine trapzd incorporate the information from previous calls and evaluate the integrand only at those new points necessary to refine the grid. The bottom line shows the totality of function evaluations after the fourth call. The routine qsimp, by weighting the intermediate results, transforms the trapezoid rule into Simpson’s rule with essentially no additional overhead. There are also formulas of higher order for this situation, but we will refrain from giving them. The semi-open formulas are just the obvious combinations of equations (4.1.11)– (4.1.14) with (4.1.15)–(4.1.18), respectively. At the closed end of the integration, use the weights from the former equations; at the open end use the weights from the latter equations. One example should give the idea, the formula with error term decreasing as 1/N 3 which is closed on the right and open on the left: xN x1 f(x)dx = h 23 12f2 + 7 12f3 + f4 + f5+ ··· + fN−2 + 13 12fN−1 + 5 12fN + O 1 N3 (4.1.20) 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 by Dover 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 Algorithms Our starting point is equation (4.1.11), the extended trapezoidal rule. There are two facts about the trapezoidal rule which make it the starting point for a variety of algorithms. 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 fixed limits a and b, one can double the number of intervals in the extended trapezoidal rule without losing the benefit of previous work. The coarsest implementation of the trapezoidal rule is to average the function at its endpoints a and b. The first stage of refinement is to add to this average the value of the function at the halfway point. 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:

138 Chapter 4.Integration of Functions Unsophisticated as it is,routine qtrap is in fact a fairly robust way of doing integrals of functions that are not very smooth.Increased sophistication will usually translate into a higher-order method whose efficiency will be greater only for sufficiently smooth integrands.qtrap is the method of choice,e.g.,for an integrand which is a function of a variable that is linearly interpolated between measured data points.Be sure that you do not require too stringent an EPS,however:If qtrap takes too many steps in trying to achieve your required accuracy,accumulated roundoff errors may start increasing,and the routine may never converge.A value 10-6 is just on the edge of trouble for most 32-bit machines;it is achievable when the convergence 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 of order 1/N2,is in fact entirely even when expressed in powers of 1/N.This follows 二分州 directly from the Euler-Maclaurin Summation Formula h=n+++…+N-1+ (4.2.1) 9 --2g--- B2h2 B2kh2k Here B2k is a Bernoulli number,defined by the generating function 9 n =∑B t (4.2.2) IENTIFIC n=0 6 with the first few even values (odd values vanish except for B1=-1/2) 1 B0=1B2=6 B4=一30 1 1 B6=2 (4.2.3) 1 10.621 B8=-30 5 691 B10= 66 B12=-2730 Numerica 431 Equation (4.2.1)is not a convergent expansion,but rather only an asymptotic 、 (outside Recipes expansion whose error when truncated at any point is always less than twice the magnitude of the first neglected term.The reason that it is not convergent is that the Bernoulli numbers become very large,e.g., North Software. B50= 495057205241079648212477525 66 The 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 84.1. For example,equation(4.1.12)has an error series beginning with O(1/N3),but continuing with all subsequent powers of N:1/N4,1/N5,etc. Suppose we evaluate(4.1.11)with N steps,getting a result SN,and then again with 2N steps,getting a result S2N.(This is done by any two consecutive calls of

138 Chapter 4. Integration of Functions Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). Unsophisticated as it is, routine qtrap is in fact a fairly robust way of doing integrals of functions that are not very smooth. Increased sophistication will usually translate into a higher-order method whose efficiency will be greater only for sufficiently smooth integrands. qtrap is the method of choice, e.g., for an integrand which is a function of a variable that is linearly interpolated between measured data points. Be sure that you do not require too stringent an EPS, however: If qtrap takes too many steps in trying to achieve your required accuracy, accumulated roundoff errors may start increasing, and the routine may never converge. A value 10 −6 is just on the edge of trouble for most 32-bit machines; it is achievable when the convergence 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 of order 1/N 2, is in fact entirely even when expressed in powers of 1/N. This follows directly from the Euler-Maclaurin Summation Formula, xN x1 f(x)dx = h 1 2 f1 + f2 + f3 + ··· + fN−1 + 1 2 fN − B2h2 2! (f N − f 1) −···− B2kh2k (2k)! (f(2k−1) N − f(2k−1) 1 ) −··· (4.2.1) Here B2k is a Bernoulli number, defined by the generating function t et − 1 = ∞ n=0 Bn tn n! (4.2.2) with the first few even values (odd values vanish except for B1 = −1/2) B0 = 1 B2 = 1 6 B4 = − 1 30 B6 = 1 42 B8 = − 1 30 B10 = 5 66 B12 = − 691 2730 (4.2.3) Equation (4.2.1) is not a convergent expansion, but rather only an asymptotic expansion whose error when truncated at any point is always less than twice the magnitude of the first neglected term. The reason that it is not convergent is that the Bernoulli numbers become very large, e.g., B50 = 495057205241079648212477525 66 The 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), but continuing with all subsequent powers of N: 1/N 4, 1/N5, etc. Suppose we evaluate (4.1.11) with N steps, getting a result S N , and then again with 2N steps, getting a result S2N . (This is done by any two consecutive calls of

4.2 Elementary Algorithms 139 trapzd)The leading error term in the second evaluation will be 1/4 the size of the error in the first evaluation.Therefore the combination 1 S= 352N-3 (4.2.4) will cancel out the leading order error term.But there is no error term of order 1/N 3, by (4.2.1).The surviving error is of order 1/N4,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 evaluating that rule, and we can write it as a routine exactly analogous to qtrap above: 8 #include #define EPS 1.0e-6 #define JMAX 20 float 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 the desired fractional accuracy and JMAX so that 2 to the power JMAX-1 is the maximum allowed RECIPES I number 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; Ameri computer University Press. one paper THE f1oats,st,0st=0.0,0s=0.0; ART 是 for (j=1;j5) Avoid spurious early convergence. st st if (fabs(s-os)<EPS*fabs(os)II (s==0.0&2os==0.0))return s; 08=S; to dir ost=st; OF SCIENTIFIC COMPUTING(ISBN nrerror("Too many steps in routine qsimp"); return 0.0; Never get here. 18881920 v@cam The routine qsimp will in general be more efficient than qtrap(i.e.,require fewer function evaluations)when the function to be integrated has a finite 4th Numerical Recipes 10-521 43108 derivative (i.e.,a continuous 3rd derivative).The combination of qsimp and its necessary workhorse trapzd is a good one for light-duty work. (outside North Software. CITED REFERENCES AND FURTHER READING: Stoer,J.,and Bulirsch,R.1980,Introduction to Numerical Analysis (New York:Springer-Verlag). 63.3. Dahlquist,G.,and Bjorck,A.1974,Numerica/Methods (Englewood Cliffs,NJ:Prentice-Hall). 887.4.1-7.4.2. Forsythe,G.E.,Malcolm,M.A.,and Moler,C.B.1977,Computer Methods for Mathematical Computations (Englewood Cliffs,NJ:Prentice-Hall),S5.3

4.2 Elementary Algorithms 139 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). trapzd.) The leading error term in the second evaluation will be 1/4 the size of the error in the first evaluation. Therefore the combination S = 4 3 S2N − 1 3 SN (4.2.4) will cancel out the leading order error term. But there is no error term of order 1/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 evaluating that rule, and we can write it as a routine exactly analogous to qtrap above: #include #define EPS 1.0e-6 #define JMAX 20 float 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 the desired fractional accuracy and JMAX so that 2 to the power JMAX-1 is the maximum allowed number 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=0.0,os=0.0; for (j=1;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., require fewer function evaluations) when the function to be integrated has a finite 4th derivative (i.e., a continuous 3rd derivative). The combination of qsimp and its necessary 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 Mathematical Computations (Englewood Cliffs, NJ: Prentice-Hall), §5.3

140 Chapter 4.Integration of Functions 4.3 Romberg Integration We can view Romberg's method as the natural generalization of the routine qsimp in the last section to integration schemes that are of higher order than Simpson's rule.The basic idea is to use the results from k successive refinements of the extended trapezoidal rule (implemented in trapzd)to remove all terms in the error series up to but not including O(1/N2).The routine qsimp is the case of k=2.This is one example of a very general idea that goes by the name of Richardson's deferred approach to the limit:Perform some numerical algorithm for various values of a parameter h,and then extrapolate the result to the continuum 81 limit h 0. Equation(4.2.4),which subtracts off the leading error term,is a special case of polynomial extrapolation.In the more general Romberg case,we can use Neville's 茶 algorithm (see 83.1)to extrapolate the successive refinements to zero stepsize. Neville's algorithm can in fact be coded very concisely within a Romberg integration routine.For clarity of the program.however,it seems better to do the extrapolation RECIPES I by function call to polint,already given in 83.1. 令 #include Press. #define EPS 1.0e-6 #define JMAX 20 #define JMAXP (JMAX+1) #define K 5 Here 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. IENTIFIC float gromb(float (*func)(float),float a,float b) Returns the integral of the function func from a to b.Integration is performed by Romberg's method 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,dssi float s [JMAXP],h[JMAXP+1]; These store the successive trapezoidal approxi- 10621 int j; mations and their relative stepsizes. h[1]=1.0; uurrggroglrion Numerical Recipes 43106 for (j=1;j=K)[ (outside polint(&h[j-K],&s [j-K],K,0.0,&ss,&dss); if (fabs(dss)<=EPS*fabs(ss))return ss; Software. h[j+1]=0.25*h[j]; Ame ying of This is a key step:The factor is 0.25 even though the stepsize is decreased by only 0.5.This makes the extrapolation a polynomial in h2 as allowed by equation (4.2.1), not just a polynomial in h d nrerror("Too many steps in routine qromb"); return 0.0; Never get here. The routine gromb,along with its required trapzd and polint,is quite powerful for sufficiently smooth(e.g.,analytic)integrands,integrated over intervals

140 Chapter 4. Integration of Functions Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 4.3 Romberg Integration We can view Romberg’s method as the natural generalization of the routine qsimp in the last section to integration schemes that are of higher order than Simpson’s rule. The basic idea is to use the results from k successive refinements of the extended trapezoidal rule (implemented in trapzd) to remove all terms in the error series up to but not including O(1/N 2k). The routine qsimp is the case of k = 2. This is one example of a very general idea that goes by the name of Richardson’s deferred approach to the limit: Perform some numerical algorithm for various values of a parameter h, and then extrapolate the result to the continuum limit h = 0. Equation (4.2.4), which subtracts off the leading error term, is a special case of polynomial extrapolation. In the more general Romberg case, we can use Neville’s algorithm (see §3.1) to extrapolate the successive refinements to zero stepsize. Neville’s algorithm can in fact be coded very concisely within a Romberg integration routine. For clarity of the program, however, it seems better to do the extrapolation by function call to polint, already given in §3.1. #include #define EPS 1.0e-6 #define JMAX 20 #define JMAXP (JMAX+1) #define K 5 Here 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’s method 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= 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 only 0.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 quite powerful for sufficiently smooth (e.g., analytic) integrands, integrated over intervals

点击进入文档下载页（PDF格式）

已到末页，全文结束

点击下载（PDF格式）

浏览记录