1 {-# LANGUAGE RebindableSyntax #-}
3 module Integration.Trapezoid
6 import Misc (partition)
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
14 -- | Use the trapezoid rule to numerically integrate @f@ over the
15 -- interval [@a@, @b@].
20 -- >>> trapezoid_1 f (-1) 1
24 -- >>> trapezoid_1 f (-1) 1
28 -- >>> trapezoid_1 f (-1) 1
32 -- >>> trapezoid_1 f (-1) 1
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@
41 (((f a) + (f b)) / 2) * coerced_interval_length
43 coerced_interval_length = fromRational' $ toRational (b - a)
45 -- | Use the composite trapezoid rule to numerically integrate @f@
46 -- over @n@ subintervals of [@a@, @b@].
50 -- >>> import Algebra.Absolute (abs)
52 -- >>> let area = trapezoid 1000 f (-1) 1
53 -- >>> abs (area - (2/3)) < 0.00001
56 -- >>> import Algebra.Absolute (abs)
57 -- >>> let area = trapezoid 1000 sin 0 pi
58 -- >>> abs (area - 2) < 0.0001
61 trapezoid :: (RealField.C a,
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@
72 sum $ map trapezoid_pairs pieces
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)