1 {-# LANGUAGE NoImplicitPrelude #-}
2 {-# LANGUAGE ScopedTypeVariables #-}
3 {-# LANGUAGE RebindableSyntax #-}
5 module Integration.Simpson (
10 import Misc ( partition )
12 import NumericPrelude hiding ( abs )
13 import qualified Algebra.RealField as RealField ( C )
14 import qualified Algebra.ToInteger as ToInteger ( C )
15 import qualified Algebra.ToRational as ToRational ( C )
17 -- | Use the Simpson's rule to numerically integrate @f@ over the
18 -- interval [@a@, @b@].
23 -- >>> simpson_1 f (-1) 1
27 -- >>> simpson_1 f (-1) 1
30 -- >>> import Algebra.Absolute (abs)
32 -- >>> let area = simpson_1 f (-1) 1
33 -- >>> abs (area - (2/3)) < 1/10^12
37 -- >>> simpson_1 f (-1) 1
41 -- >>> simpson_1 f 0 1
44 simpson_1 :: forall a b. (RealField.C a, ToRational.C a, RealField.C b)
45 => (a -> b) -- ^ The function @f@
46 -> a -- ^ The \"left\" endpoint, @a@
47 -> a -- ^ The \"right\" endpoint, @b@
50 coefficient * ((f x) + 4*(f midpoint) + (f y))
52 coefficient = fromRational' $ (toRational (y - x)) / 6 :: b
53 midpoint = (x + y) / 2
56 -- | Use the composite Simpson's rule to numerically integrate @f@
57 -- over @n@ subintervals of [@a@, @b@].
61 -- >>> import Algebra.Absolute (abs)
63 -- >>> let area = simpson 10 f (-1) 1
64 -- >>> abs (area - (2/5)) < 0.0001
67 -- Note that the convergence here is much faster than the Trapezoid
70 -- >>> import Algebra.Absolute (abs)
71 -- >>> let area = simpson 10 sin 0 pi
72 -- >>> abs (area - 2) < 0.00001
75 simpson :: (RealField.C a,
80 => c -- ^ The number of subintervals to use, @n@
81 -> (a -> b) -- ^ The function @f@
82 -> a -- ^ The \"left\" endpoint, @a@
83 -> a -- ^ The \"right\" endpoint, @b@
86 sum $ map simpson_pairs pieces
88 pieces = partition n a b
89 -- Convert the simpson_1 function into one that takes pairs
90 -- (a,b) instead of individual arguments 'a' and 'b'.
91 simpson_pairs = uncurry (simpson_1 f)