Numerical Integration (With Precision)

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:

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:

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

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

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:

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:

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:

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:

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:

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

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?

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:

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:

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

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

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

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:

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

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

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

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