18.3 Integral Equations with Singular Kernels 799 While the terms in brackets superficially appear to scale ask2,there is typically cancellation at both O()and O(). Equation (18.3.5)can be specialized to various choices of (a,b).The obvious choice is a =kh,b=(k+3)h,in which case we get a four-point quadrature rule that generalizes Simpson's 3/8 rule (equation 4.1.5).In fact,we can recover this special case by setting w(x)=1,in which case (18.3.4)becomes n+lk+3)m+1-kn+] (18.3.6) The four terms in square brackets equation (18.3.5)each become independent of k,and (18.3.5)in fact reduces to r(k+3)h kh e)k=3距， +警1+++2物+警r+A(83刀 81 Back to the case of general w(x),some other choices for a and b are also useful.For example,we may want to choose (a,b)to be (k+1h,+3h)or (k+2h,k+3h), 速 allowing us to finish off an extended rule whose number of intervals is not a multiple of three,without loss of accuracy:The integral will be estimated using the four values f(kh),...,f([+3]h).Even more useful is to choose (a,b)to be (+1]h,[+2]h),thus using four points to integrate a centered single interval.These weights,when sewed together RECIPES I into an extended formula,give quadrature schemes that have smooth coefficients,i.e.,without 令 the Simpson-like 2.4,2.4.2 alternation.(In fact,this was the technique that we used to derive equation 4.1.14,which you may now wish to reexamine. All these rules are of the same order as the extended Simpson's rule,that is,exact for f()a cubic polynomial.Rules of lower order,if desired,are similarly obtained.The Press. three point formula is u(f)=f(k+1(k+2Wo-(2k+3w+W +f(k+)-k(k+2)W+2(+1)W1-W (18.3.8) SCIENTIFIC 6 +[k+2kk+1Wo-(2k+1)W+W Here the simple special case is to take,w()=1,so that w=车7k+2+1-k h (18.3.9) Then equation (18.3.8)becomes Simpson's rule, Numerical 色 e)证=专+警rk+)+专k+2 (k+2)h (18.3.10) 43126 For nonconstant weight functions w(),however,equation (18.3.8)gives rules of one order (outside less than Simpson,since they do not benefit from the extra symmetry of the constant case. The two point formula is simply North (k+1)h w(x)f(x)da=f(kh)[(k+1)Wo-W+f(k+1]h)[-Wo+W](18.3.11) /kh Here is a routine wwghts that uses the above formulas to return an extended N-point quadrature rule for the interval (a,b)=(0,[N-1]h).Input to wwghts is a user-supplied routine,kermom,that is called to get the first four indefinite-integral moments of w(),namely 「y Fm(y）= s"w(s)ds m=0,1,2,3 (18.3.12) (The lower limit is arbitrary and can be chosen for convenience.Cautionary note:When called with N<4,wwghts returns a rule of lower order than Simpson;you should structure your problem to avoid this.18.3 Integral Equations with Singular Kernels 799 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). While the terms in brackets superficially appear to scale as k2, there is typically cancellation at both O(k2) and O(k). Equation (18.3.5) can be specialized to various choices of (a, b). The obvious choice is a = kh, b = (k + 3)h, in which case we get a four-point quadrature rule that generalizes Simpson’s 3/8 rule (equation 4.1.5). In fact, we can recover this special case by setting w(x)=1, in which case (18.3.4) becomes Wn = h n + 1 [(k + 3)n+1 − kn+1] (18.3.6) The four terms in square brackets equation (18.3.5) each become independent of k, and (18.3.5) in fact reduces to
(k+3)h kh f(x)dx = 3h 8 f(kh)+ 9h 8 f([k+1]h)+ 9h 8 f([k+2]h)+ 3h 8 f([k+3]h) (18.3.7) Back to the case of general w(x), some other choices for a and b are also useful. For example, we may want to choose (a, b) to be ([k + 1]h, [k + 3]h) or ([k + 2]h, [k + 3]h), allowing us to finish off an extended rule whose number of intervals is not a multiple of three, without loss of accuracy: The integral will be estimated using the four values f(kh),...,f([k + 3]h). Even more useful is to choose (a, b) to be ([k + 1]h, [k + 2]h), thus using four points to integrate a centered single interval. These weights, when sewed together into an extended formula, give quadrature schemes that have smooth coefficients, i.e., without the Simpson-like 2, 4, 2, 4, 2 alternation. (In fact, this was the technique that we used to derive equation 4.1.14, which you may now wish to reexamine.) All these rules are of the same order as the extended Simpson’s rule, that is, exact for f(x) a cubic polynomial. Rules of lower order, if desired, are similarly obtained. The three point formula is
b a w(x)f(x)dx = 1 2 f(kh)  (k + 1)(k + 2)W0 − (2k + 3)W1 + W2  + f([k + 1]h)  − k(k + 2)W0 + 2(k + 1)W1 − W2  + 1 2 f([k + 2]h)  k(k + 1)W0 − (2k + 1)W1 + W2  (18.3.8) Here the simple special case is to take, w(x)=1, so that Wn = h n + 1 [(k + 2)n+1 − kn+1] (18.3.9) Then equation (18.3.8) becomes Simpson’s rule,
(k+2)h kh f(x)dx = h 3 f(kh) + 4h 3 f([k + 1]h) + h 3 f([k + 2]h) (18.3.10) For nonconstant weight functions w(x), however, equation (18.3.8) gives rules of one order less than Simpson, since they do not benefit from the extra symmetry of the constant case. The two point formula is simply
(k+1)h kh w(x)f(x)dx = f(kh)[(k + 1)W0 − W1] + f([k + 1]h)[−kW0 + W1] (18.3.11) Here is a routine wwghts that uses the above formulas to return an extended N-point quadrature rule for the interval (a, b) = (0, [N − 1]h). Input to wwghts is a user-supplied routine, kermom, that is called to get the first four indefinite-integral moments of w(x), namely Fm(y) ≡
y smw(s)ds m = 0, 1, 2, 3 (18.3.12) (The lower limit is arbitrary and can be chosen for convenience.) Cautionary note: When called with N < 4, wwghts returns a rule of lower order than Simpson; you should structure your problem to avoid this