Mattias Holmqvist – Writings about software development

Inspired by the latest renaissance of functional programming, I decided to take on some real-world examples of calculations using Clojure. This is an exercise from the book Structure and Interpretation of Computer Programs.

There is a method for numeric integration called Simpson’s rule which can be used to approximate the integral of a given function in some interval . I decided to try some Clojure on this.

The rule states that this approximation can be found by calculating:

where and .

First off, we can make simplify the problem a bit. We can break out the first and last element in the series (by calculating for the 0′th and n’th index) and also realize that the odd and even elements in the series have common factors. This yields:

If we break this down into separate problems, we realize that we first need to find the indices for the first and second series, respectively. Once we have found the indices, we can apply the function to the value of for each index.

Finding the indices turns out to be quite simple. We can generate a complete sequence (1-n) and then filter out all interesting values and return that sequence. The range function generates a sequence and the filter function filters the sequence and returns all values that satisfies a supplied predicate function. Luckily for us, Clojure provides the predicate functions odd? and even?, which answers to whether or not numbers are odd or even. As an example from the REPL:

(filter odd? (range 1 10))
(1 3 5 7 9)

We can use this method in the function sub-series which will be responsible for generating our sums:

(defn sub-series [n index-filter]
  (map #(f (+ a (* h %))) (filter index-filter (range 1 n))))

The sub-series function takes as arguments the function for which we’re approximating the integral, the start of the integration interval , the granularity of the integral , the value as defined in the problem statement and the predicate function which is used to determine which indices to include in the series.

This function will apply the function for each value in the series (from 1-n) that is satisfied by the given predicate function. Once the function has been applied to all those values, it will return the calculated values as a sequence.

(defn sub-series [f a h n index-filter]
  (map #(f (+ a (* h %))) (filter index-filter (range 1 n))))

I’m using a few Clojure API functions (map and filter) but I will not explain these further here. For more info about these functions, see the Clojure API.

Since the first series we want to generate only include even indexes, we pass the even? predicate function to our sub-series function. For the second series we pass the odd? predicate function. We also multiply all values in each sequence with the proper factor (2 and 4, respectively). We use reduce to add together all values in the generated sequences. As a final step we also add and and multiply the whole sum with :

(defn simpson [f a b n]
  (def h (/ (- b a ) n))
  (* (/ h 3) (+
    (reduce + (map #(* % 2) (sub-series f a h n even?)))
    (reduce + (map #(* % 4) (sub-series f a h n odd?)))
    (f a)
    (f b))))

Phew. Let’s try this out. First let’s try it on a simple function such as . We can use the identity function from the Clojure API for this:

(simpson identity 0.0 1.0 100)
0.5

So far so good…

Second, we’ll try a more sophisticated polynomial function, . We’ll define the polynomial:

(defn poly [x]
  (- (+ (* 7 x) (* 3 (cube x))) (* 8.5 (square x))))

And approximate the integral in :

(simpson poly 0.0 2.0 100)
3.3333333333333335

which is correct.

To clean the whole thing up I’ll letting the sub-series function be a local function. This yields the final version of our Simpson-approximation:

(defn simpson [f a b n]
  (let [h (/ (- b a ) n)]
  (letfn [
      (sub-series [index-filter]
        (map #(f (+ a (* h %))) (filter index-filter (range 1 n))))]
    (* (/ h 3) (+
    (reduce + (map #(* % 2) (sub-series even?)))
    (reduce + (map #(* % 4) (sub-series odd?)))
    (f a)
    (f b))))))

I think the end result of 10 LOC is quite neat. The rich API with filter, map, reduce etc is really fun to work with.

If you ever need to make an integral calculation, you can just paste this code and you have your area!

Simplified solution

As pointed out by Chris B, we don’t need to use map in this case, since we only multiply the sub-series values with constants (2 and 4, respectively) we can simply multiply the reduced sub-series with these values. If the sub-series values should be multiplied with a function that depends on the current value in the sub-series, map would be a good approach. This simplifies the solution a bit:

(defn simpson [f a b n]
  (let [h (/ (- b a ) n)]
  (letfn [
      (sub-series [index-filter]
        (map #(f (+ a (* h %))) (filter index-filter (range 1 n))))]
    (* (/ h 3) (+
    (* 2 (reduce + (sub-series even?)))
    (* 4 (reduce + (sub-series odd?)))
    (f a)
    (f b))))))