Numerical Integration (With Precision)

Update: Changed from nested defn-  inside defn  to instead define local functions in let -bindings.

In a previous blog entry, I explained the higher-order function sum and how to use the Substitution Model to follow the execution of a function. In this entry, I will use the sum function to perform numerical integration, and in the process run into some pitfalls of Java’s BigDecimal. I will show how the language Clojure provides an elegant solution to the “exact quotient cannot be represented” problem of BigDecimal.

We will once again look at Section 1.3 of Structure and Interpretation of Computer Programs, Formulating Abstractions with Higher-Order Procedures. To recap, the sum function is an abstraction of a sum of some sequence. It takes a term function for calculating the current term, a start value a, a next function for calculating the next value of a, and an end value b:

(defn sum [term a next b]
  (if (> a b)
    0
    (+ (term a)
       (sum term (next a) next b))))

After introducing the sum function, the authors describe a simple sum for calculating the integral of a function:

Looking at the first term, we see that (a+dx/2) is the start value, ie what we will pass as the a parameter to sum. The term function for calculating the value of the current term is obviously f.

We can also see that each term has one dx added to it. The first term has no dx, the second term has one dx, the third term has 2dx, etc. A simple function adding dx to its argument will be perfect as the next function. We’ll define it as a private function using defn-, to not clutter the namespace with internal functions. However, it needs to know the dx, so we’ll define it inside the integral function, like this:

(defn integral [f a b dx]
  (let [add-dx (fn [x] (+ x dx))]
    ,,,))

Finally, we need to multiply the result of calling sum by dx. Here is the implementation of the integral function:

(defn integral [f a b dx]
  (let [add-dx (fn [x] (+ x dx))]
    (* (sum f (+ a (/ dx 2)) add-dx b)
       dx)))

We’ll need a function that we can pass as f, so we can do some interesting integration. Let’s use cube:

(defn cube [x] (* x x x))

The integral of x^3 is x^4/4. Integrating from 0 to 1 results in (1^4)/4 = 1/4 = 0.25. That’s what we will compare with.

Let’s test our integral function using a dx of, say 0.01:

user> (integral cube 0 1 0.01)
0.24998750000000042

Not that bad. However, in Exercise 1.29, the authors suggest using Simpson’s Rule, which is a more accurate method of numerical integration.

Simpson’s Rule

Simpson’s Rule for numerical integration goes like this:

OK, so how do we squeeze this into the sum function? First we note that h is used in more than one place, so we’ll begin by creating a local binding h. We’ll also multiply the result of sum with h/3:

(defn simpson-sum [f a b n]
  (let [h (/ (- b a) n)]
	,,,
    (* (sum ,,, ,,, ,,, ,,,)
       (/ h 3))))

Next, we see that the factor of each term is alternating between 4 and 2, except the first and the last one, which are 1. This can be implemented as a conditional:

(cond (zero? k) y
      (= k n)   y
      (odd? k)  (* 4 y)
      (even? k) (* 2 y))

The k will be provided as a parameter, as we’ll see in a second, but we’re missing the y. So let’s add the expression of y and package it all up in a function inside simpson-sum:

(defn simpson-sum [f a b n]
  (let [h (/ (- b a) n)
        simpson-term (fn [k]
                       (let [y (f (+ a (* k h)))]
                         (cond (zero? k) y
                               (= k n) y
                               (odd? k) (* 4 y)
                               (even? k) (* 2 y))))]
    (* (sum ,,, ,,, ,,, ,,,)
       (/ h 3))))

OK, so now we have something that, given a k, provides the corresponding term in the sum, with the correct factor. The sequence for the sum is from 0 to n, so all we need to do now is to call sum with simpson-term as the term function, zero as the start value of the sequence, inc as the next function, and n as the end of the sequence:

(defn simpson-sum [f a b n]
  (let [h (/ (- b a) n)
        simpson-term (fn [k]
                       (let [y (f (+ a (* k h)))]
                         (cond (zero? k) y
                               (= k n) y
                               (odd? k) (* 4 y)
                               (even? k) (* 2 y))))]
    (* (sum simpson-term 0 inc n)
       (/ h 3))))

The test run on cube using n=100 gives some impressive results:

user> (simpson-sum cube 0 1 100)
1/4
user> (simpson-sum cube 0 2 100)
4
user> (simpson-sum cube 3 4 100)
175/4

Not bad at all. Clojure’s use of Ratio enables us to get accurate precision. But what happens if I select a floating point interval?

user> (simpson-sum cube 0 0.5 100)
0.015624999999999995

Pretty good that too.

Using BigDecimal

Now, what if I want to use BigDecimal in order to get really high precision in my floating points? In Clojure, the letter M on a number indicates that it’s a BigDecimal literal:

