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 <yfe...@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 <jensa...@soegaard.net>
> wrote:
>>
>> 2013/11/4 Ben Duan <yfe...@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
____________________
Racket Users list:
http://lists.racket-lang.org/users
