To recap, the 3-point Simpson rule divides the interval of integration (x0, x1, x2) into 2 equal panels of width h, and the 3 ordinates at the panel boundaries are combined using simple arithmetic
|S = (h/3) (y0 + 4y1 + y2)|
to approximate the actual integral I. Next, the rule is applied to the two subintervals (x0, x1) and (x1, x2), 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 (x0, x1, x2) into (x0, x0.5, x1) and (x1, x1.5, x2), two new ordinates are needed at x0.5 and x1.5, but the old ordinates at x0, x1, and x2 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
|(2/45) h [7y0 + 32y1 + 12y2 + 32y3 + 7y4]|
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.
|error_term(S) = k * hp|
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,
|error_term(S) = (2p-1)* error_term(S next)|
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
||(2p-1 - 1)* error_term(S next)| must be ≤ (2p-1 - 1)ε|
In other words,
||S next - S| must be ≤ (2p-1 - 1)ε|
For Simpson’s rule, p = 5, so we confirm
||S next - S| must indeed be ≤ 15ε|
For Boole’s rule, p = 7, so
||S next - S| must be ≤ 63ε|
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 entire code is in quadboole.lisp.
*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