Pages

Monday, April 4, 2011

Exercise 1.29: Simpson's Rule

Use the sum procedure from the chapter to estimate the definite integral of a procedure based on Simpson's rule.

Here's the summation procedure from the chapter:
(define (sum term a next b)
  (if (> a b)
      0
      (+ (term a)
         (sum term (next a) next b))))

My function definition should look like: (define (simpsons f a b n) ... where f is the function to integrate (defined as a procedure), a is the upper bound, b is the lower bound, and n is the number of subdivisions on the interval [a, b].

To get it straight in my head, I rewrote Simpson's rule in a way that was more clearly applicable for this problem.


Writing it down this way made two things clearer:
  1. The argument of the function in the series changes in a straightforward manner. I just need to add h to get the next argument.
  2. The procedure applied is different in each term. Sometimes it's f, sometimes it's 2f, and sometimes it's 4f
First, just to make things easier on myself, I'll define h.
(define h (/ (- b a) n))
Based on the way the sum procedure is defined, I know I need to pass it two procedures. I'll call them simp-term and simp-next.

As mentioned above, it's easy to see how the the function argument changes from term to term. I define simp-next like this.

(define (simp-next x)
    (+ x h))

Figuring out how to apply one of three different procedures to the terms of the summation is a little more puzzling. I start out with some pseudocode:

(define (simp-term x)
    (decide on a procedure)(apply the right procedure to x))

How do I decide on a procedure? With a cond statement, naturally. Something like.

(cond (or (first?) (last?) (f x))
      (even? (* 4 (f x)))
      (odd? (* 2 (f x))))

Next I need to figure out how to recognize the first, last, even, and odd terms. The parameters at my disposal are bounds a, b, total number of intervals n, and the specific function argument x. Looking at my rewritten equation above, I can see the coefficient of h has the pattern 0, 1, 2, 3, ... , n. I can see how to use this to write a test.

  • If the coefficient of h = 0, it is the first term.
  • If the coefficient = n, it is the last term.
  • If the coefficient is odd, multiply the function results by 4.
  • If the coefficient is even, multiply the function results by 2.

The parameter x is equal to a + nkh. I can't get the coefficient directly, but together with a and h, I have enough elements to solve for for the coefficent, nk.


So my test becomes

(define (simp-term x)
  (cond ((or (= x a) (= x b))(f x))
        ((even? (/ (- x a) h)) (* 2 (f x)))
        ((odd? (/ (- x a) h))(* 4 (f x)))))

I should note that in reality I'm doing a lot of trial-and-error coding at this point. The DrRacket stepper was a godsend, also. In hindsight, I can see why what I came up with is correct, but I don't really approach the problem like a math proof. I approach it more of a jigsaw puzzle, trying to put things together until it works.

Putting it all together I end up with:
(define (simpsons f a b n)
  (define h
    (/ (- b a) n))
  (define (simp-next x)
    (+ x h))
  (define (simp-term x)
    (cond ((or (= x a) (= x b))(f x))
          ((even? (/ (- x a) h)) (* 2 (f x)))
          ((odd? (/ (- x a) h))(* 4 (f x)))))
  (* (/ h 3)
     (sum simp-term a simp-next b)))

I spent a lot of time on this problem. It just about killed me. Working on this blog post was hugely helpful (I started it long before I was finished with the problem). My reasoning was not quite as methodical as this blog post appears. Some things that tripped me up:

  • I was confused about data vs. procedures. At one point I was trying to make simp-next and simp-term return procedures. Wrong. The names represent a procedure, but they evaluate to a number.
  • I was confused about what simp-next should return. It needs to return the argument to f(x), not the result of f(x).
  • It wasn't clear what data I would have available when evaluating simp-term. I was getting hung up that I was defining simp-term inside simpsons, but it would be evaluated in sum. Sum doesn't have access to a, b, and n. Turns out I didn't need to worry where it would be evaluated.

I hope that reading this post help you as much as writing it helped me.