diff options
author | Alberto Ruiz <aruiz@um.es> | 2015-06-03 18:21:24 +0200 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2015-06-03 18:21:24 +0200 |
commit | b4fa1d9ca83f99f5057fcf70ed409b5038908f66 (patch) | |
tree | efedc612bfb9b27446fa8e043f314a26e95c9969 /packages/gsl | |
parent | 5481f94f5b936a2ad1947200838ea27ee860d0e7 (diff) | |
parent | 55e2326e76e7e1ee2ac0d753198b8f47e1c695ef (diff) |
Merge pull request #126 from ntfrgl/master
Generalize step control function
Diffstat (limited to 'packages/gsl')
-rw-r--r-- | packages/gsl/src/Numeric/GSL/ODE.hs | 126 | ||||
-rw-r--r-- | packages/gsl/src/Numeric/GSL/gsl-ode.c | 31 |
2 files changed, 103 insertions, 54 deletions
diff --git a/packages/gsl/src/Numeric/GSL/ODE.hs b/packages/gsl/src/Numeric/GSL/ODE.hs index 3258b83..9e52873 100644 --- a/packages/gsl/src/Numeric/GSL/ODE.hs +++ b/packages/gsl/src/Numeric/GSL/ODE.hs | |||
@@ -32,7 +32,7 @@ main = mplot (ts : toColumns sol) | |||
32 | ----------------------------------------------------------------------------- | 32 | ----------------------------------------------------------------------------- |
33 | 33 | ||
34 | module Numeric.GSL.ODE ( | 34 | module Numeric.GSL.ODE ( |
35 | odeSolve, odeSolveV, ODEMethod(..), Jacobian | 35 | odeSolve, odeSolveV, odeSolveVWith, ODEMethod(..), Jacobian, StepControl(..) |
36 | ) where | 36 | ) where |
37 | 37 | ||
38 | import Numeric.LinearAlgebra.HMatrix | 38 | import Numeric.LinearAlgebra.HMatrix |
@@ -44,9 +44,10 @@ import System.IO.Unsafe(unsafePerformIO) | |||
44 | 44 | ||
45 | ------------------------------------------------------------------------- | 45 | ------------------------------------------------------------------------- |
46 | 46 | ||
47 | type TVV = TV (TV Res) | 47 | type TVV = TV (TV Res) |
48 | type TVM = TV (TM Res) | 48 | type TVM = TV (TM Res) |
49 | type TVVM = TV (TV (TM Res)) | 49 | type TVVM = TV (TV (TM Res)) |
50 | type TVVVM = TV (TV (TV (TM Res))) | ||
50 | 51 | ||
51 | type Jacobian = Double -> Vector Double -> Matrix Double | 52 | type Jacobian = Double -> Vector Double -> Matrix Double |
52 | 53 | ||
@@ -63,68 +64,100 @@ data ODEMethod = RK2 -- ^ Embedded Runge-Kutta (2, 3) method. | |||
63 | | 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. | 64 | | 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. |
64 | | 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. | 65 | | 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. |
65 | 66 | ||
67 | -- | Adaptive step-size control functions | ||
68 | data StepControl = X Double Double -- ^ abs. and rel. tolerance for x(t) | ||
69 | | X' Double Double -- ^ abs. and rel. tolerance for x'(t) | ||
70 | | XX' Double Double Double Double -- ^ include both via rel. tolerance scaling factors a_x, a_x' | ||
71 | | ScXX' Double Double Double Double (Vector Double) -- ^ scale abs. tolerance of x(t) components | ||
66 | 72 | ||
67 | -- | A version of 'odeSolveV' with reasonable default parameters and system of equations defined using lists. | 73 | -- | A version of 'odeSolveV' with reasonable default parameters and system of equations defined using lists. |
68 | odeSolve | 74 | odeSolve |
69 | :: (Double -> [Double] -> [Double]) -- ^ xdot(t,x) | 75 | :: (Double -> [Double] -> [Double]) -- ^ x'(t,x) |
70 | -> [Double] -- ^ initial conditions | 76 | -> [Double] -- ^ initial conditions |
71 | -> Vector Double -- ^ desired solution times | 77 | -> Vector Double -- ^ desired solution times |
72 | -> Matrix Double -- ^ solution | 78 | -> Matrix Double -- ^ solution |
73 | odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) (fromList xi) ts | 79 | odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) (fromList xi) ts |
74 | where hi = (ts!1 - ts!0)/100 | 80 | where hi = (ts!1 - ts!0)/100 |
75 | epsAbs = 1.49012e-08 | 81 | epsAbs = 1.49012e-08 |
76 | epsRel = 1.49012e-08 | 82 | epsRel = epsAbs |
77 | l2v f = \t -> fromList . f t . toList | 83 | l2v f = \t -> fromList . f t . toList |
78 | 84 | ||
79 | -- | Evolution of the system with adaptive step-size control. | 85 | -- | A version of 'odeSolveVWith' with reasonable default step control. |
80 | odeSolveV | 86 | odeSolveV |
81 | :: ODEMethod | 87 | :: ODEMethod |
82 | -> Double -- ^ initial step size | 88 | -> Double -- ^ initial step size |
83 | -> Double -- ^ absolute tolerance for the state vector | 89 | -> Double -- ^ absolute tolerance for the state vector |
84 | -> Double -- ^ relative tolerance for the state vector | 90 | -> Double -- ^ relative tolerance for the state vector |
85 | -> (Double -> Vector Double -> Vector Double) -- ^ xdot(t,x) | 91 | -> (Double -> Vector Double -> Vector Double) -- ^ x'(t,x) |
86 | -> Vector Double -- ^ initial conditions | 92 | -> Vector Double -- ^ initial conditions |
87 | -> Vector Double -- ^ desired solution times | 93 | -> Vector Double -- ^ desired solution times |
88 | -> Matrix Double -- ^ solution | 94 | -> Matrix Double -- ^ solution |
89 | odeSolveV RK2 = odeSolveV' 0 Nothing | 95 | odeSolveV meth hi epsAbs epsRel = odeSolveVWith meth (XX' epsAbs epsRel 1 1) hi |
90 | odeSolveV RK4 = odeSolveV' 1 Nothing | 96 | |
91 | odeSolveV RKf45 = odeSolveV' 2 Nothing | 97 | -- | Evolution of the system with adaptive step-size control. |
92 | odeSolveV RKck = odeSolveV' 3 Nothing | 98 | odeSolveVWith |
93 | odeSolveV RK8pd = odeSolveV' 4 Nothing | 99 | :: ODEMethod |
94 | odeSolveV (RK2imp jac) = odeSolveV' 5 (Just jac) | 100 | -> StepControl |
95 | odeSolveV (RK4imp jac) = odeSolveV' 6 (Just jac) | 101 | -> Double -- ^ initial step size |
96 | odeSolveV (BSimp jac) = odeSolveV' 7 (Just jac) | 102 | -> (Double -> Vector Double -> Vector Double) -- ^ x'(t,x) |
97 | odeSolveV (RK1imp jac) = odeSolveV' 8 (Just jac) | ||
98 | odeSolveV MSAdams = odeSolveV' 9 Nothing | ||
99 | odeSolveV (MSBDF jac) = odeSolveV' 10 (Just jac) | ||
100 | |||
101 | |||
102 | odeSolveV' | ||
103 | :: CInt | ||
104 | -> Maybe (Double -> Vector Double -> Matrix Double) -- ^ optional jacobian | ||
105 | -> Double -- ^ initial step size | ||
106 | -> Double -- ^ absolute tolerance for the state vector | ||
107 | -> Double -- ^ relative tolerance for the state vector | ||
108 | -> (Double -> Vector Double -> Vector Double) -- ^ xdot(t,x) | ||
109 | -> Vector Double -- ^ initial conditions | 103 | -> Vector Double -- ^ initial conditions |
110 | -> Vector Double -- ^ desired solution times | 104 | -> Vector Double -- ^ desired solution times |
111 | -> Matrix Double -- ^ solution | 105 | -> Matrix Double -- ^ solution |
112 | odeSolveV' method mbjac h epsAbs epsRel f xiv ts = unsafePerformIO $ do | 106 | odeSolveVWith method control = odeSolveVWith' m mbj c epsAbs epsRel aX aX' mbsc |
113 | let n = size xiv | 107 | where (m, mbj) = case method of |
114 | fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t)) | 108 | RK2 -> (0 , Nothing ) |
115 | jp <- case mbjac of | 109 | RK4 -> (1 , Nothing ) |
116 | Just jac -> mkDoubleVecMatfun (\t -> aux_vTom (checkdim2 n . jac t)) | 110 | RKf45 -> (2 , Nothing ) |
117 | Nothing -> return nullFunPtr | 111 | RKck -> (3 , Nothing ) |
118 | sol <- vec xiv $ \xiv' -> | 112 | RK8pd -> (4 , Nothing ) |
119 | vec (checkTimes ts) $ \ts' -> | 113 | RK2imp jac -> (5 , Just jac) |
120 | createMIO (size ts) n | 114 | RK4imp jac -> (6 , Just jac) |
121 | (ode_c (method) h epsAbs epsRel fp jp // xiv' // ts' ) | 115 | BSimp jac -> (7 , Just jac) |
122 | "ode" | 116 | RK1imp jac -> (8 , Just jac) |
123 | freeHaskellFunPtr fp | 117 | MSAdams -> (9 , Nothing ) |
124 | return sol | 118 | MSBDF jac -> (10, Just jac) |
119 | (c, epsAbs, epsRel, aX, aX', mbsc) = case control of | ||
120 | X ea er -> (0, ea, er, 1 , 0 , Nothing) | ||
121 | X' ea er -> (0, ea, er, 0 , 1 , Nothing) | ||
122 | XX' ea er ax ax' -> (0, ea, er, ax, ax', Nothing) | ||
123 | ScXX' ea er ax ax' sc -> (1, ea, er, ax, ax', Just sc) | ||
124 | |||
125 | odeSolveVWith' | ||
126 | :: CInt -- ^ stepping function | ||
127 | -> Maybe (Double -> Vector Double -> Matrix Double) -- ^ optional jacobian | ||
128 | -> CInt -- ^ step-size control function | ||
129 | -> Double -- ^ absolute tolerance for step-size control | ||
130 | -> Double -- ^ relative tolerance for step-size control | ||
131 | -> Double -- ^ scaling factor for relative tolerance of x(t) | ||
132 | -> Double -- ^ scaling factor for relative tolerance of x'(t) | ||
133 | -> Maybe (Vector Double) -- ^ optional scaling for absolute error | ||
134 | -> Double -- ^ initial step size | ||
135 | -> (Double -> Vector Double -> Vector Double) -- ^ x'(t,x) | ||
136 | -> Vector Double -- ^ initial conditions | ||
137 | -> Vector Double -- ^ desired solution times | ||
138 | -> Matrix Double -- ^ solution | ||
139 | odeSolveVWith' method mbjac control epsAbs epsRel aX aX' mbsc h f xiv ts = | ||
140 | unsafePerformIO $ do | ||
141 | let n = size xiv | ||
142 | sc = case mbsc of | ||
143 | Just scv -> checkdim1 n scv | ||
144 | Nothing -> xiv | ||
145 | fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t)) | ||
146 | jp <- case mbjac of | ||
147 | Just jac -> mkDoubleVecMatfun (\t -> aux_vTom (checkdim2 n . jac t)) | ||
148 | Nothing -> return nullFunPtr | ||
149 | sol <- vec sc $ \sc' -> vec xiv $ \xiv' -> | ||
150 | vec (checkTimes ts) $ \ts' -> createMIO (size ts) n | ||
151 | (ode_c method control h epsAbs epsRel aX aX' fp jp | ||
152 | // sc' // xiv' // ts' ) | ||
153 | "ode" | ||
154 | freeHaskellFunPtr fp | ||
155 | return sol | ||
125 | 156 | ||
126 | foreign import ccall safe "ode" | 157 | foreign import ccall safe "ode" |
127 | ode_c :: CInt -> Double -> Double -> Double -> FunPtr (Double -> TVV) -> FunPtr (Double -> TVM) -> TVVM | 158 | ode_c :: CInt -> CInt -> Double |
159 | -> Double -> Double -> Double -> Double | ||
160 | -> FunPtr (Double -> TVV) -> FunPtr (Double -> TVM) -> TVVVM | ||
128 | 161 | ||
129 | ------------------------------------------------------- | 162 | ------------------------------------------------------- |
130 | 163 | ||
@@ -141,4 +174,3 @@ checkdim2 n m | |||
141 | checkTimes ts | size ts > 1 && all (>0) (zipWith subtract ts' (tail ts')) = ts | 174 | checkTimes ts | size ts > 1 && all (>0) (zipWith subtract ts' (tail ts')) = ts |
142 | | otherwise = error "odeSolve requires increasing times" | 175 | | otherwise = error "odeSolve requires increasing times" |
143 | where ts' = toList ts | 176 | where ts' = toList ts |
144 | |||
diff --git a/packages/gsl/src/Numeric/GSL/gsl-ode.c b/packages/gsl/src/Numeric/GSL/gsl-ode.c index 3f2771b..a6bdb55 100644 --- a/packages/gsl/src/Numeric/GSL/gsl-ode.c +++ b/packages/gsl/src/Numeric/GSL/gsl-ode.c | |||
@@ -23,10 +23,11 @@ int odejac (double t, const double y[], double *dfdy, double dfdt[], void *param | |||
23 | } | 23 | } |
24 | 24 | ||
25 | 25 | ||
26 | int ode(int method, double h, double eps_abs, double eps_rel, | 26 | int ode(int method, int control, double h, |
27 | double eps_abs, double eps_rel, double a_y, double a_dydt, | ||
27 | int f(double, int, const double*, int, double*), | 28 | int f(double, int, const double*, int, double*), |
28 | int jac(double, int, const double*, int, int, double*), | 29 | int jac(double, int, const double*, int, int, double*), |
29 | KRVEC(xi), KRVEC(ts), RMAT(sol)) { | 30 | KRVEC(sc), KRVEC(xi), KRVEC(ts), RMAT(sol)) { |
30 | 31 | ||
31 | const gsl_odeiv_step_type * T; | 32 | const gsl_odeiv_step_type * T; |
32 | 33 | ||
@@ -46,8 +47,16 @@ int ode(int method, double h, double eps_abs, double eps_rel, | |||
46 | } | 47 | } |
47 | 48 | ||
48 | gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, xin); | 49 | gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, xin); |
49 | gsl_odeiv_control * c = gsl_odeiv_control_y_new (eps_abs, eps_rel); | ||
50 | gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (xin); | 50 | gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (xin); |
51 | gsl_odeiv_control * c; | ||
52 | |||
53 | switch(control) { | ||
54 | case 0: { c = gsl_odeiv_control_standard_new | ||
55 | (eps_abs, eps_rel, a_y, a_dydt); break; } | ||
56 | case 1: { c = gsl_odeiv_control_scaled_new | ||
57 | (eps_abs, eps_rel, a_y, a_dydt, scp, scn); break; } | ||
58 | default: ERROR(BAD_CODE); | ||
59 | } | ||
51 | 60 | ||
52 | Tode P; | 61 | Tode P; |
53 | P.f = f; | 62 | P.f = f; |
@@ -112,10 +121,11 @@ int odejac (double t, const double y[], double *dfdy, double dfdt[], void *param | |||
112 | } | 121 | } |
113 | 122 | ||
114 | 123 | ||
115 | int ode(int method, double h, double eps_abs, double eps_rel, | 124 | int ode(int method, int control, double h, |
125 | double eps_abs, double eps_rel, double a_y, double a_dydt, | ||
116 | int f(double, int, const double*, int, double*), | 126 | int f(double, int, const double*, int, double*), |
117 | int jac(double, int, const double*, int, int, double*), | 127 | int jac(double, int, const double*, int, int, double*), |
118 | KRVEC(xi), KRVEC(ts), RMAT(sol)) { | 128 | KRVEC(sc), KRVEC(xi), KRVEC(ts), RMAT(sol)) { |
119 | 129 | ||
120 | const gsl_odeiv2_step_type * T; | 130 | const gsl_odeiv2_step_type * T; |
121 | 131 | ||
@@ -141,8 +151,15 @@ int ode(int method, double h, double eps_abs, double eps_rel, | |||
141 | 151 | ||
142 | gsl_odeiv2_system sys = {odefunc, odejac, xin, &P}; | 152 | gsl_odeiv2_system sys = {odefunc, odejac, xin, &P}; |
143 | 153 | ||
144 | gsl_odeiv2_driver * d = | 154 | gsl_odeiv2_driver * d; |
145 | gsl_odeiv2_driver_alloc_y_new (&sys, T, h, eps_abs, eps_rel); | 155 | |
156 | switch(control) { | ||
157 | case 0: { d = gsl_odeiv2_driver_alloc_standard_new | ||
158 | (&sys, T, h, eps_abs, eps_rel, a_y, a_dydt); break; } | ||
159 | case 1: { d = gsl_odeiv2_driver_alloc_scaled_new | ||
160 | (&sys, T, h, eps_abs, eps_rel, a_y, a_dydt, scp); break; } | ||
161 | default: ERROR(BAD_CODE); | ||
162 | } | ||
146 | 163 | ||
147 | double t = tsp[0]; | 164 | double t = tsp[0]; |
148 | 165 | ||