summaryrefslogtreecommitdiff
path: root/lib/Numeric
diff options
context:
space:
mode:
authorAlberto Ruiz <aruiz@um.es>2012-02-25 12:36:30 +0100
committerAlberto Ruiz <aruiz@um.es>2012-02-25 12:36:30 +0100
commit47ea0f25eed9ab4a9bcdb4da2cf573a26883f9d4 (patch)
treeb91e101621e6a352bfeef68a1fe0a31df8684065 /lib/Numeric
parentb67148cf9e97a259c398280493829fc3cbed6582 (diff)
odeSolveV with safe method interface
Diffstat (limited to 'lib/Numeric')
-rw-r--r--lib/Numeric/GSL/ODE.hs47
-rw-r--r--lib/Numeric/GSL/gsl-aux.c40
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)@
29----------------------------------------------------------------------------- 29-----------------------------------------------------------------------------
30 30
31module Numeric.GSL.ODE ( 31module Numeric.GSL.ODE (
32 odeSolve, odeSolveV, ODEMethod(..) 32 odeSolve, odeSolveV, ODEMethod(..), Jacobian
33) where 33) where
34 34
35import Data.Packed.Internal 35import Data.Packed.Internal
@@ -41,18 +41,21 @@ import System.IO.Unsafe(unsafePerformIO)
41 41
42------------------------------------------------------------------------- 42-------------------------------------------------------------------------
43 43
44type Jacobian = Double -> Vector Double -> Matrix Double
45
44-- | Stepping functions 46-- | Stepping functions
45data ODEMethod = RK2 -- ^ Embedded Runge-Kutta (2, 3) method. 47data ODEMethod = RK2 -- ^ Embedded Runge-Kutta (2, 3) method.
46 | 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'. 48 | 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.
47 | RKf45 -- ^ Embedded Runge-Kutta-Fehlberg (4, 5) method. This method is a good general-purpose integrator. 49 | RKf45 -- ^ Embedded Runge-Kutta-Fehlberg (4, 5) method. This method is a good general-purpose integrator.
48 | RKck -- ^ Embedded Runge-Kutta Cash-Karp (4, 5) method. 50 | RKck -- ^ Embedded Runge-Kutta Cash-Karp (4, 5) method.
49 | RK8pd -- ^ Embedded Runge-Kutta Prince-Dormand (8,9) method. 51 | RK8pd -- ^ Embedded Runge-Kutta Prince-Dormand (8,9) method.
50 | RK2imp -- ^ Implicit 2nd order Runge-Kutta at Gaussian points. 52 | RK2imp Jacobian -- ^ Implicit 2nd order Runge-Kutta at Gaussian points.
51 | RK4imp -- ^ Implicit 4th order Runge-Kutta at Gaussian points. 53 | RK4imp Jacobian -- ^ Implicit 4th order Runge-Kutta at Gaussian points.
52 | BSimp -- ^ Implicit Bulirsch-Stoer method of Bader and Deuflhard. This algorithm requires the Jacobian. 54 | BSimp Jacobian -- ^ Implicit Bulirsch-Stoer method of Bader and Deuflhard. The method is generally suitable for stiff problems.
53 | MSAdams -- ^ A variable-coefficient linear multistep Adams method in Nordsieck form. 55 | 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.
54 | MSBDF -- ^ A variable-coefficient linear multistep backward differentiation formula (BDF) method in Nordsieck form. 56 | 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.
55 deriving (Enum,Eq,Show,Bounded) 57 | 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.
58
56 59
57-- | A version of 'odeSolveV' with reasonable default parameters and system of equations defined using lists. 60-- | A version of 'odeSolveV' with reasonable default parameters and system of equations defined using lists.
58odeSolve 61odeSolve
@@ -60,7 +63,7 @@ odeSolve
60 -> [Double] -- ^ initial conditions 63 -> [Double] -- ^ initial conditions
61 -> Vector Double -- ^ desired solution times 64 -> Vector Double -- ^ desired solution times
62 -> Matrix Double -- ^ solution 65 -> Matrix Double -- ^ solution
63odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) Nothing (fromList xi) ts 66odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) (fromList xi) ts
64 where hi = (ts@>1 - ts@>0)/100 67 where hi = (ts@>1 - ts@>0)/100
65 epsAbs = 1.49012e-08 68 epsAbs = 1.49012e-08
66 epsRel = 1.49012e-08 69 epsRel = 1.49012e-08
@@ -73,11 +76,33 @@ odeSolveV
73 -> Double -- ^ absolute tolerance for the state vector 76 -> Double -- ^ absolute tolerance for the state vector
74 -> Double -- ^ relative tolerance for the state vector 77 -> Double -- ^ relative tolerance for the state vector
75 -> (Double -> Vector Double -> Vector Double) -- ^ xdot(t,x) 78 -> (Double -> Vector Double -> Vector Double) -- ^ xdot(t,x)
79 -> Vector Double -- ^ initial conditions
80 -> Vector Double -- ^ desired solution times
81 -> Matrix Double -- ^ solution
82odeSolveV RK2 = odeSolveV' 0 Nothing
83odeSolveV RK4 = odeSolveV' 1 Nothing
84odeSolveV RKf45 = odeSolveV' 2 Nothing
85odeSolveV RKck = odeSolveV' 3 Nothing
86odeSolveV RK8pd = odeSolveV' 4 Nothing
87odeSolveV (RK2imp jac) = odeSolveV' 5 (Just jac)
88odeSolveV (RK4imp jac) = odeSolveV' 6 (Just jac)
89odeSolveV (BSimp jac) = odeSolveV' 7 (Just jac)
90odeSolveV (RK1imp jac) = odeSolveV' 8 (Just jac)
91odeSolveV MSAdams = odeSolveV' 9 Nothing
92odeSolveV (MSBDF jac) = odeSolveV' 10 (Just jac)
93
94
95odeSolveV'
96 :: CInt
76 -> Maybe (Double -> Vector Double -> Matrix Double) -- ^ optional jacobian 97 -> Maybe (Double -> Vector Double -> Matrix Double) -- ^ optional jacobian
98 -> Double -- ^ initial step size
99 -> Double -- ^ absolute tolerance for the state vector
100 -> Double -- ^ relative tolerance for the state vector
101 -> (Double -> Vector Double -> Vector Double) -- ^ xdot(t,x)
77 -> Vector Double -- ^ initial conditions 102 -> Vector Double -- ^ initial conditions
78 -> Vector Double -- ^ desired solution times 103 -> Vector Double -- ^ desired solution times
79 -> Matrix Double -- ^ solution 104 -> Matrix Double -- ^ solution
80odeSolveV method h epsAbs epsRel f mbjac xiv ts = unsafePerformIO $ do 105odeSolveV' method mbjac h epsAbs epsRel f xiv ts = unsafePerformIO $ do
81 let n = dim xiv 106 let n = dim xiv
82 fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t)) 107 fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t))
83 jp <- case mbjac of 108 jp <- case mbjac of
@@ -86,7 +111,7 @@ odeSolveV method h epsAbs epsRel f mbjac xiv ts = unsafePerformIO $ do
86 sol <- vec xiv $ \xiv' -> 111 sol <- vec xiv $ \xiv' ->
87 vec (checkTimes ts) $ \ts' -> 112 vec (checkTimes ts) $ \ts' ->
88 createMIO (dim ts) n 113 createMIO (dim ts) n
89 (ode_c (fi (fromEnum method)) h epsAbs epsRel fp jp // xiv' // ts' ) 114 (ode_c (method) h epsAbs epsRel fp jp // xiv' // ts' )
90 "ode" 115 "ode"
91 freeHaskellFunPtr fp 116 freeHaskellFunPtr fp
92 return sol 117 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,
1325 case 5 : {T = gsl_odeiv2_step_rk2imp; break; } 1325 case 5 : {T = gsl_odeiv2_step_rk2imp; break; }
1326 case 6 : {T = gsl_odeiv2_step_rk4imp; break; } 1326 case 6 : {T = gsl_odeiv2_step_rk4imp; break; }
1327 case 7 : {T = gsl_odeiv2_step_bsimp; break; } 1327 case 7 : {T = gsl_odeiv2_step_bsimp; break; }
1328 case 8 : {T = gsl_odeiv2_step_msadams; break; } 1328 case 8 : {T = gsl_odeiv2_step_rk1imp; break; }
1329 case 9 : {T = gsl_odeiv2_step_msbdf; break; } 1329 case 9 : {T = gsl_odeiv2_step_msadams; break; }
1330 case 10: {T = gsl_odeiv2_step_msbdf; break; }
1330 default: ERROR(BAD_CODE); 1331 default: ERROR(BAD_CODE);
1331 } 1332 }
1332 1333
1333
1334 gsl_odeiv2_step * s = gsl_odeiv2_step_alloc (T, xin);
1335 gsl_odeiv2_control * c = gsl_odeiv2_control_y_new (eps_abs, eps_rel);
1336 gsl_odeiv2_evolve * e = gsl_odeiv2_evolve_alloc (xin);
1337
1338 Tode P; 1334 Tode P;
1339 P.f = f; 1335 P.f = f;
1340 P.j = jac; 1336 P.j = jac;
@@ -1342,10 +1338,14 @@ int ode(int method, double h, double eps_abs, double eps_rel,
1342 1338
1343 gsl_odeiv2_system sys = {odefunc, odejac, xin, &P}; 1339 gsl_odeiv2_system sys = {odefunc, odejac, xin, &P};
1344 1340
1341 gsl_odeiv2_driver * d =
1342 gsl_odeiv2_driver_alloc_y_new (&sys, T, h, eps_abs, eps_rel);
1343
1345 double t = tsp[0]; 1344 double t = tsp[0];
1346 1345
1347 double* y = (double*)calloc(xin,sizeof(double)); 1346 double* y = (double*)calloc(xin,sizeof(double));
1348 int i,j; 1347 int i,j;
1348 int status;
1349 for(i=0; i< xin; i++) { 1349 for(i=0; i< xin; i++) {
1350 y[i] = xip[i]; 1350 y[i] = xip[i];
1351 solp[i] = xip[i]; 1351 solp[i] = xip[i];
@@ -1354,22 +1354,24 @@ int ode(int method, double h, double eps_abs, double eps_rel,
1354 for (i = 1; i < tsn ; i++) 1354 for (i = 1; i < tsn ; i++)
1355 { 1355 {
1356 double ti = tsp[i]; 1356 double ti = tsp[i];
1357 while (t < ti) 1357
1358 { 1358 status = gsl_odeiv2_driver_apply (d, &t, ti, y);
1359 gsl_odeiv2_evolve_apply (e, c, s, 1359
1360 &sys, 1360 if (status != GSL_SUCCESS) {
1361 &t, ti, &h, 1361 printf ("error in ode, return value=%d\n", status);
1362 y); 1362 break;
1363 // if (h < hmin) h = hmin; 1363 }
1364 } 1364
1365// printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
1366
1365 for(j=0; j<xin; j++) { 1367 for(j=0; j<xin; j++) {
1366 solp[i*xin + j] = y[j]; 1368 solp[i*xin + j] = y[j];
1367 } 1369 }
1368 } 1370 }
1369 1371
1370 free(y); 1372 free(y);
1371 gsl_odeiv2_evolve_free (e); 1373 gsl_odeiv2_driver_free (d);
1372 gsl_odeiv2_control_free (c); 1374
1373 gsl_odeiv2_step_free (s); 1375 return status;
1374 return 0;
1375} 1376}
1377