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