]> gitweb.michael.orlitzky.com - numerical-analysis.git/blob - src/Integration/Simpson.hs
Fix compiler warnings and doctests.
[numerical-analysis.git] / src / Integration / Simpson.hs
1 {-# LANGUAGE RebindableSyntax #-}
2
3 module Integration.Simpson
4 where
5
6 import Misc (partition)
7
8 import NumericPrelude hiding (abs)
9 import qualified Algebra.RealField as RealField
10 import qualified Algebra.ToInteger as ToInteger
11 import qualified Algebra.ToRational as ToRational
12
13 -- | Use the Simpson's rule to numerically integrate @f@ over the
14 -- interval [@a@, @b@].
15 --
16 -- Examples:
17 --
18 -- >>> let f x = 1
19 -- >>> simpson_1 f (-1) 1
20 -- 2.0
21 --
22 -- >>> let f x = x
23 -- >>> simpson_1 f (-1) 1
24 -- 0.0
25 --
26 -- >>> import Algebra.Absolute (abs)
27 -- >>> let f x = x^2
28 -- >>> let area = simpson_1 f (-1) 1
29 -- >>> abs (area - (2/3)) < 1/10^12
30 -- True
31 --
32 -- >>> let f x = x^3
33 -- >>> simpson_1 f (-1) 1
34 -- 0.0
35 --
36 -- >>> let f x = x^3
37 -- >>> simpson_1 f 0 1
38 -- 0.25
39 --
40 simpson_1 :: (RealField.C a, ToRational.C a, RealField.C b)
41 => (a -> b) -- ^ The function @f@
42 -> a -- ^ The \"left\" endpoint, @a@
43 -> a -- ^ The \"right\" endpoint, @b@
44 -> b
45 simpson_1 f a b =
46 coefficient * ((f a) + 4*(f midpoint) + (f b))
47 where
48 coefficient = (fromRational' $ toRational (b - a)) / 6
49 midpoint = (a + b) / 2
50
51
52 -- | Use the composite Simpson's rule to numerically integrate @f@
53 -- over @n@ subintervals of [@a@, @b@].
54 --
55 -- Examples:
56 --
57 -- >>> import Algebra.Absolute (abs)
58 -- >>> let f x = x^4
59 -- >>> let area = simpson 10 f (-1) 1
60 -- >>> abs (area - (2/5)) < 0.0001
61 -- True
62 --
63 -- Note that the convergence here is much faster than the Trapezoid
64 -- rule!
65 --
66 -- >>> import Algebra.Absolute (abs)
67 -- >>> let area = simpson 10 sin 0 pi
68 -- >>> abs (area - 2) < 0.00001
69 -- True
70 --
71 simpson :: (RealField.C a,
72 ToRational.C a,
73 RealField.C b,
74 ToInteger.C c,
75 Enum c)
76 => c -- ^ The number of subintervals to use, @n@
77 -> (a -> b) -- ^ The function @f@
78 -> a -- ^ The \"left\" endpoint, @a@
79 -> a -- ^ The \"right\" endpoint, @b@
80 -> b
81 simpson n f a b =
82 sum $ map simpson_pairs pieces
83 where
84 pieces = partition n a b
85 -- Convert the simpson_1 function into one that takes pairs
86 -- (a,b) instead of individual arguments 'a' and 'b'.
87 simpson_pairs = uncurry (simpson_1 f)