Implementation of Simpson's Rule (SICP Exercise 1.29)

Question

Following is my code for SICP exercise 1.29. The exercise asks us to implement Simpson's Rule using higher order procedure sum. It's supposed to be more accurate than the original integral procedure. But I don't know why it's not the case in my code:

(define (simpson-integral f a b n)
  (define h (/ (- b a) n))
  (define (next x) (+ x (* 2 h)))
  (* (/ h 3) (+ (f a)
                (* 4 (sum f (+ a h) next (- b h)))
                (* 2 (sum f (+ a (* 2 h)) next (- b (* 2 h))))
                (f b))))

Some explanations of my code: As

h/3 * (y_{0} + 4*y_{1} + 2*y_{2} + 4*y_{3} + 2*y_{4} + ... + 2*y_{n-2} + 4*y_{n-1} + y_{n})

equals

h/3 * (y_{0}
       + 4 * (y_{1} + y_{3} + ... + y_{n-1})
       + 2 * (y_{2} + y_{4} + ... + y_{n-2})
       + y_{n})

I just use sum to compute y_{1} + y_{3} + ... + y_{n-1} and y_{2} + y_{4} + ... + y_{n-2}.

Complete code here:

#lang racket

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

(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))

(define (simpson-integral f a b n)
  (define h (/ (- b a) n))
  (define (next x) (+ x (* 2 h)))
  (* (/ h 3) (+ (f a)
                (* 4 (sum f (+ a h) next (- b h)))
                (* 2 (sum f (+ a (* 2 h)) next (- b (* 2 h))))
                (f b))))

Some tests(The exact value should be 0.25):

> (integral cube 0 1 0.01)
0.24998750000000042
> (integral cube 0 1 0.001)
0.249999875000001

> (simpson-integral cube 0 1.0 100)
0.23078806666666699
> (simpson-integral cube 0 1.0 1000)
0.24800798800666748
> (simpson-integral cube 0 1.0 10000)
0.2499999999999509

Answer 1

In your solution the x-values are computed as follows:

h = (b-a)/n
x1 = a+1
x3 = x1 +2*h
x5 = x3 +2*h
...

This means rounding errors slowly accumulate. It happens when (b-a)/n is not representable as floating point.

If we instead compute xi as a+ (i*(b-a))/n you will get more accurate results.

This variant of your solution uses the above method to compute the xi.

(define (simpson-integral3 f a b n)
  (define h (/ (- b a) n))
  (define (next i) (+ i 2))
  (define (f* i) (f (+ a (/ (* i (- b a)) n))))
  (* (/ h 3)
     (+ (f a)
        (* 4 (sum f* 1 next n))
        (* 2 (sum f* 2 next (- n 1)))
        (f b))))

Answer 2

There's a problem in how you're constructing the terms, the way you're alternating between even terms (multiplied by 2) and odd terms (multiplied by 4) is not correct. I solved this problem by passing an additional parameter to sum to keep track of the current term's even-or-odd nature, there are other ways but this worked for me, and the accuracy got improved:

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

(define (simpson-integral f a b n)
  (let* ((h (/ (- b a) n))
         (term (lambda (x i)
                 (if (even? i)
                     (* 2.0 (f x))
                     (* 4.0 (f x)))))
         (next (lambda (x) (+ x h))))
    (* (+ (f a)
          (sum term a next b 1)
          (f b))
       (/ h 3.0))))

(simpson-integral cube 0 1 1000)
=> 0.2510004999999994