summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--packages/gsl/src/Numeric/GSL/ODE.hs126
-rw-r--r--packages/gsl/src/Numeric/GSL/gsl-ode.c31
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
34module Numeric.GSL.ODE ( 34module Numeric.GSL.ODE (
35 odeSolve, odeSolveV, ODEMethod(..), Jacobian 35 odeSolve, odeSolveV, odeSolveVWith, ODEMethod(..), Jacobian, StepControl(..)
36) where 36) where
37 37
38import Numeric.LinearAlgebra.HMatrix 38import Numeric.LinearAlgebra.HMatrix
@@ -44,9 +44,10 @@ import System.IO.Unsafe(unsafePerformIO)
44 44
45------------------------------------------------------------------------- 45-------------------------------------------------------------------------
46 46
47type TVV = TV (TV Res) 47type TVV = TV (TV Res)
48type TVM = TV (TM Res) 48type TVM = TV (TM Res)
49type TVVM = TV (TV (TM Res)) 49type TVVM = TV (TV (TM Res))
50type TVVVM = TV (TV (TV (TM Res)))
50 51
51type Jacobian = Double -> Vector Double -> Matrix Double 52type 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
68data 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.
68odeSolve 74odeSolve
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
73odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) (fromList xi) ts 79odeSolve 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.
80odeSolveV 86odeSolveV
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
89odeSolveV RK2 = odeSolveV' 0 Nothing 95odeSolveV meth hi epsAbs epsRel = odeSolveVWith meth (XX' epsAbs epsRel 1 1) hi
90odeSolveV RK4 = odeSolveV' 1 Nothing 96
91odeSolveV RKf45 = odeSolveV' 2 Nothing 97-- | Evolution of the system with adaptive step-size control.
92odeSolveV RKck = odeSolveV' 3 Nothing 98odeSolveVWith
93odeSolveV RK8pd = odeSolveV' 4 Nothing 99 :: ODEMethod
94odeSolveV (RK2imp jac) = odeSolveV' 5 (Just jac) 100 -> StepControl
95odeSolveV (RK4imp jac) = odeSolveV' 6 (Just jac) 101 -> Double -- ^ initial step size
96odeSolveV (BSimp jac) = odeSolveV' 7 (Just jac) 102 -> (Double -> Vector Double -> Vector Double) -- ^ x'(t,x)
97odeSolveV (RK1imp jac) = odeSolveV' 8 (Just jac)
98odeSolveV MSAdams = odeSolveV' 9 Nothing
99odeSolveV (MSBDF jac) = odeSolveV' 10 (Just jac)
100
101
102odeSolveV'
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
112odeSolveV' method mbjac h epsAbs epsRel f xiv ts = unsafePerformIO $ do 106odeSolveVWith 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
125odeSolveVWith'
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
139odeSolveVWith' 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
126foreign import ccall safe "ode" 157foreign 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
141checkTimes ts | size ts > 1 && all (>0) (zipWith subtract ts' (tail ts')) = ts 174checkTimes 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
26int ode(int method, double h, double eps_abs, double eps_rel, 26int 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
115int ode(int method, double h, double eps_abs, double eps_rel, 124int 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