144 Chapter 4.Integration of Functions The differences between gromo and gromb(84.3)are so slight that it is perhaps gratuitous to list gromo in full.It,however,is an excellent driver routine for solving all the other problems of improper integrals in our first list (except the intractable fifth),as we shall now see. The basic trick for improper integrals is to make a change of variables to eliminate the singularity.or to map an infinite range of integration to a finite one. For example,the identity f(x)dr= 1/a1 ab 0 (4.4.2) 83 can be used with either b-oo and a positive,or with a--oo and b negative,and 鱼 18881892 works for any function which decreases towards infinity faster than 1/z2 100 You can make the change of variable implied by (4.4.2)either analytically and then use (e.g.)gromo and midpnt to do the numerical evaluation,or you can let from NUMERICAL RECIPES I the numerical algorithm make the change of variable for you. We prefer the latter method as being more transparent to the user.To implement equation (4.4.2)we simply write a modified version of midpnt,called midinf,which allows b to be (Nort server infinite (or,more precisely,a very large number on your particular machine,such e University Press. THE as 1 x 1030),or a to be negative and infinite. America computer, to make one paper ART #define FUNC(x)((funk)(1.0/(x))/((x)*(x))) Effects the change of variable 9 Programs float midinf(float (*funk)(float),float aa,float bb,int n) This routine is an exact replacement for midpnt,i.e.,returns the nth stage of refinement of the integral of funk from aa to bb,except that the function is evaluated at evenly spaced points in 1/rather than in This allows the upper limit bb to be as large and positive as , the computer allows,or the lower limit aa to be as large and negative,but not both.aa and bb 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 e% 1881892 OF SCIENTIFIC COMPUTING(ISBN Numerical Recipes 10-621 a=1.0/bb: 1f(n==1){ From this point on,the routine is identical to midpnt. return (s=(b-a)*FUNC(0.5*(a+b))); 43108 else for(1t=1,j=1;j<n-1;j++)1t*=3; tnm=it; (outside del=(b-a)/(3.0*tnm): ddel=del+del: x=a+0.5*de1; sum=0.0; North Amer Software. for(j=1;j<=1t:j+){ sum +FUNC(x); x +ddel; sum +FUNC(x) x +del; return (s=(s+(b-a)*sum/tnm)/3.0);144 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). The differences between qromo and qromb (§4.3) are so slight that it is perhaps gratuitous to list qromo in full. It, however, is an excellent driver routine for solving all the other problems of improper integrals in our first list (except the intractable fifth), as we shall now see. The basic trick for improper integrals is to make a change of variables to eliminate the singularity, or to map an infinite range of integration to a finite one. For example, the identity
b a f(x)dx =
1/a 1/b 1 t2 f 1 t  dt ab > 0 (4.4.2) can be used with either b → ∞ and a positive, or with a → −∞ and b negative, and works for any function which decreases towards infinity faster than 1/x2. You can make the change of variable implied by (4.4.2) either analytically and then use (e.g.) qromo and midpnt to do the numerical evaluation, or you can let the numerical algorithm make the change of variable for you. We prefer the latter method as being more transparent to the user. To implement equation (4.4.2) we simply write a modified version of midpnt, called midinf, which allows b to be infinite (or, more precisely, a very large number on your particular machine, such as 1 × 1030), or a to be negative and infinite. #define FUNC(x) ((*funk)(1.0/(x))/((x)*(x))) Effects the change of variable. float midinf(float (*funk)(float), float aa, float bb, int n) This routine is an exact replacement for midpnt, i.e., returns the nth stage of refinement of the integral of funk from aa to bb, except that the function is evaluated at evenly spaced points in 1/x rather than in x. This allows the upper limit bb to be as large and positive as the computer allows, or the lower limit aa to be as large and negative, but not both. aa and bb 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); } }