src/Integration/Simpson.hs

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