{-# OPTIONS_GHC -Wall #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE Safe #-} {- | Module : Physics.Learn.CompositeQuadrature Copyright : (c) Scott N. Walck 2012-2018 License : BSD3 (see LICENSE) Maintainer : Scott N. Walck Stability : experimental Composite Trapezoid Rule and Composite Simpson's Rule -} module Physics.Learn.CompositeQuadrature ( compositeTrapezoid , compositeSimpson ) where import Data.VectorSpace ( VectorSpace , Scalar , (^+^) , (*^) , zeroV ) -- | Composite Trapezoid Rule compositeTrapezoid :: (VectorSpace v, Fractional (Scalar v)) => Int -- ^ number of intervals (one less than the number of function evaluations) -> Scalar v -- ^ lower limit -> Scalar v -- ^ upper limit -> (Scalar v -> v) -- ^ function to be integrated -> v -- ^ definite integral compositeTrapezoid n a b f = let dt = (b - a) / fromIntegral n ts = [a + fromIntegral m * dt | m <- [0..n]] pairs = [(t,f t) | t <- ts] combine [] = error "compositeSimpson: odd number of half-intervals" -- this should never happen combine [_] = zeroV combine ((t0,f0):(t1,f1):ps) = ((t1 - t0) / 2) *^ (f0 ^+^ f1) ^+^ combine ((t1,f1):ps) in combine pairs -- | Composite Simpson's Rule compositeSimpson :: (VectorSpace v, Fractional (Scalar v)) => Int -- ^ number of half-intervals (one less than the number of function evaluations) -> Scalar v -- ^ lower limit -> Scalar v -- ^ upper limit -> (Scalar v -> v) -- ^ function to be integrated -> v -- ^ definite integral compositeSimpson n a b f = let nEven = 2 * div n 2 dt = (b - a) / fromIntegral nEven ts = [a + fromIntegral m * dt | m <- [0..nEven]] pairs = [(t,f t) | t <- ts] combine [] = error "compositeSimpson: odd number of half-intervals" -- this should never happen combine [_] = zeroV combine (_:_:[]) = error "compositeSimpson: odd number of half-intervals" -- this should never happen combine ((t0,f0):(_,f1):(t2,f2):ps) = ((t2 - t0) / 6) *^ (f0 ^+^ 4 *^ f1 ^+^ f2) ^+^ combine ((t2,f2):ps) in combine pairs