summaryrefslogtreecommitdiff
path: root/lib/Numeric/GSL/Integration.hs
diff options
context:
space:
mode:
Diffstat (limited to 'lib/Numeric/GSL/Integration.hs')
-rw-r--r--lib/Numeric/GSL/Integration.hs250
1 files changed, 0 insertions, 250 deletions
diff --git a/lib/Numeric/GSL/Integration.hs b/lib/Numeric/GSL/Integration.hs
deleted file mode 100644
index 5f0a415..0000000
--- a/lib/Numeric/GSL/Integration.hs
+++ /dev/null
@@ -1,250 +0,0 @@
1{-# OPTIONS #-}
2-----------------------------------------------------------------------------
3{- |
4Module : Numeric.GSL.Integration
5Copyright : (c) Alberto Ruiz 2006
6License : GPL-style
7
8Maintainer : Alberto Ruiz (aruiz at um dot es)
9Stability : provisional
10Portability : uses ffi
11
12Numerical integration routines.
13
14<http://www.gnu.org/software/gsl/manual/html_node/Numerical-Integration.html#Numerical-Integration>
15-}
16-----------------------------------------------------------------------------
17
18module Numeric.GSL.Integration (
19 integrateQNG,
20 integrateQAGS,
21 integrateQAGI,
22 integrateQAGIU,
23 integrateQAGIL,
24 integrateCQUAD
25) where
26
27import Foreign.C.Types
28import Foreign.Marshal.Alloc(malloc, free)
29import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr)
30import Foreign.Storable(peek)
31import Data.Packed.Internal(check,(//))
32import System.IO.Unsafe(unsafePerformIO)
33
34eps = 1e-12
35
36{- | conversion of Haskell functions into function pointers that can be used in the C side
37-}
38foreign import ccall safe "wrapper" mkfun:: (Double -> Ptr() -> Double) -> IO( FunPtr (Double -> Ptr() -> Double))
39
40--------------------------------------------------------------------
41{- | Numerical integration using /gsl_integration_qags/ (adaptive integration with singularities). For example:
42
43>>> let quad = integrateQAGS 1E-9 1000
44>>> let f a x = x**(-0.5) * log (a*x)
45>>> quad (f 1) 0 1
46(-3.999999999999974,4.871658632055187e-13)
47
48-}
49
50integrateQAGS :: Double -- ^ precision (e.g. 1E-9)
51 -> Int -- ^ size of auxiliary workspace (e.g. 1000)
52 -> (Double -> Double) -- ^ function to be integrated on the interval (a,b)
53 -> Double -- ^ a
54 -> Double -- ^ b
55 -> (Double, Double) -- ^ result of the integration and error
56integrateQAGS prec n f a b = unsafePerformIO $ do
57 r <- malloc
58 e <- malloc
59 fp <- mkfun (\x _ -> f x)
60 c_integrate_qags fp a b eps prec (fromIntegral n) r e // check "integrate_qags"
61 vr <- peek r
62 ve <- peek e
63 let result = (vr,ve)
64 free r
65 free e
66 freeHaskellFunPtr fp
67 return result
68
69foreign import ccall safe "integrate_qags" c_integrate_qags
70 :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
71 -> Double -> Double -> CInt -> Ptr Double -> Ptr Double -> IO CInt
72
73-----------------------------------------------------------------
74{- | Numerical integration using /gsl_integration_qng/ (useful for fast integration of smooth functions). For example:
75
76>>> let quad = integrateQNG 1E-6
77>>> quad (\x -> 4/(1+x*x)) 0 1
78(3.141592653589793,3.487868498008632e-14)
79
80-}
81
82integrateQNG :: Double -- ^ precision (e.g. 1E-9)
83 -> (Double -> Double) -- ^ function to be integrated on the interval (a,b)
84 -> Double -- ^ a
85 -> Double -- ^ b
86 -> (Double, Double) -- ^ result of the integration and error
87integrateQNG prec f a b = unsafePerformIO $ do
88 r <- malloc
89 e <- malloc
90 fp <- mkfun (\x _ -> f x)
91 c_integrate_qng fp a b eps prec r e // check "integrate_qng"
92 vr <- peek r
93 ve <- peek e
94 let result = (vr,ve)
95 free r
96 free e
97 freeHaskellFunPtr fp
98 return result
99
100foreign import ccall safe "integrate_qng" c_integrate_qng
101 :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
102 -> Double -> Double -> Ptr Double -> Ptr Double -> IO CInt
103
104--------------------------------------------------------------------
105{- | Numerical integration using /gsl_integration_qagi/ (integration over the infinite integral -Inf..Inf using QAGS).
106For example:
107
108>>> let quad = integrateQAGI 1E-9 1000
109>>> let f a x = exp(-a * x^2)
110>>> quad (f 0.5)
111(2.5066282746310002,6.229215880648858e-11)
112
113-}
114
115integrateQAGI :: Double -- ^ precision (e.g. 1E-9)
116 -> Int -- ^ size of auxiliary workspace (e.g. 1000)
117 -> (Double -> Double) -- ^ function to be integrated on the interval (-Inf,Inf)
118 -> (Double, Double) -- ^ result of the integration and error
119integrateQAGI prec n f = unsafePerformIO $ do
120 r <- malloc
121 e <- malloc
122 fp <- mkfun (\x _ -> f x)
123 c_integrate_qagi fp eps prec (fromIntegral n) r e // check "integrate_qagi"
124 vr <- peek r
125 ve <- peek e
126 let result = (vr,ve)
127 free r
128 free e
129 freeHaskellFunPtr fp
130 return result
131
132foreign import ccall safe "integrate_qagi" c_integrate_qagi
133 :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
134 -> CInt -> Ptr Double -> Ptr Double -> IO CInt
135
136--------------------------------------------------------------------
137{- | Numerical integration using /gsl_integration_qagiu/ (integration over the semi-infinite integral a..Inf).
138For example:
139
140>>> let quad = integrateQAGIU 1E-9 1000
141>>> let f a x = exp(-a * x^2)
142>>> quad (f 0.5) 0
143(1.2533141373155001,3.114607940324429e-11)
144
145-}
146
147integrateQAGIU :: Double -- ^ precision (e.g. 1E-9)
148 -> Int -- ^ size of auxiliary workspace (e.g. 1000)
149 -> (Double -> Double) -- ^ function to be integrated on the interval (a,Inf)
150 -> Double -- ^ a
151 -> (Double, Double) -- ^ result of the integration and error
152integrateQAGIU prec n f a = unsafePerformIO $ do
153 r <- malloc
154 e <- malloc
155 fp <- mkfun (\x _ -> f x)
156 c_integrate_qagiu fp a eps prec (fromIntegral n) r e // check "integrate_qagiu"
157 vr <- peek r
158 ve <- peek e
159 let result = (vr,ve)
160 free r
161 free e
162 freeHaskellFunPtr fp
163 return result
164
165foreign import ccall safe "integrate_qagiu" c_integrate_qagiu
166 :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
167 -> Double -> CInt -> Ptr Double -> Ptr Double -> IO CInt
168
169--------------------------------------------------------------------
170{- | Numerical integration using /gsl_integration_qagil/ (integration over the semi-infinite integral -Inf..b).
171For example:
172
173>>> let quad = integrateQAGIL 1E-9 1000
174>>> let f a x = exp(-a * x^2)
175>>> quad (f 0.5) 0
176(1.2533141373155001,3.114607940324429e-11)
177
178-}
179
180integrateQAGIL :: Double -- ^ precision (e.g. 1E-9)
181 -> Int -- ^ size of auxiliary workspace (e.g. 1000)
182 -> (Double -> Double) -- ^ function to be integrated on the interval (a,Inf)
183 -> Double -- ^ b
184 -> (Double, Double) -- ^ result of the integration and error
185integrateQAGIL prec n f b = unsafePerformIO $ do
186 r <- malloc
187 e <- malloc
188 fp <- mkfun (\x _ -> f x)
189 c_integrate_qagil fp b eps prec (fromIntegral n) r e // check "integrate_qagil"
190 vr <- peek r
191 ve <- peek e
192 let result = (vr,ve)
193 free r
194 free e
195 freeHaskellFunPtr fp
196 return result
197
198foreign import ccall safe "gsl-aux.h integrate_qagil" c_integrate_qagil
199 :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
200 -> Double -> CInt -> Ptr Double -> Ptr Double -> IO CInt
201
202
203--------------------------------------------------------------------
204{- | Numerical integration using /gsl_integration_cquad/ (quadrature
205for general integrands). From the GSL manual:
206
207@CQUAD is a new doubly-adaptive general-purpose quadrature routine
208which can handle most types of singularities, non-numerical function
209values such as Inf or NaN, as well as some divergent integrals. It
210generally requires more function evaluations than the integration
211routines in QUADPACK, yet fails less often for difficult integrands.@
212
213For example:
214
215>>> let quad = integrateCQUAD 1E-12 1000
216>>> let f a x = exp(-a * x^2)
217>>> quad (f 0.5) 2 5
218(5.7025405463957006e-2,9.678874441303705e-16,95)
219
220Unlike other quadrature methods, integrateCQUAD also returns the
221number of function evaluations required.
222
223-}
224
225integrateCQUAD :: Double -- ^ precision (e.g. 1E-9)
226 -> Int -- ^ size of auxiliary workspace (e.g. 1000)
227 -> (Double -> Double) -- ^ function to be integrated on the interval (a, b)
228 -> Double -- ^ a
229 -> Double -- ^ b
230 -> (Double, Double, Int) -- ^ result of the integration, error and number of function evaluations performed
231integrateCQUAD prec n f a b = unsafePerformIO $ do
232 r <- malloc
233 e <- malloc
234 neval <- malloc
235 fp <- mkfun (\x _ -> f x)
236 c_integrate_cquad fp a b eps prec (fromIntegral n) r e neval // check "integrate_cquad"
237 vr <- peek r
238 ve <- peek e
239 vneval <- peek neval
240 let result = (vr,ve,vneval)
241 free r
242 free e
243 free neval
244 freeHaskellFunPtr fp
245 return result
246
247foreign import ccall safe "integrate_cquad" c_integrate_cquad
248 :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
249 -> Double -> Double -> CInt -> Ptr Double -> Ptr Double -> Ptr Int -> IO CInt
250
generated by cgit v1.2.3 (git 2.39.1) at 2024-09-15 09:19:34 +0000