## 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