src/Integration/Simpson.hs

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