142 Chapter 4.Integration of Functions having an error series that is entirely even in h.Indeed there is a formula,not as well known as it ought to be,called the Second Euler-Maclaurin summation formula, N f(x)dx=h[f3/2+f5/2+fn/2+·+fN-3/2+fN-1/2] B2 4(f-)+ (4.4.1) 2T1-22+1f2-1)-f2-+ B2kh2k This equation can be derived by writing out(4.2.1)with stepsize h,then writing it out again with stepsize h/2,then subtracting the first from twice the second 18881992 It is not possible to double the number of steps in the extended midpoint rule -00 g and still have the benefit of previous function evaluations(try it!).However,it is possible to triple the number of steps and do so.Shall we do this,or double and from NUMERICAL RECIPES I accept the loss?On the average,tripling does a factor v3 of unnecessary work, since the "right"number of steps for a desired accuracy criterion may in fact fall server anywhere in the logarithmic interval implied by tripling.For doubling,the factor (Nort is only v2,but we lose an extra factor of 2 in being unable to use all the previous THE evaluations.Since 1.732 <2 x 1.414,it is better to triple. Here is the resulting routine,which is directly comparable to trapzd. ART #define FUNC(x)((*func)(x)) Progra float midpnt(float (*func)(float),float a,float b,int n) This routine computes the nth stage of refinement of an extended midpoint rule.func is input as 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 off()dr.Subsequent calls with n=2.3... to dir (in that sequential order)will improve the accuracy of s by adding(2/3)x-1additional interior points.s should not be modified between sequential calls. OF SCIENTIFIC COMPUTING(ISBN float x,tnm,sum,del,ddel; 198918920 static float s; int it,ji v@cam 10-621 1f(n==1){ return (s=(b-a)*FUNC(0.5*(a+b))); .Further reproduction, Numerical Recipes 43108 else for(1t=1,j=1;j<n-1;j++)it*=3; tnmsit; del=(b-a)/(3.0*tnm); (outside ddel=del+del; The added points alternate in spacing between Software. x=a+0.5*de1; del and ddel. sum=0.0; for (j=1;j<=it;j++){ sum +FUNC(x); x +ddel; sum +FUNC(x) x +del; s=(s+(b-a)*sum/tnm)/3.0; The new sum is combined with the old integral return s; to give a refined integral.142 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). having an error series that is entirely even in h. Indeed there is a formula, not as well known as it ought to be, called the Second Euler-Maclaurin summation formula,
xN x1 f(x)dx = h[f3/2 + f5/2 + f7/2 + ··· + fN−3/2 + fN−1/2] + B2h2 4 (f N − f 1) + ··· + B2kh2k (2k)! (1 − 2−2k+1)(f(2k−1) N − f(2k−1) 1 ) + ··· (4.4.1) This equation can be derived by writing out (4.2.1) with stepsize h, then writing it out again with stepsize h/2, then subtracting the first from twice the second. It is not possible to double the number of steps in the extended midpoint rule and still have the benefit of previous function evaluations (try it!). However, it is possible to triple the number of steps and do so. Shall we do this, or double and accept the loss? On the average, tripling does a factor √ 3 of unnecessary work, since the “right” number of steps for a desired accuracy criterion may in fact fall anywhere in the logarithmic interval implied by tripling. For doubling, the factor is only √ 2, but we lose an extra factor of 2 in being unable to use all the previous evaluations. Since 1.732 < 2 × 1.414, it is better to triple. Here is the resulting routine, which is directly comparable to trapzd. #define FUNC(x) ((*func)(x)) float midpnt(float (*func)(float), float a, float b, int n) This routine computes the nth stage of refinement of an extended midpoint rule. func is input as 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 b a f(x)dx. Subsequent calls with n=2,3,... (in that sequential order) will improve the accuracy of s by adding (2/3) × 3n-1 additional interior points. s should not be modified between sequential calls. { float x,tnm,sum,del,ddel; static float s; int it,j; if (n == 1) { 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; The added points alternate in spacing between x=a+0.5*del; del and ddel. sum=0.0; for (j=1;j<=it;j++) { sum += FUNC(x); x += ddel; sum += FUNC(x); x += del; } s=(s+(b-a)*sum/tnm)/3.0; The new sum is combined with the old integral return s; to give a refined integral. } }