diff options
author | Matthew Peddie <peddie@alum.mit.edu> | 2013-08-19 01:26:08 -0700 |
---|---|---|
committer | Matthew Peddie <peddie@alum.mit.edu> | 2013-08-19 01:26:08 -0700 |
commit | 1661785e53574a501c7b40d10fabe808f49dd241 (patch) | |
tree | 3346e16f8f1df23b2f424bad8d6f1871c5b10e53 | |
parent | 374194ee454622e66c931dce8f315c68167c7ac9 (diff) |
Add bindings for gsl_integrate_cquad doubly-adaptive quadrature for difficult integrands.
-rw-r--r-- | lib/Numeric/GSL/Integration.hs | 51 | ||||
-rw-r--r-- | lib/Numeric/GSL/gsl-aux.c | 13 |
2 files changed, 63 insertions, 1 deletions
diff --git a/lib/Numeric/GSL/Integration.hs b/lib/Numeric/GSL/Integration.hs index defbf80..d1f552d 100644 --- a/lib/Numeric/GSL/Integration.hs +++ b/lib/Numeric/GSL/Integration.hs | |||
@@ -20,7 +20,8 @@ module Numeric.GSL.Integration ( | |||
20 | integrateQAGS, | 20 | integrateQAGS, |
21 | integrateQAGI, | 21 | integrateQAGI, |
22 | integrateQAGIU, | 22 | integrateQAGIU, |
23 | integrateQAGIL | 23 | integrateQAGIL, |
24 | integrateCQUAD | ||
24 | ) where | 25 | ) where |
25 | 26 | ||
26 | import Foreign.C.Types | 27 | import Foreign.C.Types |
@@ -196,3 +197,51 @@ foreign import ccall safe "gsl-aux.h integrate_qagil" | |||
196 | c_integrate_qagil :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> CInt | 197 | c_integrate_qagil :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> CInt |
197 | -> Ptr Double -> Ptr Double -> IO CInt | 198 | -> Ptr Double -> Ptr Double -> IO CInt |
198 | 199 | ||
200 | |||
201 | -------------------------------------------------------------------- | ||
202 | {- | Numerical integration using /gsl_integration_cquad/ (quadrature | ||
203 | for general integrands). From the GSL manual: | ||
204 | |||
205 | @CQUAD is a new doubly-adaptive general-purpose quadrature routine | ||
206 | which can handle most types of singularities, non-numerical function | ||
207 | values such as Inf or NaN, as well as some divergent integrals. It | ||
208 | generally requires more function evaluations than the integration | ||
209 | routines in QUADPACK, yet fails less often for difficult integrands.@ | ||
210 | |||
211 | For example: | ||
212 | |||
213 | @\> let quad = integrateCQUAD 1E-12 1000 | ||
214 | \> let f a x = exp(-a * x^2) | ||
215 | \> quad (f 0.5) 2 5 | ||
216 | (5.7025405463957006e-2,9.678874441303705e-16,95)@ | ||
217 | |||
218 | Unlike other quadrature methods, integrateCQUAD also returns the | ||
219 | number of function evaluations required. | ||
220 | |||
221 | -} | ||
222 | |||
223 | integrateCQUAD :: Double -- ^ precision (e.g. 1E-9) | ||
224 | -> Int -- ^ size of auxiliary workspace (e.g. 1000) | ||
225 | -> (Double -> Double) -- ^ function to be integrated on the interval (a, b) | ||
226 | -> Double -- ^ a | ||
227 | -> Double -- ^ b | ||
228 | -> (Double, Double, Int) -- ^ result of the integration, error and number of function evaluations performed | ||
229 | integrateCQUAD prec n f a b = unsafePerformIO $ do | ||
230 | r <- malloc | ||
231 | e <- malloc | ||
232 | neval <- malloc | ||
233 | fp <- mkfun (\x _ -> f x) | ||
234 | c_integrate_cquad fp a b prec (fromIntegral n) r e neval // check "integrate_cquad" | ||
235 | vr <- peek r | ||
236 | ve <- peek e | ||
237 | vneval <- peek neval | ||
238 | let result = (vr,ve,vneval) | ||
239 | free r | ||
240 | free e | ||
241 | free neval | ||
242 | freeHaskellFunPtr fp | ||
243 | return result | ||
244 | |||
245 | foreign import ccall safe "gsl-aux.h integrate_cquad" | ||
246 | c_integrate_cquad :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> Double -> CInt | ||
247 | -> Ptr Double -> Ptr Double -> Ptr Int -> IO CInt | ||
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c index e727c91..756edf1 100644 --- a/lib/Numeric/GSL/gsl-aux.c +++ b/lib/Numeric/GSL/gsl-aux.c | |||
@@ -802,6 +802,19 @@ int integrate_qagil(double f(double,void*), double b, double prec, int w, | |||
802 | OK | 802 | OK |
803 | } | 803 | } |
804 | 804 | ||
805 | int integrate_cquad(double f(double,void*), double a, double b, double prec, | ||
806 | int w, double *result, double* error, int *neval) { | ||
807 | DEBUGMSG("integrate_cquad"); | ||
808 | gsl_integration_cquad_workspace * wk = gsl_integration_cquad_workspace_alloc (w); | ||
809 | gsl_function F; | ||
810 | F.function = f; | ||
811 | F.params = NULL; | ||
812 | int res = gsl_integration_cquad (&F, a, b, 0, prec, wk, result, error, neval); | ||
813 | CHECK(res,res); | ||
814 | gsl_integration_cquad_workspace_free (wk); | ||
815 | OK | ||
816 | } | ||
817 | |||
805 | 818 | ||
806 | int polySolve(KRVEC(a), CVEC(z)) { | 819 | int polySolve(KRVEC(a), CVEC(z)) { |
807 | DEBUGMSG("polySolve"); | 820 | DEBUGMSG("polySolve"); |