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

Interestingly, your code gives the expected result with exact arithmetic:

> (simpson-integral cube 0 1 100)
1/4

and less of an error for a smaller number of steps:

> (simpson-integral cube 0 1.0 10)
0.24999999999999994

Suggestion: Modify the code to print out the intermediate steps in the
calculation and compare inexact arithmetic with exact in detail for the n =
100 / h = 0.01 case.


Dan



On Mon, Nov 4, 2013 at 7:14 AM, Ben Duan <yfefyf at gmail.com> wrote:

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


-- 
*Daniel Prager*
Agile/Lean Coaching, Software Development and Leadership
Twitter: @agilejitsu <https://twitter.com/agilejitsu>
Blog: agile-jitsu.blogspot.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.racket-lang.org/users/archive/attachments/20131104/e6a03bad/attachment.html>