Add the Integration.Trapezoid module; import it in the .ghci file.

diff --git a/.ghci b/.ghci
index 82c4f29164c9880f5f83ecc6de63ce3b988e1da0..c87711e61eb2ae70b9a19a2d3c522150166f52b4 100644 (file)
--- a/.ghci
+++ b/.ghci
@@ -2,13 +2,14 @@
 :set -isrc
 
 -- Load everything.
-:l src/Roots/Simple.hs src/Aliases.hs src/TwoTuple.hs
+:l src/Roots/Simple.hs src/Aliases.hs src/TwoTuple.hs src/Integration/Trapezoid.hs
 
 -- Just for convenience.
 import Data.Number.BigFloat
 
 import Aliases
 import TwoTuple
+import Integration.Trapezoid
 
 -- Use a calmer prompt.
 :set prompt "numerical-analysis> "

diff --git a/src/Integration/Trapezoid.hs b/src/Integration/Trapezoid.hs
new file mode 100644 (file)
index 0000000..9598001
--- /dev/null
+++ b/src/Integration/Trapezoid.hs
@@ -0,0 +1,81 @@
+module Integration.Trapezoid
+where
+
+
+-- | Partition the interval [@a@, @b@] into @n@ subintervals, which we
+--   then return as a list of pairs.
+partition :: (RealFrac a, Integral b)
+          => b -- ^ The number of subintervals to use, @n@
+          -> a -- ^ The \"left\" endpoint of the interval, @a@
+          -> a -- ^ The \"right\" endpoint of the interval, @b@
+          -> [(a,a)]
+ -- Somebody asked for zero subintervals? Ok.
+partition 0 _ _ = []
+partition n a b
+  | n < 0 = error "partition: asked for a negative number of subintervals"
+  | otherwise =
+    [ (xi, xj) | k <- [0..n-1],
+                 let k' = fromIntegral k,
+                 let xi = a + k'*h,
+                 let xj = a + (k'+1)*h ]
+    where
+      h = fromRational $ (toRational (b-a))/(toRational n)
+
+
+-- | Use the trapezoid rule to numerically integrate @f@ over the
+--   interval [@a@, @b@].
+--
+--   Examples:
+--
+--   >>> let f x = x
+--   >>> trapezoid_1 f (-1) 1
+--   0.0
+--
+--   >>> let f x = x^3
+--   >>> trapezoid_1 f (-1) 1
+--   0.0
+--
+--   >>> let f x = 1
+--   >>> trapezoid_1 f (-1) 1
+--   2.0
+--
+--   >>> let f x = x^2
+--   >>> trapezoid_1 f (-1) 1
+--   2.0
+--
+trapezoid_1 :: (RealFrac a, Fractional b, Num b)
+            => (a -> b) -- ^ The function @f@
+            -> a       -- ^ The \"left\" endpoint, @a@
+            -> a       -- ^ The \"right\" endpoint, @b@
+            -> b
+trapezoid_1 f a b =
+  (((f a) + (f b)) / 2) * (fromRational $ toRational (b - a))
+
+
+-- | Use the composite trapezoid tule to numerically integrate @f@
+--   over @n@ subintervals of [@a@, @b@].
+--
+--   Examples:
+--
+--   >>> let f x = x^2
+--   >>> let area = trapezoid 1000 f (-1) 1
+--   abs (area - (2/3)) < 0.00001
+--   True
+--
+--   >>> let area = trapezoid 1000 sin (-1) 1
+--   >>> abs (area - 2) < 0.00001
+--   True
+--
+trapezoid :: (RealFrac a, Fractional b, Num b, Integral c)
+          => c -- ^ The number of subintervals to use, @n@
+          -> (a -> b) -- ^ The function @f@
+          -> a       -- ^ The \"left\" endpoint, @a@
+          -> a       -- ^ The \"right\" endpoint, @b@
+          -> b
+trapezoid n f a b =
+  sum $ map trapezoid_pairs pieces
+  where
+    pieces = partition n a b
+    -- Convert the trapezoid_1 function into one that takes pairs
+    -- (a,b) instead of individual arguments 'a' and 'b'.
+    trapezoid_pairs = uncurry (trapezoid_1 f)