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

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