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 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 trapezoid rule to numerically integrate @f@ over the
  17 --   interval [@a@, @b@].
  18 --
  19 --   Examples:
  20 --
  21 --   >>> let f x = x
  22 --   >>> trapezoid_1 f (-1) 1
  23 --   0.0
  24 --
  25 --   >>> let f x = x^3
  26 --   >>> trapezoid_1 f (-1) 1
  27 --   0.0
  28 --
  29 --   >>> let f x = 1
  30 --   >>> trapezoid_1 f (-1) 1
  31 --   2.0
  32 --
  33 --   >>> let f x = x^2
  34 --   >>> trapezoid_1 f (-1) 1
  35 --   2.0
  36 --
  37 trapezoid_1 :: (Field.C a, ToRational.C a, Field.C b)
  38             => (a -> b) -- ^ The function @f@
  39             -> a       -- ^ The \"left\" endpoint, @a@
  40             -> a       -- ^ The \"right\" endpoint, @b@
  41             -> b
  42 trapezoid_1 f a b =
  43   (((f a) + (f b)) / 2) * (fromRational' $ toRational (b - a))
  44
  45
  46 -- | Use the composite trapezoid rule to numerically integrate @f@
  47 --   over @n@ subintervals of [@a@, @b@].
  48 --
  49 --   Examples:
  50 --
  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 --   >>> let area = trapezoid 1000 sin 0 pi
  57 --   >>> abs (area - 2) < 0.0001
  58 --   True
  59 --
  60 trapezoid :: (RealField.C a,
  61               ToRational.C a,
  62               RealField.C b,
  63               ToInteger.C c,
  64               Enum c)
  65           => c -- ^ The number of subintervals to use, @n@
  66           -> (a -> b) -- ^ The function @f@
  67           -> a       -- ^ The \"left\" endpoint, @a@
  68           -> a       -- ^ The \"right\" endpoint, @b@
  69           -> b
  70 trapezoid n f a b =
  71   sum $ map trapezoid_pairs pieces
  72   where
  73     pieces = partition n a b
  74     -- Convert the trapezoid_1 function into one that takes pairs
  75     -- (a,b) instead of individual arguments 'a' and 'b'.
  76     trapezoid_pairs = uncurry (trapezoid_1 f)