diff options
author | Alberto Ruiz <aruiz@um.es> | 2012-02-08 09:26:46 -0800 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2012-02-08 09:26:46 -0800 |
commit | 746427e5ec5d1f3fec2ec2ddec4f041f5563596a (patch) | |
tree | f38528b173720ed9b20f8c9adcf9a28b227170ab | |
parent | 3ac7bc8155e472b2c85966bd1d521d0739e67562 (diff) | |
parent | 1319c4ddae89c41633a4cd894ba28896d3445e79 (diff) |
Merge pull request #9 from cdhf/master
Support for integration over infinite intervals
-rw-r--r-- | lib/Numeric/GSL/Integration.hs | 104 | ||||
-rw-r--r-- | lib/Numeric/GSL/gsl-aux.c | 42 |
2 files changed, 145 insertions, 1 deletions
diff --git a/lib/Numeric/GSL/Integration.hs b/lib/Numeric/GSL/Integration.hs index a0e922b..2330fc6 100644 --- a/lib/Numeric/GSL/Integration.hs +++ b/lib/Numeric/GSL/Integration.hs | |||
@@ -17,7 +17,10 @@ Numerical integration routines. | |||
17 | 17 | ||
18 | module Numeric.GSL.Integration ( | 18 | module Numeric.GSL.Integration ( |
19 | integrateQNG, | 19 | integrateQNG, |
20 | integrateQAGS | 20 | integrateQAGS, |
21 | integrateQAGI, | ||
22 | integrateQAGIU, | ||
23 | integrateQAGIL | ||
21 | ) where | 24 | ) where |
22 | 25 | ||
23 | import Foreign.C.Types | 26 | import Foreign.C.Types |
@@ -94,3 +97,102 @@ integrateQNG prec f a b = unsafePerformIO $ do | |||
94 | foreign import ccall "gsl-aux.h integrate_qng" | 97 | foreign import ccall "gsl-aux.h integrate_qng" |
95 | c_integrate_qng :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> Double | 98 | c_integrate_qng :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> Double |
96 | -> Ptr Double -> Ptr Double -> IO CInt | 99 | -> Ptr Double -> Ptr Double -> IO CInt |
100 | |||
101 | -------------------------------------------------------------------- | ||
102 | {- | Numerical integration using /gsl_integration_qagi/ (integration over the infinite integral -Inf..Inf using QAGS). | ||
103 | For example: | ||
104 | |||
105 | @\> let quad = integrateQAGI 1E-9 1000 | ||
106 | \> let f a x = exp(-a * x^2) | ||
107 | \> quad (f 0.5) | ||
108 | (2.5066282746310002,6.229215880648858e-11)@ | ||
109 | |||
110 | -} | ||
111 | |||
112 | integrateQAGI :: Double -- ^ precision (e.g. 1E-9) | ||
113 | -> Int -- ^ size of auxiliary workspace (e.g. 1000) | ||
114 | -> (Double -> Double) -- ^ function to be integrated on the interval (-Inf,Inf) | ||
115 | -> (Double, Double) -- ^ result of the integration and error | ||
116 | integrateQAGI prec n f = unsafePerformIO $ do | ||
117 | r <- malloc | ||
118 | e <- malloc | ||
119 | fp <- mkfun (\x _ -> f x) | ||
120 | c_integrate_qagi fp prec (fromIntegral n) r e // check "integrate_qagi" | ||
121 | vr <- peek r | ||
122 | ve <- peek e | ||
123 | let result = (vr,ve) | ||
124 | free r | ||
125 | free e | ||
126 | freeHaskellFunPtr fp | ||
127 | return result | ||
128 | |||
129 | foreign import ccall "gsl-aux.h integrate_qagi" | ||
130 | c_integrate_qagi :: FunPtr (Double-> Ptr() -> Double) -> Double -> CInt | ||
131 | -> Ptr Double -> Ptr Double -> IO CInt | ||
132 | |||
133 | -------------------------------------------------------------------- | ||
134 | {- | Numerical integration using /gsl_integration_qagiu/ (integration over the semi-infinite integral a..Inf). | ||
135 | For example: | ||
136 | |||
137 | @\> let quad = integrateQAGIU 1E-9 1000 | ||
138 | \> let f a x = exp(-a * x^2) | ||
139 | \> quad (f 0.5) 0 | ||
140 | (1.2533141373155001,3.114607940324429e-11)@ | ||
141 | |||
142 | -} | ||
143 | |||
144 | integrateQAGIU :: Double -- ^ precision (e.g. 1E-9) | ||
145 | -> Int -- ^ size of auxiliary workspace (e.g. 1000) | ||
146 | -> (Double -> Double) -- ^ function to be integrated on the interval (a,Inf) | ||
147 | -> Double -- ^ a | ||
148 | -> (Double, Double) -- ^ result of the integration and error | ||
149 | integrateQAGIU prec n f a = unsafePerformIO $ do | ||
150 | r <- malloc | ||
151 | e <- malloc | ||
152 | fp <- mkfun (\x _ -> f x) | ||
153 | c_integrate_qagiu fp a prec (fromIntegral n) r e // check "integrate_qagiu" | ||
154 | vr <- peek r | ||
155 | ve <- peek e | ||
156 | let result = (vr,ve) | ||
157 | free r | ||
158 | free e | ||
159 | freeHaskellFunPtr fp | ||
160 | return result | ||
161 | |||
162 | foreign import ccall "gsl-aux.h integrate_qagiu" | ||
163 | c_integrate_qagiu :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> CInt | ||
164 | -> Ptr Double -> Ptr Double -> IO CInt | ||
165 | |||
166 | -------------------------------------------------------------------- | ||
167 | {- | Numerical integration using /gsl_integration_qagil/ (integration over the semi-infinite integral -Inf..b). | ||
168 | For example: | ||
169 | |||
170 | @\> let quad = integrateQAGIL 1E-9 1000 | ||
171 | \> let f a x = exp(-a * x^2) | ||
172 | \> quad (f 0.5) 0 | ||
173 | (1.2533141373155001,3.114607940324429e-11)@ | ||
174 | |||
175 | -} | ||
176 | |||
177 | integrateQAGIL :: Double -- ^ precision (e.g. 1E-9) | ||
178 | -> Int -- ^ size of auxiliary workspace (e.g. 1000) | ||
179 | -> (Double -> Double) -- ^ function to be integrated on the interval (a,Inf) | ||
180 | -> Double -- ^ b | ||
181 | -> (Double, Double) -- ^ result of the integration and error | ||
182 | integrateQAGIL prec n f b = unsafePerformIO $ do | ||
183 | r <- malloc | ||
184 | e <- malloc | ||
185 | fp <- mkfun (\x _ -> f x) | ||
186 | c_integrate_qagil fp b prec (fromIntegral n) r e // check "integrate_qagil" | ||
187 | vr <- peek r | ||
188 | ve <- peek e | ||
189 | let result = (vr,ve) | ||
190 | free r | ||
191 | free e | ||
192 | freeHaskellFunPtr fp | ||
193 | return result | ||
194 | |||
195 | foreign import ccall "gsl-aux.h integrate_qagil" | ||
196 | c_integrate_qagil :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> CInt | ||
197 | -> Ptr Double -> Ptr Double -> IO CInt | ||
198 | |||
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c index 33d7dab..5c17836 100644 --- a/lib/Numeric/GSL/gsl-aux.c +++ b/lib/Numeric/GSL/gsl-aux.c | |||
@@ -761,6 +761,48 @@ int integrate_qags(double f(double,void*), double a, double b, double prec, int | |||
761 | OK | 761 | OK |
762 | } | 762 | } |
763 | 763 | ||
764 | int integrate_qagi(double f(double,void*), double prec, int w, | ||
765 | double *result, double* error) { | ||
766 | DEBUGMSG("integrate_qagi"); | ||
767 | gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w); | ||
768 | gsl_function F; | ||
769 | F.function = f; | ||
770 | F.params = NULL; | ||
771 | int res = gsl_integration_qagi (&F, 0, prec, w,wk, result, error); | ||
772 | CHECK(res,res); | ||
773 | gsl_integration_workspace_free (wk); | ||
774 | OK | ||
775 | } | ||
776 | |||
777 | |||
778 | int integrate_qagiu(double f(double,void*), double a, double prec, int w, | ||
779 | double *result, double* error) { | ||
780 | DEBUGMSG("integrate_qagiu"); | ||
781 | gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w); | ||
782 | gsl_function F; | ||
783 | F.function = f; | ||
784 | F.params = NULL; | ||
785 | int res = gsl_integration_qagiu (&F, a, 0, prec, w,wk, result, error); | ||
786 | CHECK(res,res); | ||
787 | gsl_integration_workspace_free (wk); | ||
788 | OK | ||
789 | } | ||
790 | |||
791 | |||
792 | int integrate_qagil(double f(double,void*), double b, double prec, int w, | ||
793 | double *result, double* error) { | ||
794 | DEBUGMSG("integrate_qagil"); | ||
795 | gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w); | ||
796 | gsl_function F; | ||
797 | F.function = f; | ||
798 | F.params = NULL; | ||
799 | int res = gsl_integration_qagil (&F, b, 0, prec, w,wk, result, error); | ||
800 | CHECK(res,res); | ||
801 | gsl_integration_workspace_free (wk); | ||
802 | OK | ||
803 | } | ||
804 | |||
805 | |||
764 | int polySolve(KRVEC(a), CVEC(z)) { | 806 | int polySolve(KRVEC(a), CVEC(z)) { |
765 | DEBUGMSG("polySolve"); | 807 | DEBUGMSG("polySolve"); |
766 | REQUIRES(an>1,BAD_SIZE); | 808 | REQUIRES(an>1,BAD_SIZE); |