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 <math.h> 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<=JMAX;++){ s[j]=trapzd(func,a,b,j); 1f(0>=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 intervals140 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 machine￾readable 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 <math.h> #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 approxi￾int 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 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