src/Integration/Trapezoid.hs

   1 module Integration.Trapezoid
   2 where
   3
   4
   5 -- | Partition the interval [@a@, @b@] into @n@ subintervals, which we
   6 --   then return as a list of pairs.
   7 partition :: (RealFrac a, Integral b)
   8           => b -- ^ The number of subintervals to use, @n@
   9           -> a -- ^ The \"left\" endpoint of the interval, @a@
  10           -> a -- ^ The \"right\" endpoint of the interval, @b@
  11           -> [(a,a)]
  12  -- Somebody asked for zero subintervals? Ok.
  13 partition 0 _ _ = []
  14 partition n a b
  15   | n < 0 = error "partition: asked for a negative number of subintervals"
  16   | otherwise =
  17     [ (xi, xj) | k <- [0..n-1],
  18                  let k' = fromIntegral k,
  19                  let xi = a + k'*h,
  20                  let xj = a + (k'+1)*h ]
  21     where
  22       h = fromRational $ (toRational (b-a))/(toRational n)
  23
  24
  25 -- | Use the trapezoid rule to numerically integrate @f@ over the
  26 --   interval [@a@, @b@].
  27 --
  28 --   Examples:
  29 --
  30 --   >>> let f x = x
  31 --   >>> trapezoid_1 f (-1) 1
  32 --   0.0
  33 --
  34 --   >>> let f x = x^3
  35 --   >>> trapezoid_1 f (-1) 1
  36 --   0.0
  37 --
  38 --   >>> let f x = 1
  39 --   >>> trapezoid_1 f (-1) 1
  40 --   2.0
  41 --
  42 --   >>> let f x = x^2
  43 --   >>> trapezoid_1 f (-1) 1
  44 --   2.0
  45 --
  46 trapezoid_1 :: (RealFrac a, Fractional b, Num b)
  47             => (a -> b) -- ^ The function @f@
  48             -> a       -- ^ The \"left\" endpoint, @a@
  49             -> a       -- ^ The \"right\" endpoint, @b@
  50             -> b
  51 trapezoid_1 f a b =
  52   (((f a) + (f b)) / 2) * (fromRational $ toRational (b - a))
  53
  54
  55 -- | Use the composite trapezoid tule to numerically integrate @f@
  56 --   over @n@ subintervals of [@a@, @b@].
  57 --
  58 --   Examples:
  59 --
  60 --   >>> let f x = x^2
  61 --   >>> let area = trapezoid 1000 f (-1) 1
  62 --   abs (area - (2/3)) < 0.00001
  63 --   True
  64 --
  65 --   >>> let area = trapezoid 1000 sin (-1) 1
  66 --   >>> abs (area - 2) < 0.00001
  67 --   True
  68 --
  69 trapezoid :: (RealFrac a, Fractional b, Num b, Integral c)
  70           => c -- ^ The number of subintervals to use, @n@
  71           -> (a -> b) -- ^ The function @f@
  72           -> a       -- ^ The \"left\" endpoint, @a@
  73           -> a       -- ^ The \"right\" endpoint, @b@
  74           -> b
  75 trapezoid n f a b =
  76   sum $ map trapezoid_pairs pieces
  77   where
  78     pieces = partition n a b
  79     -- Convert the trapezoid_1 function into one that takes pairs
  80     -- (a,b) instead of individual arguments 'a' and 'b'.
  81     trapezoid_pairs = uncurry (trapezoid_1 f)