user> (simpson-sum cube 0 0.5M 100)
Non-terminating decimal expansion; no exact representable decimal result.
  [Thrown class java.lang.ArithmeticException]
Backtrace:
  0: java.math.BigDecimal.divide(BigDecimal.java:1603)
  1: clojure.lang.Numbers$BigDecimalOps.divide(Numbers.java:1067)
  2: clojure.lang.Numbers.divide(Numbers.java:139)
  3: user$simpson_sum.invoke(NO_SOURCE_FILE:1)

What happened here? Actually, this is something that can potentially happen whenever using BigDecimal divison, whether it’s from Java, Clojure or something else. The JavaDoc for BigDecimal.divide states:

divide(BigDecimal divisor)
Returns a BigDecimal whose value is (this / divisor), and whose preferred scale is (this.scale() – divisor.scale()); if the exact quotient cannot be represented (because it has a non-terminating decimal expansion) an ArithmeticException is thrown.

Could it be that division of h by 3 that messes things up? Division by 3 is always troublesome when using floating point numbers. Anyway, in order to solve this, we need to use a version of BigDecimal.divide that takes a scale (we’ll use 10) and a rounding mode. Only line 10 needs to be changed:

(defn simpson-sum [f a b n]
  (let [h (/ (- b a) n)
        simpson-term (fn [k]
                       (let [y (f (+ a (* k h)))]
                         (cond (zero? k) y
                               (= k n) y
                               (odd? k) (* 4 y)
                               (even? k) (* 2 y))))]
    (* (sum simpson-term 0 inc n)
       (.divide (bigdec h) 3M 10 java.math.RoundingMode/HALF_UP))))

Ugly as hell, but what can you do? (Plenty, actually, as we’ll see soon) Let’s test it:

user> (simpson-sum cube 0 0.5M 100)
0.0156250031250000000000M

OK, so it works. Are we happy now? Not really. Look what happens if we use an n that is a multiple of 3:

user> (simpson-sum cube 0 0.5M 6)
Non-terminating decimal expansion; no exact representable decimal result.
  [Thrown class java.lang.ArithmeticException]
Backtrace:
  0: java.math.BigDecimal.divide(BigDecimal.java:1603)
  1: clojure.lang.Numbers$BigDecimalOps.divide(Numbers.java:1067)
  2: clojure.lang.Numbers.divide(Numbers.java:139)
  3: user$simpson_sum.invoke(NO_SOURCE_FILE:1)

This time it’s the other division that is the problem. We’ll change line 2 as well:

(defn simpson-sum [f a b n]
  (let [h (.divide (bigdec (- b a)) (bigdec n) 10 java.math.RoundingMode/HALF_UP)
        simpson-term (fn [k]
                       (let [y (f (+ a (* k h)))]
                         (cond (zero? k) y
                               (= k n) y
                               (odd? k) (* 4 y)
                               (even? k) (* 2 y))))]
    (* (sum simpson-term 0 inc n)
       (.divide (bigdec h) 3M 10 java.math.RoundingMode/HALF_UP))))
user> (simpson-sum cube 0 0.5M 6)
0.0156249999937499999925000000049999999992M

So, we’re good now? No, Clojure has actually a much better way of handling this than resorting to the above. Let’s look at the documentation for the function with-precision:

user> (doc with-precision)
-------------------------
clojure.core/with-precision
([precision & exprs])
Macro
  Sets the precision and rounding mode to be used for BigDecimal operations.

  Usage: (with-precision 10 (/ 1M 3))
  or:    (with-precision 10 :rounding HALF_DOWN (/ 1M 3))

  The rounding mode is one of CEILING, FLOOR, HALF_UP, HALF_DOWN,
  HALF_EVEN, UP, DOWN and UNNECESSARY; it defaults to HALF_UP.

OK, so we simply call with-precision with the expression we want precision on. Pretty elegant. Here is the final implementation:

(defn simpson-sum [f a b n]
  (with-precision 10
    (let [h (/ (- b a) n)
          simpson-term (fn [k]
                         (let [y (f (+ a (* k h)))]
                           (cond (zero? k) y
                                 (= k n) y
                                 (odd? k) (* 4 y)
                                 (even? k) (* 2 y))))]
      (* (sum simpson-term 0 inc n)
         (/ h 3)))))

We’re now getting even better results than when manually providing a scale to BigDecimal.divide. Here is a run with n=6:

user> (simpson-sum cube 0 0.5M 6)
0.01562500000M

Compare that with the integral function that we implemented first. That would have given us this result:

user> (integral cube 0 0.5 0.01)
0.015621875000000023

The exact integral of x^3 between 0 and 0.5 is (0.5^4)/4 = 0.015625, so it’s clear that Simpson’s Rule is very good for numerical integration. The code becomes really elegant even when using BigDecimal arithmetic, thanks to Clojure’s support for BigDecimal precision.

References

Leave a Reply

Close Menu