Integral using Simpson's Rule

Question

As an answer to this problem:

Exercise 1.29

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 / 3) * (y_0 + 4y_1 + 2y_2 + 4y_3 + 2y_4 + ... + 2y_(n-2) + 4y_(n-1) + y_n

where h = (b - a)/n, for some even integer n, and yk = 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.

I wrote the following solution:

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

(define (simpsons-rule f a b n)
  (let ((h (/ (- b a) n)))
    (define (y_k k) (f (+ a (* k h))))
    (define (even n) (= (remainder n 2) 0))
    (define (term n) (* (if (even n) 2.0 4.0) (y_k n)))
    (define (next n) (+ n 1))
    (* (/ h 3.0) (+ (y_k 0.0) (sum term 0.0 next (- n 1.0)) (y_k n)))))

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

What do you think?

Answer 1

When you call sum, make sure you start with 1 (and not 0). Otherwise, you get an error of + 2 y_0 h / 3. In case of (simpsons-rule cube 1 2 1000), this error is 0.000666....

Another way to rewrite the series is to group even and odd terms together, excluding the first and last terms.

(h / 3) * (y_0 + y_n + 4 * (y_1 + y_3 + ... + y_n-1) + 2 * (y_2 + y_4 + ... + y_n-2))

This gives us another possible implementation:

(define (simpsons-rule f a b n)
  (define h (/ (- b a) n))
  (define (y-k k) (f (+ a (* k h))))
  (define (next n) (+ n 2))
  (* (/ h 3.0) (+ (y-k 0) (y-k n) (* 4 (sum y-k 1 next (- n 1))) (* 2 (sum y-k 2 next (- n 2))))))