Using Boole’s rule for numerical quadrature in Lisp

Numerical quadrature is usually programmed using the 2- or 3-point Newton-Cotes formulas, better known as the trapezoidal and Simpson rules respectively. As simple as they are, these rules are computationally convenient, because they let one refine the result by the simple technique of repeated halving of intervals, and they allow reuse of ordinates already calculated.

To recap, the 3-point Simpson rule divides the interval of integration (x₀, x₁, x₂) into 2 equal panels of width h, and the 3 ordinates at the panel boundaries are combined using simple arithmetic

to approximate the actual integral I. Next, the rule is applied to the two subintervals (x₀, x₁) and (x₁, x₂), and the sum of these results, S _next, i.e., the new best value, is checked to see if it is “close enough” to S, the most recent old best value. If so, S _next is returned. If not, further subdivision is called for. This is the so-called “adaptive quadrature”, as the frenzy of subdividing adapts to those sections of the interval that need it most. (Note that when subdividing the interval (x₀, x₁, x₂) into (x₀, x_0.5, x₁) and (x₁, x_1.5, x₂), two new ordinates are needed at x_0.5 and x_1.5, but the old ordinates at x₀, x₁, and x₂ can be reused.)

Newton-Cotes formulas with more points give better results. Simpson’s rule (3-point) has an error term that is proportial to the 5th power of the interval times the 4th derivative of the integrand at some point ξ within the interval. What this means is that if the integrand does not have a 4th derivative (i.e., it is a polynomial of degree ≤ 3), Simpson’s rule will give accurate quadratures at the first try! If one chooses the 5-point Boole’s rule, the error term is proportional to the 7th power of the interval times the 6th derivative of the integrand. So polynomials of degree ≤ 5 would get accurate quadratures at first try. In general, a Newton-Cotes formula has an error term proportional to a power p of the interval, with p being 5 for Simpson, and 7 for Boole.

The ordinate summation for Boole’s rule is

which is no tougher than Simpson’s rule to program. (See Abramowitz & Stegun, Handbook of Mathematical Functions, sec. 25.4.14, Bode’s [sic] rule.)

Furthermore, Boole’s rule starts out with 4 panels, and thus is convenient for interval-halving and ordinate sharing, just like Simpson’s rule (which starts with 2 panels).

A crucial part of any adaptive Newton-Cotes quadrature is determining when the difference between S and S _next, |S _next - S|, is small enough. It stands to reason that the rule at every interval-halving should also halve the error spec ε, so that the subintervals equitably share the burden. However, what does that mean for the |S _next - S| check? Ward & Cheney (Numerical Mathematics and Computing, 2/e, sec. 5.4, An Adaptive Simpson’s Scheme) show how |S _next - S| must be ≤ 15ε for Simpson’s rule. Let’s try the same steps without specifying which Newton-Cotes formula is in use.

Let’s say at a certain step, for an interval of width h, the actual integral is I; the calculated integral is S; and the error spec (for that interval) is ε. As we’ve noted, the error term in S, error_term(S) = I - S is proportial to the p’th power of h.

Now, on halving the interval, we get S _next = the sum of two quadratures, each of interval h/2. So, we have

In other words,

Obviously, | error_term(S _next)| is quite a bit smaller than | error_term(S)|, which is good, but what we really want is for | error_term(S _next)| to shrink to at least ε. Let’s now investigate the difference between S and S _next.

So, if | error_term(S _next)| is to be ≤ ε, then simple multiplication gives

In other words,

For Simpson’s rule, p = 5, so we confirm

For Boole’s rule, p = 7, so

One final wrinkle is due to the fact that in the course of our subdividing, we may go past the smallest representable non-zero float in our Lisp. Thus, if we find that either h or ε has decayed to machine 0, we will return whatever result we got at that point.

The code

Boole’s rule requires the ordinates to be finite. For some improper integrals where the integrand diverges at either or both ends, one can program a wrapper that makes the integral more palatable. The wrapper changes the variable of integration from x to u, where x = log u (if the lower limit is - ∞) or x = - log u (if the upper limit is ∞).

The entire code is in quadboole.lisp.

Compare quadsimpson.lisp.

Examples

All of these use a *quadrature-error* (i.e. the initial ε) of 1e-12, i.e., 11 correct places after the decimal.

? (quadrature (lambda (x) (/ (+ 1 (* x x)))) 0 1) ; = π/4
0.7853981633974341

? (/ pi 4)
0.7853981633974483

? (defun normpdf (x)
(/ (exp (* -1/2 x x))
   #.(sqrt (* 2L0 pi))))
NORMPDF

? (quadrature #'normpdf 0 1)
0.3413447460685443

;octave -q --eval 'printf("%.20f", normcdf(1) - 1/2)' gives
;0.34134474606854292578

? (quadrature (lambda (x) (exp (- x))) 0 nil)
1.0

? (quadrature #'exp nil 0)
1.0

? (defun gamma (n)
(cond ((< n 1) (/ (gamma (+ n 1)) n))
      ((= n 1) 1)
      ((< 1 n 2)
       (let ((n-1 (- n 1)))
         (quadrature (lambda (x)
                       (* (expt x n-1)
                          (exp (- x))))
                     0 nil)))
      (t (let ((n-1 (- n 1)))
           (* n-1 (gamma n-1))))))
GAMMA

? (gamma 3/2) ; = (1/2)√π
0.8862269254527402

? (* 1/2 (sqrt pi))
0.8862269254527579