Programming Praxis - Numerical Integration

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

About these ads

Tags: bonsai, code, Haskell, integration, kata, numerical, praxis, programming