<div dir="ltr"><div>Thank you, Jens. I didn't know that the inexactness of floating point numbers could make such a big difference.</div><div> </div></div><div class="gmail_extra"> <div class="gmail_quote">On Tue, Nov 5, 2013 at 10:01 PM, Jens Axel Søgaard <<a href="mailto:jensaxel@soegaard.net" target="_blank">jensaxel@soegaard.net</a>> wrote: <blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">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) <div class="im"> (define h (/ (- b a) n)) </div> (define (next i) (+ i 2)) (define (f* i) (f (+ a (/ (* i (- b a)) n)))) <div class="im"> (* (/ h 3) (+ (f a) </div> (* 4 (sum f* 1 next n)) (* 2 (sum f* 2 next (- n 1))) (f b)))) 2013/11/5 Ben Duan <<a href="mailto:yfefyf@gmail.com">yfefyf@gmail.com</a>>: <div class="HOEnZb"><div class="h5">> 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 <<a href="mailto:jensaxel@soegaard.net">jensaxel@soegaard.net</a>> > wrote: >> >> 2013/11/4 Ben Duan <<a href="mailto:yfefyf@gmail.com">yfefyf@gmail.com</a>>: >> > 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: >> >> <a href="http://www.math.uconn.edu/~heffernan/math1132s13/files/trapezoid_rule.pdf" target="_blank">http://www.math.uconn.edu/~heffernan/math1132s13/files/trapezoid_rule.pdf</a> >> >> >> -- >> Jens Axel Søgaard > > </div></div>-- -- Jens Axel Søgaard </blockquote></div> </div>