diff options
Diffstat (limited to 'packages/gsl/src/Numeric/GSL/ODE.hs')
-rw-r--r-- | packages/gsl/src/Numeric/GSL/ODE.hs | 136 |
1 files changed, 86 insertions, 50 deletions
diff --git a/packages/gsl/src/Numeric/GSL/ODE.hs b/packages/gsl/src/Numeric/GSL/ODE.hs index 7549a65..9e52873 100644 --- a/packages/gsl/src/Numeric/GSL/ODE.hs +++ b/packages/gsl/src/Numeric/GSL/ODE.hs | |||
@@ -1,3 +1,6 @@ | |||
1 | {-# LANGUAGE FlexibleContexts #-} | ||
2 | |||
3 | |||
1 | {- | | 4 | {- | |
2 | Module : Numeric.GSL.ODE | 5 | Module : Numeric.GSL.ODE |
3 | Copyright : (c) Alberto Ruiz 2010 | 6 | Copyright : (c) Alberto Ruiz 2010 |
@@ -29,10 +32,10 @@ main = mplot (ts : toColumns sol) | |||
29 | ----------------------------------------------------------------------------- | 32 | ----------------------------------------------------------------------------- |
30 | 33 | ||
31 | module Numeric.GSL.ODE ( | 34 | module Numeric.GSL.ODE ( |
32 | odeSolve, odeSolveV, ODEMethod(..), Jacobian | 35 | odeSolve, odeSolveV, odeSolveVWith, ODEMethod(..), Jacobian, StepControl(..) |
33 | ) where | 36 | ) where |
34 | 37 | ||
35 | import Data.Packed | 38 | import Numeric.LinearAlgebra.HMatrix |
36 | import Numeric.GSL.Internal | 39 | import Numeric.GSL.Internal |
37 | 40 | ||
38 | import Foreign.Ptr(FunPtr, nullFunPtr, freeHaskellFunPtr) | 41 | import Foreign.Ptr(FunPtr, nullFunPtr, freeHaskellFunPtr) |
@@ -41,9 +44,10 @@ import System.IO.Unsafe(unsafePerformIO) | |||
41 | 44 | ||
42 | ------------------------------------------------------------------------- | 45 | ------------------------------------------------------------------------- |
43 | 46 | ||
44 | type TVV = TV (TV Res) | 47 | type TVV = TV (TV Res) |
45 | type TVM = TV (TM Res) | 48 | type TVM = TV (TM Res) |
46 | type TVVM = TV (TV (TM Res)) | 49 | type TVVM = TV (TV (TM Res)) |
50 | type TVVVM = TV (TV (TV (TM Res))) | ||
47 | 51 | ||
48 | type Jacobian = Double -> Vector Double -> Matrix Double | 52 | type Jacobian = Double -> Vector Double -> Matrix Double |
49 | 53 | ||
@@ -60,73 +64,105 @@ data ODEMethod = RK2 -- ^ Embedded Runge-Kutta (2, 3) method. | |||
60 | | 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. |
61 | | 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. |
62 | 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 | ||
63 | 72 | ||
64 | -- | 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. |
65 | odeSolve | 74 | odeSolve |
66 | :: (Double -> [Double] -> [Double]) -- ^ xdot(t,x) | 75 | :: (Double -> [Double] -> [Double]) -- ^ x'(t,x) |
67 | -> [Double] -- ^ initial conditions | 76 | -> [Double] -- ^ initial conditions |
68 | -> Vector Double -- ^ desired solution times | 77 | -> Vector Double -- ^ desired solution times |
69 | -> Matrix Double -- ^ solution | 78 | -> Matrix Double -- ^ solution |
70 | 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 |
71 | where hi = (ts@>1 - ts@>0)/100 | 80 | where hi = (ts!1 - ts!0)/100 |
72 | epsAbs = 1.49012e-08 | 81 | epsAbs = 1.49012e-08 |
73 | epsRel = 1.49012e-08 | 82 | epsRel = epsAbs |
74 | l2v f = \t -> fromList . f t . toList | 83 | l2v f = \t -> fromList . f t . toList |
75 | 84 | ||
76 | -- | Evolution of the system with adaptive step-size control. | 85 | -- | A version of 'odeSolveVWith' with reasonable default step control. |
77 | odeSolveV | 86 | odeSolveV |
78 | :: ODEMethod | 87 | :: ODEMethod |
79 | -> Double -- ^ initial step size | 88 | -> Double -- ^ initial step size |
80 | -> Double -- ^ absolute tolerance for the state vector | 89 | -> Double -- ^ absolute tolerance for the state vector |
81 | -> Double -- ^ relative tolerance for the state vector | 90 | -> Double -- ^ relative tolerance for the state vector |
82 | -> (Double -> Vector Double -> Vector Double) -- ^ xdot(t,x) | 91 | -> (Double -> Vector Double -> Vector Double) -- ^ x'(t,x) |
83 | -> Vector Double -- ^ initial conditions | 92 | -> Vector Double -- ^ initial conditions |
84 | -> Vector Double -- ^ desired solution times | 93 | -> Vector Double -- ^ desired solution times |
85 | -> Matrix Double -- ^ solution | 94 | -> Matrix Double -- ^ solution |
86 | odeSolveV RK2 = odeSolveV' 0 Nothing | 95 | odeSolveV meth hi epsAbs epsRel = odeSolveVWith meth (XX' epsAbs epsRel 1 1) hi |
87 | odeSolveV RK4 = odeSolveV' 1 Nothing | 96 | |
88 | odeSolveV RKf45 = odeSolveV' 2 Nothing | 97 | -- | Evolution of the system with adaptive step-size control. |
89 | odeSolveV RKck = odeSolveV' 3 Nothing | 98 | odeSolveVWith |
90 | odeSolveV RK8pd = odeSolveV' 4 Nothing | 99 | :: ODEMethod |
91 | odeSolveV (RK2imp jac) = odeSolveV' 5 (Just jac) | 100 | -> StepControl |
92 | odeSolveV (RK4imp jac) = odeSolveV' 6 (Just jac) | 101 | -> Double -- ^ initial step size |
93 | odeSolveV (BSimp jac) = odeSolveV' 7 (Just jac) | 102 | -> (Double -> Vector Double -> Vector Double) -- ^ x'(t,x) |
94 | odeSolveV (RK1imp jac) = odeSolveV' 8 (Just jac) | ||
95 | odeSolveV MSAdams = odeSolveV' 9 Nothing | ||
96 | odeSolveV (MSBDF jac) = odeSolveV' 10 (Just jac) | ||
97 | |||
98 | |||
99 | odeSolveV' | ||
100 | :: CInt | ||
101 | -> Maybe (Double -> Vector Double -> Matrix Double) -- ^ optional jacobian | ||
102 | -> Double -- ^ initial step size | ||
103 | -> Double -- ^ absolute tolerance for the state vector | ||
104 | -> Double -- ^ relative tolerance for the state vector | ||
105 | -> (Double -> Vector Double -> Vector Double) -- ^ xdot(t,x) | ||
106 | -> Vector Double -- ^ initial conditions | 103 | -> Vector Double -- ^ initial conditions |
107 | -> Vector Double -- ^ desired solution times | 104 | -> Vector Double -- ^ desired solution times |
108 | -> Matrix Double -- ^ solution | 105 | -> Matrix Double -- ^ solution |
109 | odeSolveV' method mbjac h epsAbs epsRel f xiv ts = unsafePerformIO $ do | 106 | odeSolveVWith method control = odeSolveVWith' m mbj c epsAbs epsRel aX aX' mbsc |
110 | let n = dim xiv | 107 | where (m, mbj) = case method of |
111 | fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t)) | 108 | RK2 -> (0 , Nothing ) |
112 | jp <- case mbjac of | 109 | RK4 -> (1 , Nothing ) |
113 | Just jac -> mkDoubleVecMatfun (\t -> aux_vTom (checkdim2 n . jac t)) | 110 | RKf45 -> (2 , Nothing ) |
114 | Nothing -> return nullFunPtr | 111 | RKck -> (3 , Nothing ) |
115 | sol <- vec xiv $ \xiv' -> | 112 | RK8pd -> (4 , Nothing ) |
116 | vec (checkTimes ts) $ \ts' -> | 113 | RK2imp jac -> (5 , Just jac) |
117 | createMIO (dim ts) n | 114 | RK4imp jac -> (6 , Just jac) |
118 | (ode_c (method) h epsAbs epsRel fp jp // xiv' // ts' ) | 115 | BSimp jac -> (7 , Just jac) |
119 | "ode" | 116 | RK1imp jac -> (8 , Just jac) |
120 | freeHaskellFunPtr fp | 117 | MSAdams -> (9 , Nothing ) |
121 | 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 | ||
122 | 156 | ||
123 | foreign import ccall safe "ode" | 157 | foreign import ccall safe "ode" |
124 | 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 | ||
125 | 161 | ||
126 | ------------------------------------------------------- | 162 | ------------------------------------------------------- |
127 | 163 | ||
128 | checkdim1 n v | 164 | checkdim1 n v |
129 | | dim v == n = v | 165 | | size v == n = v |
130 | | otherwise = error $ "Error: "++ show n | 166 | | otherwise = error $ "Error: "++ show n |
131 | ++ " components expected in the result of the function supplied to odeSolve" | 167 | ++ " components expected in the result of the function supplied to odeSolve" |
132 | 168 | ||
@@ -135,6 +171,6 @@ checkdim2 n m | |||
135 | | otherwise = error $ "Error: "++ show n ++ "x" ++ show n | 171 | | otherwise = error $ "Error: "++ show n ++ "x" ++ show n |
136 | ++ " Jacobian expected in odeSolve" | 172 | ++ " Jacobian expected in odeSolve" |
137 | 173 | ||
138 | checkTimes ts | dim ts > 1 && all (>0) (zipWith subtract ts' (tail ts')) = ts | 174 | checkTimes ts | size ts > 1 && all (>0) (zipWith subtract ts' (tail ts')) = ts |
139 | | otherwise = error "odeSolve requires increasing times" | 175 | | otherwise = error "odeSolve requires increasing times" |
140 | where ts' = toList ts | 176 | where ts' = toList ts |