summaryrefslogtreecommitdiff
path: root/packages/hmatrix/src/Numeric/GSL/Differentiation.hs
diff options
context:
space:
mode:
authorAlberto Ruiz <aruiz@um.es>2014-05-08 08:48:12 +0200
committerAlberto Ruiz <aruiz@um.es>2014-05-08 08:48:12 +0200
commit1925c123d7d8184a1d2ddc0a413e0fd2776e1083 (patch)
treefad79f909d9c3be53d68e6ebd67202650536d387 /packages/hmatrix/src/Numeric/GSL/Differentiation.hs
parenteb3f702d065a4a967bb754977233e6eec408fd1f (diff)
empty hmatrix-base
Diffstat (limited to 'packages/hmatrix/src/Numeric/GSL/Differentiation.hs')
-rw-r--r--packages/hmatrix/src/Numeric/GSL/Differentiation.hs87
1 files changed, 87 insertions, 0 deletions
diff --git a/packages/hmatrix/src/Numeric/GSL/Differentiation.hs b/packages/hmatrix/src/Numeric/GSL/Differentiation.hs
new file mode 100644
index 0000000..93c5007
--- /dev/null
+++ b/packages/hmatrix/src/Numeric/GSL/Differentiation.hs
@@ -0,0 +1,87 @@
1{-# OPTIONS #-}
2-----------------------------------------------------------------------------
3{- |
4Module : Numeric.GSL.Differentiation
5Copyright : (c) Alberto Ruiz 2006
6License : GPL-style
7
8Maintainer : Alberto Ruiz (aruiz at um dot es)
9Stability : provisional
10Portability : uses ffi
11
12Numerical differentiation.
13
14<http://www.gnu.org/software/gsl/manual/html_node/Numerical-Differentiation.html#Numerical-Differentiation>
15
16From the GSL manual: \"The functions described in this chapter compute numerical derivatives by finite differencing. An adaptive algorithm is used to find the best choice of finite difference and to estimate the error in the derivative.\"
17-}
18-----------------------------------------------------------------------------
19module Numeric.GSL.Differentiation (
20 derivCentral,
21 derivForward,
22 derivBackward
23) where
24
25import Foreign.C.Types
26import Foreign.Marshal.Alloc(malloc, free)
27import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr)
28import Foreign.Storable(peek)
29import Data.Packed.Internal(check,(//))
30import System.IO.Unsafe(unsafePerformIO)
31
32derivGen ::
33 CInt -- ^ type: 0 central, 1 forward, 2 backward
34 -> Double -- ^ initial step size
35 -> (Double -> Double) -- ^ function
36 -> Double -- ^ point where the derivative is taken
37 -> (Double, Double) -- ^ result and error
38derivGen c h f x = unsafePerformIO $ do
39 r <- malloc
40 e <- malloc
41 fp <- mkfun (\y _ -> f y)
42 c_deriv c fp x h r e // check "deriv"
43 vr <- peek r
44 ve <- peek e
45 let result = (vr,ve)
46 free r
47 free e
48 freeHaskellFunPtr fp
49 return result
50
51foreign import ccall safe "gsl-aux.h deriv"
52 c_deriv :: CInt -> FunPtr (Double -> Ptr () -> Double) -> Double -> Double
53 -> Ptr Double -> Ptr Double -> IO CInt
54
55
56{- | Adaptive central difference algorithm, /gsl_deriv_central/. For example:
57
58>>> let deriv = derivCentral 0.01
59>>> deriv sin (pi/4)
60(0.7071067812000676,1.0600063101654055e-10)
61>>> cos (pi/4)
620.7071067811865476
63
64-}
65derivCentral :: Double -- ^ initial step size
66 -> (Double -> Double) -- ^ function
67 -> Double -- ^ point where the derivative is taken
68 -> (Double, Double) -- ^ result and absolute error
69derivCentral = derivGen 0
70
71-- | Adaptive forward difference algorithm, /gsl_deriv_forward/. The function is evaluated only at points greater than x, and never at x itself. The derivative is returned in result and an estimate of its absolute error is returned in abserr. This function should be used if f(x) has a discontinuity at x, or is undefined for values less than x. A backward derivative can be obtained using a negative step.
72derivForward :: Double -- ^ initial step size
73 -> (Double -> Double) -- ^ function
74 -> Double -- ^ point where the derivative is taken
75 -> (Double, Double) -- ^ result and absolute error
76derivForward = derivGen 1
77
78-- | Adaptive backward difference algorithm, /gsl_deriv_backward/.
79derivBackward ::Double -- ^ initial step size
80 -> (Double -> Double) -- ^ function
81 -> Double -- ^ point where the derivative is taken
82 -> (Double, Double) -- ^ result and absolute error
83derivBackward = derivGen 2
84
85{- | conversion of Haskell functions into function pointers that can be used in the C side
86-}
87foreign import ccall safe "wrapper" mkfun:: (Double -> Ptr() -> Double) -> IO( FunPtr (Double -> Ptr() -> Double))