In today’s Programming Praxis exercise we have to do some numerical integration. Let’s get started.

Since the first three algorithms all work roughly the same, and as a programmer I hate repeating myself, we’re going to make a generic integration function to abstract out the common stuff.

int combine f a b n = sum $ map g [0..n - 1] where w = (b - a) / n lo i = a + w * i g i = w * combine (f $ lo i) (f $ lo i + w/2) (f $ lo i + w)

Now we can just define the three simple integration methods in terms of what we want to do with the low, mid and high points.

intRect = int (\_ m _ -> m) intTrap = int (\l _ h -> (l + h) / 2) intSimp = int (\l m h -> (l + 4 * m + h) / 6)

The intAdapt function is pretty much the same as the Scheme version.

intAdapt m f a b epsilon = if abs (g10 - g5) < e then g10 else intAdapt m f a mid (Just e) + intAdapt m f mid b (Just e) where g5 = m f a b 5 g10 = m f a b 10 mid = (a + b) / 2 e = maybe 1e-7 id epsilon

Using adaptive integration for prime counting:

approxPi n = round $ intAdapt intSimp (recip . log) 2 n Nothing

And finally a test to show that everything’s working properly.

main = do print $ intRect cube 0 1 10000 print $ intTrap cube 0 1 10000 print $ intSimp cube 0 1 10000 print $ approxPi 1e21 where cube x = x * x * x