src/Integration/Simpson.hs

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