src/Integration/Trapezoid.hs

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