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