From 1319c4ddae89c41633a4cd894ba28896d3445e79 Mon Sep 17 00:00:00 2001 From: Danny Chan Date: Mon, 6 Feb 2012 21:24:10 +0100 Subject: add support for gsl_integration_qagi, gsl_integration_qagiu, gsl_integration_qagil --- lib/Numeric/GSL/Integration.hs | 104 ++++++++++++++++++++++++++++++++++++++++- lib/Numeric/GSL/gsl-aux.c | 42 +++++++++++++++++ 2 files changed, 145 insertions(+), 1 deletion(-) 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. module Numeric.GSL.Integration ( integrateQNG, - integrateQAGS + integrateQAGS, + integrateQAGI, + integrateQAGIU, + integrateQAGIL ) where import Foreign.C.Types @@ -94,3 +97,102 @@ integrateQNG prec f a b = unsafePerformIO $ do foreign import ccall "gsl-aux.h integrate_qng" c_integrate_qng :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> Double -> Ptr Double -> Ptr Double -> IO CInt + +-------------------------------------------------------------------- +{- | Numerical integration using /gsl_integration_qagi/ (integration over the infinite integral -Inf..Inf using QAGS). +For example: + +@\> let quad = integrateQAGI 1E-9 1000 +\> let f a x = exp(-a * x^2) +\> quad (f 0.5) +(2.5066282746310002,6.229215880648858e-11)@ + +-} + +integrateQAGI :: Double -- ^ precision (e.g. 1E-9) + -> Int -- ^ size of auxiliary workspace (e.g. 1000) + -> (Double -> Double) -- ^ function to be integrated on the interval (-Inf,Inf) + -> (Double, Double) -- ^ result of the integration and error +integrateQAGI prec n f = unsafePerformIO $ do + r <- malloc + e <- malloc + fp <- mkfun (\x _ -> f x) + c_integrate_qagi fp prec (fromIntegral n) r e // check "integrate_qagi" + vr <- peek r + ve <- peek e + let result = (vr,ve) + free r + free e + freeHaskellFunPtr fp + return result + +foreign import ccall "gsl-aux.h integrate_qagi" + c_integrate_qagi :: FunPtr (Double-> Ptr() -> Double) -> Double -> CInt + -> Ptr Double -> Ptr Double -> IO CInt + +-------------------------------------------------------------------- +{- | Numerical integration using /gsl_integration_qagiu/ (integration over the semi-infinite integral a..Inf). +For example: + +@\> let quad = integrateQAGIU 1E-9 1000 +\> let f a x = exp(-a * x^2) +\> quad (f 0.5) 0 +(1.2533141373155001,3.114607940324429e-11)@ + +-} + +integrateQAGIU :: Double -- ^ precision (e.g. 1E-9) + -> Int -- ^ size of auxiliary workspace (e.g. 1000) + -> (Double -> Double) -- ^ function to be integrated on the interval (a,Inf) + -> Double -- ^ a + -> (Double, Double) -- ^ result of the integration and error +integrateQAGIU prec n f a = unsafePerformIO $ do + r <- malloc + e <- malloc + fp <- mkfun (\x _ -> f x) + c_integrate_qagiu fp a prec (fromIntegral n) r e // check "integrate_qagiu" + vr <- peek r + ve <- peek e + let result = (vr,ve) + free r + free e + freeHaskellFunPtr fp + return result + +foreign import ccall "gsl-aux.h integrate_qagiu" + c_integrate_qagiu :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> CInt + -> Ptr Double -> Ptr Double -> IO CInt + +-------------------------------------------------------------------- +{- | Numerical integration using /gsl_integration_qagil/ (integration over the semi-infinite integral -Inf..b). +For example: + +@\> let quad = integrateQAGIL 1E-9 1000 +\> let f a x = exp(-a * x^2) +\> quad (f 0.5) 0 +(1.2533141373155001,3.114607940324429e-11)@ + +-} + +integrateQAGIL :: Double -- ^ precision (e.g. 1E-9) + -> Int -- ^ size of auxiliary workspace (e.g. 1000) + -> (Double -> Double) -- ^ function to be integrated on the interval (a,Inf) + -> Double -- ^ b + -> (Double, Double) -- ^ result of the integration and error +integrateQAGIL prec n f b = unsafePerformIO $ do + r <- malloc + e <- malloc + fp <- mkfun (\x _ -> f x) + c_integrate_qagil fp b prec (fromIntegral n) r e // check "integrate_qagil" + vr <- peek r + ve <- peek e + let result = (vr,ve) + free r + free e + freeHaskellFunPtr fp + return result + +foreign import ccall "gsl-aux.h integrate_qagil" + c_integrate_qagil :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> CInt + -> Ptr Double -> Ptr Double -> IO CInt + 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 OK } +int integrate_qagi(double f(double,void*), double prec, int w, + double *result, double* error) { + DEBUGMSG("integrate_qagi"); + gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w); + gsl_function F; + F.function = f; + F.params = NULL; + int res = gsl_integration_qagi (&F, 0, prec, w,wk, result, error); + CHECK(res,res); + gsl_integration_workspace_free (wk); + OK +} + + +int integrate_qagiu(double f(double,void*), double a, double prec, int w, + double *result, double* error) { + DEBUGMSG("integrate_qagiu"); + gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w); + gsl_function F; + F.function = f; + F.params = NULL; + int res = gsl_integration_qagiu (&F, a, 0, prec, w,wk, result, error); + CHECK(res,res); + gsl_integration_workspace_free (wk); + OK +} + + +int integrate_qagil(double f(double,void*), double b, double prec, int w, + double *result, double* error) { + DEBUGMSG("integrate_qagil"); + gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w); + gsl_function F; + F.function = f; + F.params = NULL; + int res = gsl_integration_qagil (&F, b, 0, prec, w,wk, result, error); + CHECK(res,res); + gsl_integration_workspace_free (wk); + OK +} + + int polySolve(KRVEC(a), CVEC(z)) { DEBUGMSG("polySolve"); REQUIRES(an>1,BAD_SIZE); -- cgit v1.2.3 From 1d29051bc2ee1838a326d2d73f4bcf5907d4ebfd Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Wed, 8 Feb 2012 18:35:03 +0100 Subject: thanks, changes, version --- CHANGES.md | 5 +++++ THANKS.md | 2 ++ hmatrix.cabal | 2 +- 3 files changed, 8 insertions(+), 1 deletion(-) diff --git a/CHANGES.md b/CHANGES.md index 99b8e96..de466e5 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -1,3 +1,8 @@ +0.13.2.0 +-------- + +- integration over infinite intervals + 0.13.0.0 -------- diff --git a/THANKS.md b/THANKS.md index e2eaa97..23c93bb 100644 --- a/THANKS.md +++ b/THANKS.md @@ -107,3 +107,5 @@ module reorganization, monadic mapVectorM, and many other improvements. - Daniel Fischer reported some Haddock markup errors. +- Danny Chan added support for integration over infinite intervals. + diff --git a/hmatrix.cabal b/hmatrix.cabal index de1d6f7..b169f26 100644 --- a/hmatrix.cabal +++ b/hmatrix.cabal @@ -1,5 +1,5 @@ Name: hmatrix -Version: 0.13.1.0 +Version: 0.13.2.0 License: GPL License-file: LICENSE Author: Alberto Ruiz -- cgit v1.2.3 From 9322df983400894625105f66e36e5718329c0053 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Thu, 23 Feb 2012 19:08:16 +0100 Subject: minor doc fix --- lib/Numeric/LinearAlgebra/Algorithms.hs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs index e2ecd4d..435cc5a 100644 --- a/lib/Numeric/LinearAlgebra/Algorithms.hs +++ b/lib/Numeric/LinearAlgebra/Algorithms.hs @@ -209,7 +209,7 @@ rightSV :: Field t => Matrix t -> (Vector Double, Matrix t) rightSV m | vertical m = let (_,s,v) = thinSVD m in (s,v) | otherwise = let (_,s,v) = svd m in (s,v) --- | Singular values and all right singular vectors. +-- | Singular values and all left singular vectors. leftSV :: Field t => Matrix t -> (Matrix t, Vector Double) leftSV m | vertical m = let (u,s,_) = svd m in (u,s) | otherwise = let (u,s,_) = thinSVD m in (u,s) -- cgit v1.2.3 From eab231907c3d11651e48e15510f4510fd6c77450 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Sat, 25 Feb 2012 12:40:49 +0100 Subject: merge changelog --- CHANGES.md | 3 +++ examples/ode.hs | 16 ++++++++++++ lib/Numeric/GSL/ODE.hs | 47 ++++++++++++++++++++++++++--------- lib/Numeric/GSL/gsl-aux.c | 62 ++++++++++++++++++++++++----------------------- 4 files changed, 87 insertions(+), 41 deletions(-) diff --git a/CHANGES.md b/CHANGES.md index de466e5..10be500 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -3,6 +3,9 @@ - integration over infinite intervals +- msadams and msbdf methods for ode + + 0.13.0.0 -------- diff --git a/examples/ode.hs b/examples/ode.hs index 082c46c..dc6e0ec 100644 --- a/examples/ode.hs +++ b/examples/ode.hs @@ -1,6 +1,9 @@ +{-# LANGUAGE ViewPatterns #-} import Numeric.GSL.ODE import Numeric.LinearAlgebra import Graphics.Plot +import Debug.Trace(trace) +debug x = trace (show x) x vanderpol mu = do let xdot mu t [x,v] = [v, -x + mu * v * (1-x^2)] @@ -32,3 +35,16 @@ main = do harmonic 1 0.1 kepler 0.3 60 kepler 0.4 70 + vanderpol' 2 + +-- example of odeSolveV with jacobian +vanderpol' mu = do + let xdot mu t (toList->[x,v]) = fromList [v, -x + mu * v * (1-x^2)] + jac t (toList->[x,v]) = (2><2) [ 0 , 1 + , -1-2*x*v*mu, mu*(1-x**2) ] + ts = linspace 1000 (0,50) + hi = (ts@>1 - ts@>0)/100 + sol = toColumns $ odeSolveV (MSBDF jac) hi 1E-8 1E-8 (xdot mu) (fromList [1,0]) ts + mplot sol + + diff --git a/lib/Numeric/GSL/ODE.hs b/lib/Numeric/GSL/ODE.hs index 2251acd..c243f4b 100644 --- a/lib/Numeric/GSL/ODE.hs +++ b/lib/Numeric/GSL/ODE.hs @@ -29,7 +29,7 @@ main = mplot (ts : toColumns sol)@ ----------------------------------------------------------------------------- module Numeric.GSL.ODE ( - odeSolve, odeSolveV, ODEMethod(..) + odeSolve, odeSolveV, ODEMethod(..), Jacobian ) where import Data.Packed.Internal @@ -41,18 +41,21 @@ import System.IO.Unsafe(unsafePerformIO) ------------------------------------------------------------------------- +type Jacobian = Double -> Vector Double -> Matrix Double + -- | Stepping functions data ODEMethod = RK2 -- ^ Embedded Runge-Kutta (2, 3) method. - | RK4 -- ^ 4th order (classical) Runge-Kutta. The error estimate is obtained by halving the step-size. For more efficient estimate of the error, use 'RKf45'. + | RK4 -- ^ 4th order (classical) Runge-Kutta. The error estimate is obtained by halving the step-size. For more efficient estimate of the error, use the embedded methods. | RKf45 -- ^ Embedded Runge-Kutta-Fehlberg (4, 5) method. This method is a good general-purpose integrator. | RKck -- ^ Embedded Runge-Kutta Cash-Karp (4, 5) method. | RK8pd -- ^ Embedded Runge-Kutta Prince-Dormand (8,9) method. - | RK2imp -- ^ Implicit 2nd order Runge-Kutta at Gaussian points. - | RK4imp -- ^ Implicit 4th order Runge-Kutta at Gaussian points. - | BSimp -- ^ Implicit Bulirsch-Stoer method of Bader and Deuflhard. This algorithm requires the Jacobian. - | Gear1 -- ^ M=1 implicit Gear method. - | Gear2 -- ^ M=2 implicit Gear method. - deriving (Enum,Eq,Show,Bounded) + | RK2imp Jacobian -- ^ Implicit 2nd order Runge-Kutta at Gaussian points. + | RK4imp Jacobian -- ^ Implicit 4th order Runge-Kutta at Gaussian points. + | BSimp Jacobian -- ^ Implicit Bulirsch-Stoer method of Bader and Deuflhard. The method is generally suitable for stiff problems. + | RK1imp Jacobian -- ^ Implicit Gaussian first order Runge-Kutta. Also known as implicit Euler or backward Euler method. Error estimation is carried out by the step doubling method. + | MSAdams -- ^ A variable-coefficient linear multistep Adams method in Nordsieck form. This stepper uses explicit Adams-Bashforth (predictor) and implicit Adams-Moulton (corrector) methods in P(EC)^m functional iteration mode. Method order varies dynamically between 1 and 12. + | MSBDF Jacobian -- ^ A variable-coefficient linear multistep backward differentiation formula (BDF) method in Nordsieck form. This stepper uses the explicit BDF formula as predictor and implicit BDF formula as corrector. A modified Newton iteration method is used to solve the system of non-linear equations. Method order varies dynamically between 1 and 5. The method is generally suitable for stiff problems. + -- | A version of 'odeSolveV' with reasonable default parameters and system of equations defined using lists. odeSolve @@ -60,7 +63,7 @@ odeSolve -> [Double] -- ^ initial conditions -> Vector Double -- ^ desired solution times -> Matrix Double -- ^ solution -odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) Nothing (fromList xi) ts +odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) (fromList xi) ts where hi = (ts@>1 - ts@>0)/100 epsAbs = 1.49012e-08 epsRel = 1.49012e-08 @@ -73,11 +76,33 @@ odeSolveV -> Double -- ^ absolute tolerance for the state vector -> Double -- ^ relative tolerance for the state vector -> (Double -> Vector Double -> Vector Double) -- ^ xdot(t,x) + -> Vector Double -- ^ initial conditions + -> Vector Double -- ^ desired solution times + -> Matrix Double -- ^ solution +odeSolveV RK2 = odeSolveV' 0 Nothing +odeSolveV RK4 = odeSolveV' 1 Nothing +odeSolveV RKf45 = odeSolveV' 2 Nothing +odeSolveV RKck = odeSolveV' 3 Nothing +odeSolveV RK8pd = odeSolveV' 4 Nothing +odeSolveV (RK2imp jac) = odeSolveV' 5 (Just jac) +odeSolveV (RK4imp jac) = odeSolveV' 6 (Just jac) +odeSolveV (BSimp jac) = odeSolveV' 7 (Just jac) +odeSolveV (RK1imp jac) = odeSolveV' 8 (Just jac) +odeSolveV MSAdams = odeSolveV' 9 Nothing +odeSolveV (MSBDF jac) = odeSolveV' 10 (Just jac) + + +odeSolveV' + :: CInt -> Maybe (Double -> Vector Double -> Matrix Double) -- ^ optional jacobian + -> Double -- ^ initial step size + -> Double -- ^ absolute tolerance for the state vector + -> Double -- ^ relative tolerance for the state vector + -> (Double -> Vector Double -> Vector Double) -- ^ xdot(t,x) -> Vector Double -- ^ initial conditions -> Vector Double -- ^ desired solution times -> Matrix Double -- ^ solution -odeSolveV method h epsAbs epsRel f mbjac xiv ts = unsafePerformIO $ do +odeSolveV' method mbjac h epsAbs epsRel f xiv ts = unsafePerformIO $ do let n = dim xiv fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t)) jp <- case mbjac of @@ -86,7 +111,7 @@ odeSolveV method h epsAbs epsRel f mbjac xiv ts = unsafePerformIO $ do sol <- vec xiv $ \xiv' -> vec (checkTimes ts) $ \ts' -> createMIO (dim ts) n - (ode_c (fi (fromEnum method)) h epsAbs epsRel fp jp // xiv' // ts' ) + (ode_c (method) h epsAbs epsRel fp jp // xiv' // ts' ) "ode" freeHaskellFunPtr fp return sol diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c index 5c17836..7f1cf68 100644 --- a/lib/Numeric/GSL/gsl-aux.c +++ b/lib/Numeric/GSL/gsl-aux.c @@ -32,7 +32,7 @@ #include #include #include -#include +#include #include #include #include @@ -1356,38 +1356,38 @@ int ode(int method, double h, double eps_abs, double eps_rel, int jac(double, int, const double*, int, int, double*), KRVEC(xi), KRVEC(ts), RMAT(sol)) { - const gsl_odeiv_step_type * T; + const gsl_odeiv2_step_type * T; switch(method) { - case 0 : {T = gsl_odeiv_step_rk2; break; } - case 1 : {T = gsl_odeiv_step_rk4; break; } - case 2 : {T = gsl_odeiv_step_rkf45; break; } - case 3 : {T = gsl_odeiv_step_rkck; break; } - case 4 : {T = gsl_odeiv_step_rk8pd; break; } - case 5 : {T = gsl_odeiv_step_rk2imp; break; } - case 6 : {T = gsl_odeiv_step_rk4imp; break; } - case 7 : {T = gsl_odeiv_step_bsimp; break; } - case 8 : {T = gsl_odeiv_step_gear1; break; } - case 9 : {T = gsl_odeiv_step_gear2; break; } + case 0 : {T = gsl_odeiv2_step_rk2; break; } + case 1 : {T = gsl_odeiv2_step_rk4; break; } + case 2 : {T = gsl_odeiv2_step_rkf45; break; } + case 3 : {T = gsl_odeiv2_step_rkck; break; } + case 4 : {T = gsl_odeiv2_step_rk8pd; break; } + case 5 : {T = gsl_odeiv2_step_rk2imp; break; } + case 6 : {T = gsl_odeiv2_step_rk4imp; break; } + case 7 : {T = gsl_odeiv2_step_bsimp; break; } + case 8 : {T = gsl_odeiv2_step_rk1imp; break; } + case 9 : {T = gsl_odeiv2_step_msadams; break; } + case 10: {T = gsl_odeiv2_step_msbdf; break; } default: ERROR(BAD_CODE); } - - gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, xin); - gsl_odeiv_control * c = gsl_odeiv_control_y_new (eps_abs, eps_rel); - gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (xin); - Tode P; P.f = f; P.j = jac; P.n = xin; - gsl_odeiv_system sys = {odefunc, odejac, xin, &P}; + gsl_odeiv2_system sys = {odefunc, odejac, xin, &P}; + + gsl_odeiv2_driver * d = + gsl_odeiv2_driver_alloc_y_new (&sys, T, h, eps_abs, eps_rel); double t = tsp[0]; double* y = (double*)calloc(xin,sizeof(double)); int i,j; + int status; for(i=0; i< xin; i++) { y[i] = xip[i]; solp[i] = xip[i]; @@ -1396,22 +1396,24 @@ int ode(int method, double h, double eps_abs, double eps_rel, for (i = 1; i < tsn ; i++) { double ti = tsp[i]; - while (t < ti) - { - gsl_odeiv_evolve_apply (e, c, s, - &sys, - &t, ti, &h, - y); - // if (h < hmin) h = hmin; - } + + status = gsl_odeiv2_driver_apply (d, &t, ti, y); + + if (status != GSL_SUCCESS) { + printf ("error in ode, return value=%d\n", status); + break; + } + +// printf ("%.5e %.5e %.5e\n", t, y[0], y[1]); + for(j=0; j