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 Algebra.Absolute (abs)
  10 import qualified Algebra.Field as Field
  11 import qualified Algebra.RealField as RealField
  12 import qualified Algebra.RealRing as RealRing
  13 import qualified Algebra.ToInteger as ToInteger
  14 import qualified Algebra.ToRational as ToRational
  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 --   >>> let f x = x^2
  30 --   >>> let area = simpson_1 f (-1) 1
  31 --   >>> abs (area - (2/3)) < 1/10^12
  32 --   True
  33 --
  34 --   >>> let f x = x^3
  35 --   >>> simpson_1 f (-1) 1
  36 --   0.0
  37 --
  38 --   >>> let f x = x^3
  39 --   >>> simpson_1 f 0 1
  40 --   0.25
  41 --
  42 simpson_1 :: (RealField.C a, ToRational.C a, RealField.C b)
  43             => (a -> b) -- ^ The function @f@
  44             -> a       -- ^ The \"left\" endpoint, @a@
  45             -> a       -- ^ The \"right\" endpoint, @b@
  46             -> b
  47 simpson_1 f a b =
  48   coefficient * ((f a) + 4*(f midpoint) + (f b))
  49   where
  50     coefficient = (fromRational' $ toRational (b - a)) / 6
  51     midpoint = (a + b) / 2
  52
  53
  54 -- | Use the composite Simpson's rule to numerically integrate @f@
  55 --   over @n@ subintervals of [@a@, @b@].
  56 --
  57 --   Examples:
  58 --
  59 --   >>> let f x = x^4
  60 --   >>> let area = simpson 10 f (-1) 1
  61 --   >>> abs (area - (2/5)) < 0.0001
  62 --   True
  63 --
  64 --   Note that the convergence here is much faster than the Trapezoid
  65 --   rule!
  66 --
  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)