1 module Integration.Simpson
4 import Misc (partition)
7 -- | Use the Simpson's rule to numerically integrate @f@ over the
8 -- interval [@a@, @b@].
13 -- >>> simpson_1 f (-1) 1
17 -- >>> simpson_1 f (-1) 1
21 -- >>> let area = simpson_1 f (-1) 1
22 -- >>> abs (area - (2/3)) < 1/10^12
26 -- >>> simpson_1 f (-1) 1
30 -- >>> simpson_1 f 0 1
33 simpson_1 :: (RealFrac a, Fractional b, Num b)
34 => (a -> b) -- ^ The function @f@
35 -> a -- ^ The \"left\" endpoint, @a@
36 -> a -- ^ The \"right\" endpoint, @b@
39 coefficient * ((f a) + 4*(f midpoint) + (f b))
41 coefficient = (realToFrac (b - a)) / 6
42 midpoint = (a + b) / 2
45 -- | Use the composite Simpson's rule to numerically integrate @f@
46 -- over @n@ subintervals of [@a@, @b@].
51 -- >>> let area = simpson 10 f (-1) 1
52 -- >>> abs (area - (2/5)) < 0.0001
55 -- Note that the convergence here is much faster than the Trapezoid
58 -- >>> let area = simpson 10 sin 0 pi
59 -- >>> abs (area - 2) < 0.00001
62 simpson :: (RealFrac a, Fractional b, Num b, Integral c)
63 => c -- ^ The number of subintervals to use, @n@
64 -> (a -> b) -- ^ The function @f@
65 -> a -- ^ The \"left\" endpoint, @a@
66 -> a -- ^ The \"right\" endpoint, @b@
69 sum $ map simpson_pairs pieces
71 pieces = partition n a b
72 -- Convert the simpson_1 function into one that takes pairs
73 -- (a,b) instead of individual arguments 'a' and 'b'.
74 simpson_pairs = uncurry (simpson_1 f)