src/Integration/Simpson.hs

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