Sunday, April 25, 2010

SICP Exercise 1.29: Integration using Simpson's Rule

From SICP section 1.3.1 Procedures as Arguments

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(y0 + 4y1 + 2y2 + 4y3 + 2y4 + ... + 2yn-2 + 4yn-1 + yn) / 3

where h = (b - a)/n, for some even integer n, and yk = 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 yk, 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 yk 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:

pbewig said...

I solved this problem at my blog Programming Praxis.

frenchy64 said...

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))

frenchy64 said...

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.

Bill the Lizard said...

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!

frenchy64 said...

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)

Bill the Lizard said...

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.

Chris Bowdon said...

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)))))))

Chris Bowdon said...

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

Anonymous said...

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.

Bill the Lizard said...

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.