A Modification of Filon's Method of Numerical Integration

• Authors: E. A. Flinn
• Venue: Journal of the ACM (JACM)
• Year: 1960

Filon's method of numerical integration was developed to deal with integrals of the form I = ∫BA ƒ(x) cos px dx (1) (Filon, 1928; Tranter, 1951). This method, most useful when p is large, is a modified Simpson's rule using an interval no larger than is required to integrate ∫BA ƒ(x) dx alone to the desired accuracy. The derivation proceeds as follows: The range of integration is divided into panels of width 2h, and a second-order curve is fitted to the middle and end ordinates of one panel. After twice integrating by parts over the width of the panel and summing over all the panels, the result is I = h{&agr;}[ƒ(B) sin pB - ƒ(A) sin pA] + &bgr;Ce + &ggr;Co}, (2) where h is the interval, Ce is the sum of all the even ordinates of ƒ(x) cos px less half the end ordinates, Co is the sum of all the odd ordinates of ƒ(x) cos px less half the end ordinates, Co is the sum of all the odd ordinates of ƒ(x) cos px, and &thgr; = hp &thgr;3&agr; = &thgr;2 + &thgr; sin &thgr; cos &thgr; - 2 sin2 &thgr; &thgr;3&bgr; = 2[&thgr;(1 + cos2) - 2 sin &thgr; cos &thgr;] &thgr;3&ggr; = 4[sin &thgr; - &thgr; cos &thgr;]. (3) In the limit as p approaches zero, (2) reduces to Simpson's rule.The present modification was developed to evaluate functions of the form F(T) = ∫T0 ƒ(x) cos px dx using a larger interval for a permissible error than is possible with Filon's formula. A fifth-order instead of a second-order curve may be fitted to the middle and end points of a panel. Substituting the first five terms of a Stirling approximation into (1), integrating over the width of the panel by parts five times, and summing over all the panels, we obtain I = ∫BA ƒ(x) cos px dx = h{S[ƒ(B) sin pB - ƒ(A) sin pA] + hP[ƒ′(B) cos pB - ƒ′(A) cos pA] + RCce + hQ C′se + NCco + hM C′so (4) where primes denote differentiation with respect to the argument, Cco = sum of odd ordinates of ƒ(x) cos px; C′so = sum of odd ordinates of ƒ′(x) sin px; Cce = sum of even ordinates of ƒ(x) cos px, less half the end ordinates; C′se = sum of even ordinates of ƒ′(x) sin px, less half the end ordinates; and &thgr; = hp &thgr;6M = 16&thgr;(15 - &thgr;2) cos &thgr; + 48(2&thgr;2 - 5) sin &thgr; &thgr;6N = 16&thgr;(3 - &thgr;2) sin &thgr; - 48&thgr;2 cos &thgr; &thgr;6P = 2&thgr;(&thgr;2 - 24) sin &thgr; cos &thgr; + 15(&thgr;2 - 4) cos2 &thgr; + &thgr;4 - 27&thgr;2 + 60 &thgr;6Q = 2[&thgr;(12 - 5&thgr;2) + 15(thgr;2 - 4) sin &thgr; cos &thgr; + 2&thgr;(24 - &thgr;2) cos2 &thgr;] &thgr;2R = 2[&thgr;(156 - 7&thgr;2) sin &thgr; cos &thgr; + 3(60 - 17&thgr;2) cos2 &thgr; - 15(12 - 5&thgr;2)] &thgr;6S = &thgr; (&thgr;4 + 8&thgr;2 - 24) + &thgr;(7&thgr;2 - 156) cos2 &thgr; + 3(60 - 17&thgr;2) sin &thgr; cos &thgr;. (5) For &thgr; less than about 0.9, it is better to expand equations (5) in powers of &thgr;: M = -16/105&thgr; + 8/945&thgr;3 - 2/10395&thgr;5 + 1/405405&thgr;7 - 1/48648600&thgr;9 N = 16/15 - 8/105&thgr;2 + 2/945&thgr;4 - 1/31185&thgr;6 + 1/3243240&thgr;8 P = -1/15 + 2/105&thgr;2 - 1/315&thgr;4 + 2/7425&thgr;6 - 62/4729725&thgr;8 Q = -8/105&thgr; + 16/945&thgr;3 - 104/51975&thgr;5 + 256/2027025&thgr;7 - 16/3274425&thgr;9 R = 14/15 - 16/105&thgr;2 + 22/945&thgr;4 - 608/311850&thgr;6 + 268/2837835&thgr;8 S = 19/105&thgr; - 2/63&thgr;3 + 1/275&thgr;5 - 2/8775&thgr;7 + 34/3869775&thgr;9. (6) For sin px instead of cos px in the integrand, the result is ∫BA ƒ(x) sin px dx= h{S[ƒ(A) cos pA - ƒ (B) cos pB] + hP[ƒ′(B) sin pB - ƒ′(A) sin pA] + RCse - hQC′co + NCso - hMC′co} with an obvious change of notation for the C's.In the limit as p approaches zero, equation (4) reduces to the modified Simpson's rule described by Lanczos (1957).Analytical estimation of the error involved in this method has not been carried out. As an empirical check on the error—and on the expansions (6)—the integral ∫1.50.5 ex cos &pgr;x dx was evaluated to nine decimal places by Filon's method and by the modification. The errors, compared with the true value of -1.7718441, were: Interval Number of points Error, Filon's method Error, Modification 0.1 11 .00000141 -8 0.25 5 .00070660 .00000016 0.5 3 .00051522 .00008785 As Lanczos (1957) points out, the validity of these integration procedures depends on the convergence of the Stirling approximation of ƒ(x). For Filon's method this necessitates the smooth behaviour of differences up to the fourth order, and for the modification presented here, of differences up to the sixth order.Used with Hitchcock's (1957) approximations, Filon's method, or this modification of it, is also useful with integrals containing J0(px) or J1(px) instead of the trigonometric functions.Luke (1954) considered the fitting of an nth order curve to the middle and end points of a panel, and gave a detailed discussion of the error involved—the results obtained here are special cases of his equations. Luke's formulae are, however, suitable primarily for hand computation using tabulated functions, whereas the results given here are intended for use with an electronic computer.The author wishes to thank the referee for drawing his attention to the paper by Luke.