Posts Tagged ‘integration’

Programming Praxis – Numerical Integration

February 9, 2010

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