4.4 Improper Integrals 143 The routine midpnt can exactly replace trapzd in a driver routine like qtrap (84.2);one simply changes trapzd(func,a,b,j)to midpnt(func,a,b,j),and perhaps also decreases the parameter JMAX since 3JMAX-1(from step tripling)is a much larger number than 2JMAX-1(step doubling). The open formula implementation analogous to Simpson's rule(qsimp in 84.2) substitutes midpnt for trapzd and decreases JMAX as above,but now also changes the extrapolation step to be s=(9.0*st-0st)/8.0; since,when the number of steps is tripled,the error decreases to 1/9th its size,not 1/4th as with step doubling. Either the modified qtrap or the modified gsimp will fix the first problem nted for on the list at the beginning of this section.Yet more sophisticated is to generalize Romberg integration in like manner: #include <math.h> #define EPS 1.0e-6 #define JMAX 14 #define JMAXP (JMAX+1) #define K 5 America computer, float qromo(float (*func)(float),float a,float b, ART float (*choose)(float(*)(float),float,float,int)) Romberg integration on an open interval.Returns the integral of the function func from a to b. using any specified integrating function choose and Romberg's method.Normally choose will Progra be an open formula,not evaluating the function at the endpoints.It is assumed that choose triples the number of steps on each call,and that its error series contains only even powers of the number of steps.The routines midpnt,midinf,midsql,midsqu,midexp,are possible choices for choose.The parameters have the same meaning as in gromb. void polint(float xa[],float ya[],int n,float x,float *y,float *dy); void nrerror(char error_text[]); int ji OF SCIENTIFIC COMPUTING(ISBN float ss,dss,h[JMAXP+1],s[JMAXP] 1881892 h[1]=1.0; for (j=1;j<=JMAX;j++){ 10-621 s[j]=(*choose)(func,a,b,j); 1f(j>=K)[ polint(&h[j-K],&s[j-K],K,0.0,&ss,&dss); Numerical Recipes 43108 if (fabs(dss)<=EPS*fabs(ss))return ss; h[j+1]=h[j]/9.0 This is where the assumption of step tripling and an even (outside error series is used. nrerror("Too many steps in routing gromo"); Software. return 0.0; Never get here. North Ame Don't be put off by qromo's complicated ANSI declaration.A typical invocation (integrating the Bessel function Yo(z)from 0 to 2)is simply #include "nr.h" float answer; answer=gromo(bessy0,0.0,2.0,midpnt);4.4 Improper Integrals 143 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). The routine midpnt can exactly replace trapzd in a driver routine like qtrap (§4.2); one simply changes trapzd(func,a,b,j) to midpnt(func,a,b, j), and perhaps also decreases the parameter JMAX since 3JMAX−1 (from step tripling) is a much larger number than 2JMAX−1 (step doubling). The open formula implementation analogous to Simpson’s rule (qsimp in §4.2) substitutes midpnt for trapzd and decreases JMAX as above, but now also changes the extrapolation step to be s=(9.0*st-ost)/8.0; since, when the number of steps is tripled, the error decreases to 1/9th its size, not 1/4th as with step doubling. Either the modified qtrap or the modified qsimp will fix the first problem on the list at the beginning of this section. Yet more sophisticated is to generalize Romberg integration in like manner: #include <math.h> #define EPS 1.0e-6 #define JMAX 14 #define JMAXP (JMAX+1) #define K 5 float qromo(float (*func)(float), float a, float b, float (*choose)(float(*)(float), float, float, int)) Romberg integration on an open interval. Returns the integral of the function func from a to b, using any specified integrating function choose and Romberg’s method. Normally choose will be an open formula, not evaluating the function at the endpoints. It is assumed that choose triples the number of steps on each call, and that its error series contains only even powers of the number of steps. The routines midpnt, midinf, midsql, midsqu, midexp, are possible choices for choose. The parameters have the same meaning as in qromb. { void polint(float xa[], float ya[], int n, float x, float *y, float *dy); void nrerror(char error_text[]); int j; float ss,dss,h[JMAXP+1],s[JMAXP]; h[1]=1.0; for (j=1;j<=JMAX;j++) { s[j]=(*choose)(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]=h[j]/9.0; This is where the assumption of step tripling and an even } error series is used. nrerror("Too many steps in routing qromo"); return 0.0; Never get here. } Don’t be put off by qromo’s complicated ANSI declaration. A typical invocation (integrating the Bessel function Y0(x) from 0 to 2) is simply #include "nr.h" float answer; ... answer=qromo(bessy0,0.0,2.0,midpnt);