From: Michael Orlitzky Date: Wed, 21 Nov 2012 03:02:27 +0000 (-0500) Subject: Add the Integration.Trapezoid module; import it in the .ghci file. X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=57a982c464095201f7635977aac937ad6a08c0b0;p=numerical-analysis.git Add the Integration.Trapezoid module; import it in the .ghci file. --- diff --git a/.ghci b/.ghci index 82c4f29..c87711e 100644 --- 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 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)