src/Integration/Trapezoid.hs

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