src/Integration/Trapezoid.hs

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