SICP Exercise 1.29: Approximating Integrals with Simpson's Rule

Simpson's Rule is a more accurate method of numerical integration than the method illustrated above. Using Simpson's Rule, the integral of a function f between a and b is approximated as

  h
  - [y0 + 4y1 + 2y2 + 4y3 + 2y4 + ... + 2yn-2 + 4yn-1 + yn]
  3

where h = (b - a)/n, for some even integer n, and y_k = f(a + kh). (Increasing n increases the accuracy of the approximation.) Define a procedure that takes as arguments f, a, b, and n and returns the value of the integral, computed using Simpson's Rule. Use your procedure to integrate cube between 0 and 1 (with n = 100 and n = 1000), and compare the results to those of the integral procedure shown above.

Okay, so we want a procedure that takes parameters f, a, b and n and uses sum to perform the sum portion of the equation above. So it's going to have the general form:

(define (simpsons-rule f a b n)
  <definitions of h and kth-term go here>
  (* (sum kth-term 0 inc n)
     (/ h 3.0)))

Producing a definition of h is straightforward. This is a constant value in the context of a particular set of values for a, b and n, and so can be defined simply as:

(define h (/ (- b a) n))

We can also produce a straightforward definition of y_k:

(define (y k) (f (+ a (* k h))))

However to produce the sum we need to know what to multiply y_k by! We need a procedure that multiplies the results of y_k by 1, 2 or 4 depending upon whether k is 0 or n, is even, or is odd respectively:

(define (kth-multiple k)
  (cond ((or (= k 0) (= k n)) 1.0)
        ((even? k) 2.0)
        (else 4.0)))

Using these we can produce kth-term as:

(define (kth-term k)
  (* (kth-multiple k)
     (kth-term k)))

And finally we can put all these together to produce:

(define (simpsons-rule f a b n)
  (define h (/ (- b a) n))
  (define (y k) (f (+ a (* k h))))
  (define (kth-multiple k)
    (cond ((or (= k 0) (= k n)) 1.0)
          ((even? k) 2.0)
          (else 4.0)))
  (define (kth-term k)
    (* (kth-multiple k) (y k)))
  (* (sum kth-term 0 inc n)
     (/ h 3.0)))

Now we could compress this a bit more. We can combine y, kth-multiple and kth-term to produce:

(define (simpsons-rule f a b n)
  (define h (/ (- b a) n))
  (define (kth-term k)
     (* (f (+ a (* k h)))
        (cond ((or (= k 0) (= k n)) 1.0)
              ((even? k) 2.0)
              (else 4.0))))
  (* (sum kth-term 0 inc n)
     (/ h 3.0)))

This is slightly more compact and perhaps a bit more readable.

Right, let's compare Simpson's rule with integral:

> (integral cube 0 1 0.01)
0.24998750000000042
> (simpsons-rule cube 0 1 100)
0.24999999999999992
> (integral cube 0 1 0.001)
0.249999875000001
> (simpsons-rule cube 0 1 1000)
0.2500000000000002

As noted right at the end of section 1.3.1, the exact value of the integral of cube between 0 and 1 is 1/4. As we can see, Simpson's rule provides a more accurate approximation than integral and also increasing n does indeed improve the accuracy of the approximation.