odeSolveV with safe method interface

2 files changed, 57 insertions, 30 deletions

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<xin; j++) {
                solp[i*xin + j] = y[j];
            }
          }
 
     free(y);
-    gsl_odeiv2_evolve_free (e);
-    gsl_odeiv2_control_free (c);
-    gsl_odeiv2_step_free (s);
-    return 0;
+    gsl_odeiv2_driver_free (d);
+    
+    return status;
 }
+