[racket] Implementation of Simpson's Rule (SICP Exercise 1.29)

[Ben Duan]

I've asked the question on Stack
Overflow<http://stackoverflow.com/questions/19706893/implementation-of-simpsons-rule-sicp-exercise-1-29>[1].
Óscar López gave me an
answer. But I still don't understand it. So I re-post it here:

Following is my code for SICP exercise 1.29
<http://mitpress.mit.edu/sicp/full-text/book/book-Z-H-12.html>[2]. 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

[1]
http://stackoverflow.com/questions/19706893/implementation-of-simpsons-rule-sicp-exercise-1-29
[2] http://mitpress.mit.edu/sicp/full-text/book/book-Z-H-12.html

Thanks

____________________
  Racket Users list:
  http://lists.racket-lang.org/users