diff options
author | Alberto Ruiz <aruiz@um.es> | 2012-02-25 12:36:30 +0100 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2012-02-25 12:36:30 +0100 |
commit | 47ea0f25eed9ab4a9bcdb4da2cf573a26883f9d4 (patch) | |
tree | b91e101621e6a352bfeef68a1fe0a31df8684065 /lib | |
parent | b67148cf9e97a259c398280493829fc3cbed6582 (diff) |
odeSolveV with safe method interface
Diffstat (limited to 'lib')
-rw-r--r-- | lib/Numeric/GSL/ODE.hs | 47 | ||||
-rw-r--r-- | lib/Numeric/GSL/gsl-aux.c | 40 |
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 | ||
31 | module Numeric.GSL.ODE ( | 31 | module Numeric.GSL.ODE ( |
32 | odeSolve, odeSolveV, ODEMethod(..) | 32 | odeSolve, odeSolveV, ODEMethod(..), Jacobian |
33 | ) where | 33 | ) where |
34 | 34 | ||
35 | import Data.Packed.Internal | 35 | import Data.Packed.Internal |
@@ -41,18 +41,21 @@ import System.IO.Unsafe(unsafePerformIO) | |||
41 | 41 | ||
42 | ------------------------------------------------------------------------- | 42 | ------------------------------------------------------------------------- |
43 | 43 | ||
44 | type Jacobian = Double -> Vector Double -> Matrix Double | ||
45 | |||
44 | -- | Stepping functions | 46 | -- | Stepping functions |
45 | data ODEMethod = RK2 -- ^ Embedded Runge-Kutta (2, 3) method. | 47 | data 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. |
58 | odeSolve | 61 | odeSolve |
@@ -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 |
63 | odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) Nothing (fromList xi) ts | 66 | odeSolve 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 | ||
82 | odeSolveV RK2 = odeSolveV' 0 Nothing | ||
83 | odeSolveV RK4 = odeSolveV' 1 Nothing | ||
84 | odeSolveV RKf45 = odeSolveV' 2 Nothing | ||
85 | odeSolveV RKck = odeSolveV' 3 Nothing | ||
86 | odeSolveV RK8pd = odeSolveV' 4 Nothing | ||
87 | odeSolveV (RK2imp jac) = odeSolveV' 5 (Just jac) | ||
88 | odeSolveV (RK4imp jac) = odeSolveV' 6 (Just jac) | ||
89 | odeSolveV (BSimp jac) = odeSolveV' 7 (Just jac) | ||
90 | odeSolveV (RK1imp jac) = odeSolveV' 8 (Just jac) | ||
91 | odeSolveV MSAdams = odeSolveV' 9 Nothing | ||
92 | odeSolveV (MSBDF jac) = odeSolveV' 10 (Just jac) | ||
93 | |||
94 | |||
95 | odeSolveV' | ||
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 |
80 | odeSolveV method h epsAbs epsRel f mbjac xiv ts = unsafePerformIO $ do | 105 | odeSolveV' 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 | |||