SICP Exercise 1.29: Integration using Simpson's Rule

From SICP section 1.3.1 Procedures as Arguments

Exercise 1.29 asks us to define a procedure that uses Simpson's rule to approximate the value of the integral of a function between two values. (Informally, you may remember that the integral of a function is the area under the curve of that function.)

We're given a big head start earlier in the section when an integral procedure is defined using a different method.

(define (sum term a next b)
  (if (> a b)
      0
      (+ (term a)
         (sum term (next a) next b))))

(define (integral f a b dx)
  (define (add-dx x) (+ x dx))
  (* (sum f (+ a (/ dx 2.0)) add-dx b)
     dx))

We're told to use this procedure to check our results by integrating cube between 0 and 1.

Simpson's rule states that the integral of a function f between a and b is approximated as

h(y₀ + 4y₁ + 2y₂ + 4y₃ + 2y₄ + ... + 2y_n-2 + 4y_n-1 + y_n) / 3

where h = (b - a)/n, for some even integer n, and y_k = f(a + kh).

This is the sum of a series, so we'll still be defining our new procedure in terms of the sum procedure used before. Using Simpson's rule, the sum of a series is multiplied by h and divided by 3, so we'll start with that and a little of the "wishful thinking" that Professor Abelson spoke so highly of in lecture 2A.

(define (simpson f a b n)
   (/ (* h (sum term 0 inc n)) 3))

Now that we know how Simpson's rule can be defined in terms of sum, we just need to fill in the pieces that are missing.

The variable h is pretty easy to define from the description.

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

The procedure used to get from one term of the series to the next is even simpler. n is just incremented by one at each step, so we just need to define a procedure to do that.

(define (inc x) (+ x 1))

We know that the sum procedure takes two functions, term and next, and two values a and b, and computes the sum of the terms of the function from a to b. Defining the terms of the series in Simpson's rule is a two-step process. First we have to define the function for computing y_k, which is given.

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

Next we have to define a rule for computing the coefficient for each of the k terms. Once we know the coefficient we'll just multiply it by y_k to get the complete term. The rules for defining the coefficients are pretty simple. Notice that if k is odd, then the coefficient is always 4. If k is even, then the coefficient is usually 2, except for the first (0th) and last (nth) terms, where the coefficient is 1.

(define (term k)
   (* (cond ((odd? k) 4)
            ((or (= k 0) (= k n)) 1)
            ((even? k) 2))
      (y k)))

Putting this all together, we have the complete procedure:

(define (simpson f a b n)
  (define h (/ (- b a) n))
  (define (inc x) (+ x 1))
  (define (y k)
    (f (+ a (* k h))))
  (define (term k)
    (* (cond ((odd? k) 4)
             ((or (= k 0) (= k n)) 1)
             ((even? k) 2))
       (y k)))
  (/ (* h (sum term 0 inc n)) 3))

The only thing left to do is to define a cube procedure and compare the results of simpson with those of the old integral procedure that we were given.

(define (cube x) (* x x x))

> (integral cube 0 1 0.01)
0.24998750000000042
> (simpson cube 0 1 100.0)
0.24999999999999992
> (integral cube 0 1 0.001)
0.249999875000001
> (simpson cube 0 1 1000.0)
0.2500000000000003

As you can see from these results, Simpson's rule gives us a much better approximation to the integral when computing the same number of terms.


Related:

For links to all of the SICP lecture notes and exercises that I've done so far, see The SICP Challenge.