src/Integration/Trapezoid.hs

   1 module Integration.Trapezoid
   2 where
   3
   4 import Misc (partition)
   5
   6 -- | Use the trapezoid rule to numerically integrate @f@ over the
   7 --   interval [@a@, @b@].
   8 --
   9 --   Examples:
  10 --
  11 --   >>> let f x = x
  12 --   >>> trapezoid_1 f (-1) 1
  13 --   0.0
  14 --
  15 --   >>> let f x = x^3
  16 --   >>> trapezoid_1 f (-1) 1
  17 --   0.0
  18 --
  19 --   >>> let f x = 1
  20 --   >>> trapezoid_1 f (-1) 1
  21 --   2.0
  22 --
  23 --   >>> let f x = x^2
  24 --   >>> trapezoid_1 f (-1) 1
  25 --   2.0
  26 --
  27 trapezoid_1 :: (RealFrac a, Fractional b, Num b)
  28             => (a -> b) -- ^ The function @f@
  29             -> a       -- ^ The \"left\" endpoint, @a@
  30             -> a       -- ^ The \"right\" endpoint, @b@
  31             -> b
  32 trapezoid_1 f a b =
  33   (((f a) + (f b)) / 2) * (fromRational $ toRational (b - a))
  34
  35
  36 -- | Use the composite trapezoid rule to numerically integrate @f@
  37 --   over @n@ subintervals of [@a@, @b@].
  38 --
  39 --   Examples:
  40 --
  41 --   >>> let f x = x^2
  42 --   >>> let area = trapezoid 1000 f (-1) 1
  43 --   >>> abs (area - (2/3)) < 0.00001
  44 --   True
  45 --
  46 --   >>> let area = trapezoid 1000 sin 0 pi
  47 --   >>> abs (area - 2) < 0.0001
  48 --   True
  49 --
  50 trapezoid :: (RealFrac a, Fractional b, Num b, Integral c)
  51           => c -- ^ The number of subintervals to use, @n@
  52           -> (a -> b) -- ^ The function @f@
  53           -> a       -- ^ The \"left\" endpoint, @a@
  54           -> a       -- ^ The \"right\" endpoint, @b@
  55           -> b
  56 trapezoid n f a b =
  57   sum $ map trapezoid_pairs pieces
  58   where
  59     pieces = partition n a b
  60     -- Convert the trapezoid_1 function into one that takes pairs
  61     -- (a,b) instead of individual arguments 'a' and 'b'.
  62     trapezoid_pairs = uncurry (trapezoid_1 f)