Exercise 1.29 asks us to define a procedure that uses Simpson's rule to approximate the value of the integral of a function between two values. (Informally, you may remember that the integral of a function is the area under the curve of that function.)

We're given a big head start earlier in the section when an integral procedure is defined using a different method.

`(define (sum term a next b)`

(if (> a b)

0

(+ (term a)

(sum term (next a) next b))))

(define (integral f a b dx)

(define (add-dx x) (+ x dx))

(* (sum f (+ a (/ dx 2.0)) add-dx b)

dx))

We're told to use this procedure to check our results by integrating

`cube`

between 0 and 1.Simpson's rule states that the integral of a function f between a and b is approximated as

h(y

_{0}+ 4y_{1}+ 2y_{2}+ 4y_{3}+ 2y_{4}+ ... + 2y_{n-2}+ 4y_{n-1}+ y_{n}) / 3where h = (b - a)/n, for some even integer n, and y

_{k}= f(a + kh).

This is the sum of a series, so we'll still be defining our new procedure in terms of the

`sum`

procedure used before. Using Simpson's rule, the sum of a series is multiplied by h and divided by 3, so we'll start with that and a little of the "wishful thinking" that Professor Abelson spoke so highly of in lecture 2A.`(define (simpson f a b n)`

(/ (* h (sum term 0 inc n)) 3))

Now that we know how Simpson's rule can be defined in terms of

`sum`

, we just need to fill in the pieces that are missing.The variable h is pretty easy to define from the description.

`(define h (/ (- b a) n))`

The procedure used to get from one term of the series to the next is even simpler. n is just incremented by one at each step, so we just need to define a procedure to do that.

`(define (inc x) (+ x 1))`

We know that the

`sum`

procedure takes two functions, `term`

and `next`

, and two values `a`

and `b`

, and computes the sum of the terms of the function from `a`

to `b`

. Defining the terms of the series in Simpson's rule is a two-step process. First we have to define the function for computing y_{k}, which is given.

`(define (y k)`

(f (+ a (* k h))))

Next we have to define a rule for computing the coefficient for each of the k terms. Once we know the coefficient we'll just multiply it by y

_{k}to get the complete term. The rules for defining the coefficients are pretty simple. Notice that if k is odd, then the coefficient is always 4. If k is even, then the coefficient is usually 2, except for the first (0th) and last (nth) terms, where the coefficient is 1.

`(define (term k)`

(* (cond ((odd? k) 4)

((or (= k 0) (= k n)) 1)

((even? k) 2))

(y k)))

Putting this all together, we have the complete procedure:

`(define (simpson f a b n)`

(define h (/ (- b a) n))

(define (inc x) (+ x 1))

(define (y k)

(f (+ a (* k h))))

(define (term k)

(* (cond ((odd? k) 4)

((or (= k 0) (= k n)) 1)

((even? k) 2))

(y k)))

(/ (* h (sum term 0 inc n)) 3))

The only thing left to do is to define a

`cube`

procedure and compare the results of `simpson`

with those of the old `integral`

procedure that we were given.`(define (cube x) (* x x x))`

> (integral cube 0 1 0.01)

0.24998750000000042

> (simpson cube 0 1 100.0)

0.24999999999999992

> (integral cube 0 1 0.001)

0.249999875000001

> (simpson cube 0 1 1000.0)

0.2500000000000003

As you can see from these results, Simpson's rule gives us a much better approximation to the integral when computing the same number of terms.

Related:

For links to all of the SICP lecture notes and exercises that I've done so far, see The SICP Challenge.

## 10 comments:

I solved this problem at my blog Programming Praxis.

Hi Bill, Thanks for doing this series, it's been a fantastic help!

On the last line of the simpson function, shouldn't it be

(* (/ h 3)

(sum term 0 inc n))

Whoops, don't mind me, just brushing up on some high school algebra... ;) Just had a "doh" moment.

I guess my alternative seems more natural to me.

frenchy64,

You're right, that would give exactly the same result, and might be a little more clear. The way I wrote it came out of the way I wrote the formula for Simpson's rule earlier in the article:

h(y0 + 4y1 + 2y2 + 4y3 + 2y4 + ... + 2yn-2 + 4yn-1 + yn) / 3

That was out of pure laziness on my part. The original Wikipedia article where I got the formula has it as h/3[y0...] and I just didn't want to steal their image or create my own. :)

I'm really glad to learn that someone is reading these articles closely enough to catch details like this. I should be back to update the series very soon. Thanks for reading!

I think anyone doing SICP seriously will read this series closely too!

This exercise was an interesting one, because initially I didn't see a way of defining 'h'. Later I realized I had forgotten to apply applicative-order evaluation.

So

(define h (/ (- b a) n))

becomes

(define h (/ (- 10 1) 10))

(define h 0.9)

as it is immediately evaluated at the first recursive call.

Theoretically, if we were to use Normal-Order evaluation, 'a' 'b' and 'n' would refer to different free variables in each recursive depth, right?

So at depth 1, h would be substituted:

(define h (/ (- 10 2) 10))

This is an interesting side effect of normal-order eval, I don't think it has been mentioned in the book so far.

(I'm trying to use the correct terminology, feel free to correct any misunderstandings I have)

frenchy64,

That's right, with applicative order evaluation 'h' is assigned a value right away and it never changes.

Things get a little bit twisty when you try to figure out what happens when you use normal order evaluation. I don't think 'b' and 'n' ever do change, but 'a' certainly does. The entire function body for 'h' gets passed in to 'y', which is used in term, which is evaluated in the recursive calls to sum. The different value of 'a' at each depth should result in a different value of 'h' and on up the chain until you get a different result as well.

I don't think these exact side-effects are mentioned in the text yet either. The only time I can remember the text saying that the result of a procedure can be different under the two different evaluation rules was in exercise 1.5.

Hi Bill,

Great post! And the others on SICP, massively useful!

There's also another way to do this one. Instead of having a procedure to determine the coefficient, you can just split the sum into two, a sum of the odd terms and a sum of the even terms. :)

(define (simpsons-rule f a b n)

(define h (/ (- b a) n))

(define (inc x) (+ x (* 2 h)))

(* (/ h 3) (+ a

(+ a (* n h))

(* 2 (sum f (inc a) inc (+ a (* (- n 2) h))))

(* 4 (sum f (+ a h) inc (+ a (* (- n 1) h)))))))

PS. Guess that way uses 2x space though? Coefficient procedure is probably most efficient.

Bill, I am trying to use Dr Scheme to follow this book and lectures. I keep running into differences between LISP and Scheme syntax (such as nested defines do not seem to be legal, etc.). Do you have any tips? Which (lisp) environment is best for following this book? I am running latest Ubuntu.

Anonymous,

I don't have any experience with Lisp myself, but for solutions to the SICP exercises in Lisp you can check out Eli Bendersky's Website. Eli wrote most of his solutions in Common Lisp.

Post a Comment