]> gitweb.michael.orlitzky.com - numerical-analysis.git/blob - src/Integration/Trapezoid.hs
Add the Integration.Trapezoid module; import it in the .ghci file.
[numerical-analysis.git] / src / Integration / Trapezoid.hs
1 module Integration.Trapezoid
2 where
3
4
5 -- | Partition the interval [@a@, @b@] into @n@ subintervals, which we
6 -- then return as a list of pairs.
7 partition :: (RealFrac a, Integral b)
8 => b -- ^ The number of subintervals to use, @n@
9 -> a -- ^ The \"left\" endpoint of the interval, @a@
10 -> a -- ^ The \"right\" endpoint of the interval, @b@
11 -> [(a,a)]
12 -- Somebody asked for zero subintervals? Ok.
13 partition 0 _ _ = []
14 partition n a b
15 | n < 0 = error "partition: asked for a negative number of subintervals"
16 | otherwise =
17 [ (xi, xj) | k <- [0..n-1],
18 let k' = fromIntegral k,
19 let xi = a + k'*h,
20 let xj = a + (k'+1)*h ]
21 where
22 h = fromRational $ (toRational (b-a))/(toRational n)
23
24
25 -- | Use the trapezoid rule to numerically integrate @f@ over the
26 -- interval [@a@, @b@].
27 --
28 -- Examples:
29 --
30 -- >>> let f x = x
31 -- >>> trapezoid_1 f (-1) 1
32 -- 0.0
33 --
34 -- >>> let f x = x^3
35 -- >>> trapezoid_1 f (-1) 1
36 -- 0.0
37 --
38 -- >>> let f x = 1
39 -- >>> trapezoid_1 f (-1) 1
40 -- 2.0
41 --
42 -- >>> let f x = x^2
43 -- >>> trapezoid_1 f (-1) 1
44 -- 2.0
45 --
46 trapezoid_1 :: (RealFrac a, Fractional b, Num b)
47 => (a -> b) -- ^ The function @f@
48 -> a -- ^ The \"left\" endpoint, @a@
49 -> a -- ^ The \"right\" endpoint, @b@
50 -> b
51 trapezoid_1 f a b =
52 (((f a) + (f b)) / 2) * (fromRational $ toRational (b - a))
53
54
55 -- | Use the composite trapezoid tule to numerically integrate @f@
56 -- over @n@ subintervals of [@a@, @b@].
57 --
58 -- Examples:
59 --
60 -- >>> let f x = x^2
61 -- >>> let area = trapezoid 1000 f (-1) 1
62 -- abs (area - (2/3)) < 0.00001
63 -- True
64 --
65 -- >>> let area = trapezoid 1000 sin (-1) 1
66 -- >>> abs (area - 2) < 0.00001
67 -- True
68 --
69 trapezoid :: (RealFrac a, Fractional b, Num b, Integral c)
70 => c -- ^ The number of subintervals to use, @n@
71 -> (a -> b) -- ^ The function @f@
72 -> a -- ^ The \"left\" endpoint, @a@
73 -> a -- ^ The \"right\" endpoint, @b@
74 -> b
75 trapezoid n f a b =
76 sum $ map trapezoid_pairs pieces
77 where
78 pieces = partition n a b
79 -- Convert the trapezoid_1 function into one that takes pairs
80 -- (a,b) instead of individual arguments 'a' and 'b'.
81 trapezoid_pairs = uncurry (trapezoid_1 f)