Add Integration.Simpson implementing Simpson's rule for numerical integration.

diff --git a/.ghci b/.ghci
index c87711e61eb2ae70b9a19a2d3c522150166f52b4..a04730cc12b26ffccd4bf1bba0dd54f5dea864df 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 src/Integration/Trapezoid.hs
+:l src/Roots/Simple.hs src/Aliases.hs src/TwoTuple.hs src/Integration/Simpson.hs src/Integration/Trapezoid.hs
 
 -- Just for convenience.
 import Data.Number.BigFloat
 
 import Aliases
 import TwoTuple
+import Integration.Simpson
 import Integration.Trapezoid
 
 -- Use a calmer prompt.

diff --git a/src/Integration/Simpson.hs b/src/Integration/Simpson.hs
new file mode 100644 (file)
index 0000000..b700022
--- /dev/null
+++ b/src/Integration/Simpson.hs
@@ -0,0 +1,74 @@
+module Integration.Simpson
+where
+
+import Misc (partition)
+
+
+-- | Use the Simpson's rule to numerically integrate @f@ over the
+--   interval [@a@, @b@].
+--
+--   Examples:
+--
+--   >>> let f x = 1
+--   >>> simpson_1 f (-1) 1
+--   2.0
+--
+--   >>> let f x = x
+--   >>> simpson_1 f (-1) 1
+--   0.0
+--
+--   >>> let f x = x^2
+--   >>> let area = simpson_1 f (-1) 1
+--   >>> abs (area - (2/3)) < 1/10^12
+--   True
+--
+--   >>> let f x = x^3
+--   >>> simpson_1 f (-1) 1
+--   0.0
+--
+--   >>> let f x = x^3
+--   >>> simpson_1 f 0 1
+--   0.25
+--
+simpson_1 :: (RealFrac a, Fractional b, Num b)
+            => (a -> b) -- ^ The function @f@
+            -> a       -- ^ The \"left\" endpoint, @a@
+            -> a       -- ^ The \"right\" endpoint, @b@
+            -> b
+simpson_1 f a b =
+  coefficient * ((f a) + 4*(f midpoint) + (f b))
+  where
+    coefficient = (fromRational $ toRational (b - a)) / 6
+    midpoint = (a + b) / 2
+
+
+-- | Use the composite Simpson's rule to numerically integrate @f@
+--   over @n@ subintervals of [@a@, @b@].
+--
+--   Examples:
+--
+--   >>> let f x = x^4
+--   >>> let area = simpson 10 f (-1) 1
+--   >>> abs (area - (2/5)) < 0.0001
+--   True
+--
+--   Note that the convergence here is much faster than the Trapezoid
+--   rule!
+--
+--   >>> let area = simpson 10 sin 0 pi
+--   >>> abs (area - 2) < 0.00001
+--   True
+--
+simpson :: (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
+simpson n f a b =
+  sum $ map simpson_pairs pieces
+  where
+    pieces = partition n a b
+    -- Convert the simpson_1 function into one that takes pairs
+    -- (a,b) instead of individual arguments 'a' and 'b'.
+    simpson_pairs = uncurry (simpson_1 f)