From: Michael Orlitzky Date: Fri, 23 Nov 2012 18:40:23 +0000 (-0500) Subject: Add Integration.Simpson implementing Simpson's rule for numerical integration. X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=4b5ff37009f1070059697a5f396987bb109a882d;p=numerical-analysis.git Add Integration.Simpson implementing Simpson's rule for numerical integration. --- diff --git a/.ghci b/.ghci index c87711e..a04730c 100644 --- 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 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)