SICP exercise 1.29

From Drewiki

Problem

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

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 in the text in section 1.3.1.

Solution

Here is the integral procedure as given in the text:

(define (integral f a b dx)
  (define (add-dx x) (+ x dx))
  (* (sum f (+ a (/ dx 2.0)) add-dx b)
     dx))
 
(define (sum term a next b)
  (if (> a b)
      0
      (+ (term a)
         (sum term (next a) next b))))

Here is a procedure which performs numerical integration using Simpson's Rule. Note that the coefficient for each y_k is 1 when k = 0 or k = n, else 4 when k is odd or 2 when k is even. The helper procedure is useful so that the constant h is only evaluated once. (A let special form that defines h would make the function more readable, but the text doesn't introduce let until the next section.)

(define (inc k) (+ k 1))
 
(define (simpsons-rule-integral f a b n)
  (define (helper h)
    (define (y k)
      (f (+ a (* k h))))
    (define (term k)
      (cond ((or (= k 0) (= k n)) (y k))
            ((even? k) (* 2 (y k)))
            (else (* 4 (y k)))))
    (* (sum term 0 inc n)
       (/ h 3)))
  (helper (/ (- b a) n)))

Here is the cube procedure:

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

The results of using integral with cube, using Chicken Scheme 3.1 on a MacBook Pro using OS X 10.5:

(integral cube 0 1 0.01)

Output:

0.2499875

(integral cube 0 1 0.001)

Output:

0.249999875000001

And the results of using simpsons-rule-integral with cube:

(simpsons-rule-integral cube 0 1 100)

Output:

0.25

(simpsons-rule-integral cube 0 1 1000)

Output:

0.25

Note that the exact answer is, in fact, 0.25. This implementation of simpsons-rule-integral converges to the exact answer given a sufficient number of terms and the limited precision of both Chicken Scheme and the hardware.