Implementation of Simpson method of numerical integration.

Integration.Simpson.hs

-- file Integration.Simpson.hs xdCrafts 06.11.2012
module
	Integration.Simpson (
		integrate,
		integrateFromN,
		integrateWithAccuracy
	)
where

-- function that splits tuple by odd and even components
splitOddsAndEvens :: [a] -> ([a], [a])
splitOddsAndEvens [] = ([], [])
splitOddsAndEvens [x] = ([x], [])
splitOddsAndEvens (x:y:xs) = (x:xp, y:yp) where (xp, yp) = splitOddsAndEvens xs

-- http://en.wikipedia.org/wiki/Simpson%27s_rule
-- f - function
-- a - begining of interval
-- b - end of inteval
-- n - number of intervals of size h
integrate :: (Fractional f) => (f -> f) -> f -> f -> Int -> f
integrate f a b n =
	h / 3 * (fx 0 + 4 * (sum $ map fx odds) + 2 * (sum $ map fx evens) +  fx (fromIntegral (2 * n + 2)))
	where
		h = (b - a) / fromIntegral (2 * n)
		x i = a + h * i
		fx = f . x
		odds = fst oddsAndEvens
		evens = snd  oddsAndEvens
		oddsAndEvens = splitOddsAndEvens (map toDouble [1..2 * n - 1])
		toDouble int = fromIntegral int

-- http://en.wikipedia.org/wiki/Simpson%27s_rule
-- f - function
-- a - begining of interval
-- b - end of inteval
-- n - number of intervals of size h
-- e - accuracy
integrateFromN :: (Fractional f, Ord f) => (f -> f) -> f -> f -> Int -> f -> f
integrateFromN f a b n e
	| accuracy < e = solution2
	| otherwise = solution1
	where
		accuracy = abs (solution2 - solution1)
		solution1 = integrate f a b n
		solution2 = integrate f a b (2 * n)

-- http://en.wikipedia.org/wiki/Simpson%27s_rule
-- f - function
-- a - begining of interval
-- b - end of inteval
-- e - acuracy
integrateWithAccuracy :: (Fractional f, Ord f) => (f -> f) -> f -> f -> f -> f
integrateWithAccuracy f a b e = integrateFromN f a b 10 e