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 <jensa...@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 <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
