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