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

Thank you, Jens. I didn't know that the inexactness of floating point
numbers could make such a big difference.



On Tue, Nov 5, 2013 at 10:01 PM, Jens Axel Søgaard <jensaxel at soegaard.net>wrote:

> I think I got it now.
>
> 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))))
>
>
>
> 2013/11/5 Ben Duan <yfefyf at gmail.com>:
> > Thanks, Jens and Daniel.
> >
> > There's no problem with Simpson's Rule. I got another way to implement
> > Simpson's Rule from Eli Bendersky. And re-wrote it in Racket:
> >
> >     (define (simpson-integral2 f a b n)
> >       (define h (/ (- b a) n))
> >       (define (next counter) (+ counter 1))
> >       (define (term counter)
> >         (* (f (+ a (* h counter))) (cond
> >                                      ((= counter 0) 1)
> >                                      ((= counter n) 1)
> >                                      ((odd? counter) 4)
> >                                      ((even? counter) 2))))
> >       (* (/ h 3) (sum term 0 next n)))
> >
> > The result is very accurate:
> >
> >     > (simpson-integral2 cube 0 1.0 10)
> >     0.25
> >     > (simpson-integral2 cube 0 1.0 100)
> >     0.24999999999999992
> >     > (simpson-integral2 cube 0 1.0 1000)
> >     0.2500000000000003
> >     > (simpson-integral2 cube 0 1.0 10000)
> >     0.2500000000000011
> >
> > Eli's code and mine are two implementations for the same Rule, but mine
> is
> > far less accurate. So I think there's some wrong with my implementation.
> But
> > I still can't find it out.
> >
> >
> > On Tue, Nov 5, 2013 at 1:59 AM, Jens Axel Søgaard <jensaxel at soegaard.net
> >
> > wrote:
> >>
> >> 2013/11/4 Ben Duan <yfefyf at gmail.com>:
> >> > I've asked the question on Stack Overflow [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 [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:
> >>
> >> A couple of points:
> >>
> >>  - to compare to methods M1 and M2 we need to specify the class of
> >> functions,
> >>    we will use the methods on
> >>  - when comparing the results of the two methods on a single function f,
> >>    we must use the same number of evaluations of f
> >>  - the most accurate method for a given number of evaluations n,
> >>    is the method which has the best worst case behaviour
> >>
> >> Thus when SICP says:
> >> "Simpson's Rule is a more accurate method of numerical integration
> >> than the method illustrated above."
> >> it doesn't mean that the Simpson rule will work better than, say, the
> >> midpoint method for all
> >> function. Some methods give very good results for polynomials of a small
> >> degree.
> >> Try your methods on some non-polynomials: say 1/x , ln(x), sqrt(x) and
> >> sin(x).
> >>
> >> See also some interesting examples here:
> >>
> >>
> http://www.math.uconn.edu/~heffernan/math1132s13/files/trapezoid_rule.pdf
> >>
> >>
> >> --
> >> Jens Axel Søgaard
> >
> >
>
>
>
> --
> --
> Jens Axel Søgaard
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.racket-lang.org/users/archive/attachments/20131106/91123934/attachment.html>