From 47ea0f25eed9ab4a9bcdb4da2cf573a26883f9d4 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Sat, 25 Feb 2012 12:36:30 +0100 Subject: odeSolveV with safe method interface --- lib/Numeric/GSL/ODE.hs | 47 ++++++++++++++++++++++++++++++++++++----------- lib/Numeric/GSL/gsl-aux.c | 40 +++++++++++++++++++++------------------- 2 files changed, 57 insertions(+), 30 deletions(-) (limited to 'lib') diff --git a/lib/Numeric/GSL/ODE.hs b/lib/Numeric/GSL/ODE.hs index d4f83aa..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. - | MSAdams -- ^ A variable-coefficient linear multistep Adams method in Nordsieck form. - | MSBDF -- ^ A variable-coefficient linear multistep backward differentiation formula (BDF) method in Nordsieck form. - 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 bf3f684..f74e2e3 100644 --- a/lib/Numeric/GSL/gsl-aux.c +++ b/lib/Numeric/GSL/gsl-aux.c @@ -1325,16 +1325,12 @@ int ode(int method, double h, double eps_abs, double eps_rel, 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_msadams; break; } - case 9 : {T = gsl_odeiv2_step_msbdf; 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_odeiv2_step * s = gsl_odeiv2_step_alloc (T, xin); - gsl_odeiv2_control * c = gsl_odeiv2_control_y_new (eps_abs, eps_rel); - gsl_odeiv2_evolve * e = gsl_odeiv2_evolve_alloc (xin); - Tode P; P.f = f; P.j = jac; @@ -1342,10 +1338,14 @@ int ode(int method, double h, double eps_abs, double eps_rel, 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]; @@ -1354,22 +1354,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_odeiv2_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