{-# 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 <walck@lvc.edu>
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 :: forall v.
(VectorSpace v, Fractional (Scalar v)) =>
Int -> Scalar v -> Scalar v -> (Scalar v -> v) -> v
compositeTrapezoid Int
n Scalar v
a Scalar v
b Scalar v -> v
f
    = let dt :: Scalar v
dt = (Scalar v
b forall a. Num a => a -> a -> a
- Scalar v
a) forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
          ts :: [Scalar v]
ts = [Scalar v
a forall a. Num a => a -> a -> a
+ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
m forall a. Num a => a -> a -> a
* Scalar v
dt | Int
m <- [Int
0..Int
n]]
          pairs :: [(Scalar v, v)]
pairs = [(Scalar v
t,Scalar v -> v
f Scalar v
t) | Scalar v
t <- [Scalar v]
ts]
          combine :: [(Scalar b, b)] -> b
combine [] = forall a. HasCallStack => [Char] -> a
error [Char]
"compositeSimpson: odd number of half-intervals" -- this should never happen
          combine [(Scalar b, b)
_] = forall v. AdditiveGroup v => v
zeroV
          combine ((Scalar b
t0,b
f0):(Scalar b
t1,b
f1):[(Scalar b, b)]
ps) = ((Scalar b
t1 forall a. Num a => a -> a -> a
- Scalar b
t0) forall a. Fractional a => a -> a -> a
/ Scalar b
2) forall v. VectorSpace v => Scalar v -> v -> v
*^ (b
f0 forall v. AdditiveGroup v => v -> v -> v
^+^ b
f1) forall v. AdditiveGroup v => v -> v -> v
^+^ [(Scalar b, b)] -> b
combine ((Scalar b
t1,b
f1)forall a. a -> [a] -> [a]
:[(Scalar b, b)]
ps)
      in forall {b}.
(VectorSpace b, Fractional (Scalar b)) =>
[(Scalar b, b)] -> b
combine [(Scalar v, v)]
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 :: forall v.
(VectorSpace v, Fractional (Scalar v)) =>
Int -> Scalar v -> Scalar v -> (Scalar v -> v) -> v
compositeSimpson Int
n Scalar v
a Scalar v
b Scalar v -> v
f
    = let nEven :: Int
nEven = Int
2 forall a. Num a => a -> a -> a
* forall a. Integral a => a -> a -> a
div Int
n Int
2
          dt :: Scalar v
dt = (Scalar v
b forall a. Num a => a -> a -> a
- Scalar v
a) forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
nEven
          ts :: [Scalar v]
ts = [Scalar v
a forall a. Num a => a -> a -> a
+ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
m forall a. Num a => a -> a -> a
* Scalar v
dt | Int
m <- [Int
0..Int
nEven]]
          pairs :: [(Scalar v, v)]
pairs = [(Scalar v
t,Scalar v -> v
f Scalar v
t) | Scalar v
t <- [Scalar v]
ts]
          combine :: [(Scalar b, b)] -> b
combine [] = forall a. HasCallStack => [Char] -> a
error [Char]
"compositeSimpson: odd number of half-intervals" -- this should never happen
          combine [(Scalar b, b)
_] = forall v. AdditiveGroup v => v
zeroV
          combine ((Scalar b, b)
_:(Scalar b, b)
_:[]) = forall a. HasCallStack => [Char] -> a
error [Char]
"compositeSimpson: odd number of half-intervals" -- this should never happen
          combine ((Scalar b
t0,b
f0):(Scalar b
_,b
f1):(Scalar b
t2,b
f2):[(Scalar b, b)]
ps) = ((Scalar b
t2 forall a. Num a => a -> a -> a
- Scalar b
t0) forall a. Fractional a => a -> a -> a
/ Scalar b
6) forall v. VectorSpace v => Scalar v -> v -> v
*^ (b
f0 forall v. AdditiveGroup v => v -> v -> v
^+^ Scalar b
4 forall v. VectorSpace v => Scalar v -> v -> v
*^ b
f1 forall v. AdditiveGroup v => v -> v -> v
^+^ b
f2) forall v. AdditiveGroup v => v -> v -> v
^+^ [(Scalar b, b)] -> b
combine ((Scalar b
t2,b
f2)forall a. a -> [a] -> [a]
:[(Scalar b, b)]
ps)
      in forall {b}.
(VectorSpace b, Fractional (Scalar b)) =>
[(Scalar b, b)] -> b
combine [(Scalar v, v)]
pairs