X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FIntegration%2FSimpson.hs;fp=src%2FIntegration%2FSimpson.hs;h=b7000224764e92095398837dedf3d56eaf4ee96c;hb=4b5ff37009f1070059697a5f396987bb109a882d;hp=0000000000000000000000000000000000000000;hpb=92d61dc6ab8dc3e8c9e5b4c869134f1e8079013f;p=numerical-analysis.git

diff --git a/src/Integration/Simpson.hs b/src/Integration/Simpson.hs
new file mode 100644
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)