module TBit.Numerical.Derivative (diff, diff4) where

import Numeric.LinearAlgebra.HMatrix

-- |Takes the element-wise first derivative of the given vector/matrix-valued
--  function of a single variable. The first argument gives the spacing epsilon
--  as used in the difference formula, /i.e./ that thing which limit is taken
--  to zero. The second argument is the matrix function, and the third
--  argument is the point at which to evaluate the derivative. 
--
--  The derivative is taken using the finite difference method to second order.
diff :: (Num (c e),  Container c e) => e -> (e -> c e) -> e -> c e
diff e m x = sum
           [ scale (2.0/3.0/e)  (m (x+e)   - m (x-e))
           , scale (1.0/12.0/e) (m (x-2*e) - m (x+2*e)) ]

-- |As 'diff', but uses fourth order finite difference method. This means that
--  the matrix argument will need to be evaluated at twice as many points, which
--  increases the runtime twofold until we can implement a full-gridded calculation.
diff4 :: (Num (c e),  Container c e) => e -> (e -> c e) -> e -> c e
diff4 e m x = sum
            [ scale (4.0/5.0/e)   (m (x+e)   - m (x-e))
            , scale (1.0/5.0/e)   (m (x-2*e) - m (x+2*e))
            , scale (4.0/105.0/e) (m (x+3*e) - m (x-3*e))
            , scale (1.0/280.0/e) (m (x-4*e) - m (x+4*e))]