summaryrefslogtreecommitdiff
path: root/packages/hmatrix/src/Numeric/GSL/ODE.hs
diff options
context:
space:
mode:
authorAlberto Ruiz <aruiz@um.es>2014-05-21 10:30:55 +0200
committerAlberto Ruiz <aruiz@um.es>2014-05-21 10:30:55 +0200
commit197e88c3b56d28840217010a2871c6ea3a4dd1a4 (patch)
tree825be9d6c9d87d23f7e5497c0133d11d52c63535 /packages/hmatrix/src/Numeric/GSL/ODE.hs
parente07c3dee7235496b71a89233106d93f6cc94ada1 (diff)
update dependencies, move examples etc
Diffstat (limited to 'packages/hmatrix/src/Numeric/GSL/ODE.hs')
-rw-r--r--packages/hmatrix/src/Numeric/GSL/ODE.hs140
1 files changed, 0 insertions, 140 deletions
diff --git a/packages/hmatrix/src/Numeric/GSL/ODE.hs b/packages/hmatrix/src/Numeric/GSL/ODE.hs
deleted file mode 100644
index 7549a65..0000000
--- a/packages/hmatrix/src/Numeric/GSL/ODE.hs
+++ /dev/null
@@ -1,140 +0,0 @@
1{- |
2Module : Numeric.GSL.ODE
3Copyright : (c) Alberto Ruiz 2010
4License : GPL
5Maintainer : Alberto Ruiz
6Stability : provisional
7
8Solution of ordinary differential equation (ODE) initial value problems.
9
10<http://www.gnu.org/software/gsl/manual/html_node/Ordinary-Differential-Equations.html>
11
12A simple example:
13
14@
15import Numeric.GSL.ODE
16import Numeric.LinearAlgebra
17import Numeric.LinearAlgebra.Util(mplot)
18
19xdot t [x,v] = [v, -0.95*x - 0.1*v]
20
21ts = linspace 100 (0,20 :: Double)
22
23sol = odeSolve xdot [10,0] ts
24
25main = mplot (ts : toColumns sol)
26@
27
28-}
29-----------------------------------------------------------------------------
30
31module Numeric.GSL.ODE (
32 odeSolve, odeSolveV, ODEMethod(..), Jacobian
33) where
34
35import Data.Packed
36import Numeric.GSL.Internal
37
38import Foreign.Ptr(FunPtr, nullFunPtr, freeHaskellFunPtr)
39import Foreign.C.Types
40import System.IO.Unsafe(unsafePerformIO)
41
42-------------------------------------------------------------------------
43
44type TVV = TV (TV Res)
45type TVM = TV (TM Res)
46type TVVM = TV (TV (TM Res))
47
48type Jacobian = Double -> Vector Double -> Matrix Double
49
50-- | Stepping functions
51data ODEMethod = RK2 -- ^ Embedded Runge-Kutta (2, 3) method.
52 | 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.
53 | RKf45 -- ^ Embedded Runge-Kutta-Fehlberg (4, 5) method. This method is a good general-purpose integrator.
54 | RKck -- ^ Embedded Runge-Kutta Cash-Karp (4, 5) method.
55 | RK8pd -- ^ Embedded Runge-Kutta Prince-Dormand (8,9) method.
56 | RK2imp Jacobian -- ^ Implicit 2nd order Runge-Kutta at Gaussian points.
57 | RK4imp Jacobian -- ^ Implicit 4th order Runge-Kutta at Gaussian points.
58 | BSimp Jacobian -- ^ Implicit Bulirsch-Stoer method of Bader and Deuflhard. The method is generally suitable for stiff problems.
59 | 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.
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.
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.
62
63
64-- | A version of 'odeSolveV' with reasonable default parameters and system of equations defined using lists.
65odeSolve
66 :: (Double -> [Double] -> [Double]) -- ^ xdot(t,x)
67 -> [Double] -- ^ initial conditions
68 -> Vector Double -- ^ desired solution times
69 -> Matrix Double -- ^ solution
70odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) (fromList xi) ts
71 where hi = (ts@>1 - ts@>0)/100
72 epsAbs = 1.49012e-08
73 epsRel = 1.49012e-08
74 l2v f = \t -> fromList . f t . toList
75
76-- | Evolution of the system with adaptive step-size control.
77odeSolveV
78 :: ODEMethod
79 -> Double -- ^ initial step size
80 -> Double -- ^ absolute tolerance for the state vector
81 -> Double -- ^ relative tolerance for the state vector
82 -> (Double -> Vector Double -> Vector Double) -- ^ xdot(t,x)
83 -> Vector Double -- ^ initial conditions
84 -> Vector Double -- ^ desired solution times
85 -> Matrix Double -- ^ solution
86odeSolveV RK2 = odeSolveV' 0 Nothing
87odeSolveV RK4 = odeSolveV' 1 Nothing
88odeSolveV RKf45 = odeSolveV' 2 Nothing
89odeSolveV RKck = odeSolveV' 3 Nothing
90odeSolveV RK8pd = odeSolveV' 4 Nothing
91odeSolveV (RK2imp jac) = odeSolveV' 5 (Just jac)
92odeSolveV (RK4imp jac) = odeSolveV' 6 (Just jac)
93odeSolveV (BSimp jac) = odeSolveV' 7 (Just jac)
94odeSolveV (RK1imp jac) = odeSolveV' 8 (Just jac)
95odeSolveV MSAdams = odeSolveV' 9 Nothing
96odeSolveV (MSBDF jac) = odeSolveV' 10 (Just jac)
97
98
99odeSolveV'
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
107 -> Vector Double -- ^ desired solution times
108 -> Matrix Double -- ^ solution
109odeSolveV' method mbjac h epsAbs epsRel f xiv ts = unsafePerformIO $ do
110 let n = dim xiv
111 fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t))
112 jp <- case mbjac of
113 Just jac -> mkDoubleVecMatfun (\t -> aux_vTom (checkdim2 n . jac t))
114 Nothing -> return nullFunPtr
115 sol <- vec xiv $ \xiv' ->
116 vec (checkTimes ts) $ \ts' ->
117 createMIO (dim ts) n
118 (ode_c (method) h epsAbs epsRel fp jp // xiv' // ts' )
119 "ode"
120 freeHaskellFunPtr fp
121 return sol
122
123foreign import ccall safe "ode"
124 ode_c :: CInt -> Double -> Double -> Double -> FunPtr (Double -> TVV) -> FunPtr (Double -> TVM) -> TVVM
125
126-------------------------------------------------------
127
128checkdim1 n v
129 | dim v == n = v
130 | otherwise = error $ "Error: "++ show n
131 ++ " components expected in the result of the function supplied to odeSolve"
132
133checkdim2 n m
134 | rows m == n && cols m == n = m
135 | otherwise = error $ "Error: "++ show n ++ "x" ++ show n
136 ++ " Jacobian expected in odeSolve"
137
138checkTimes ts | dim ts > 1 && all (>0) (zipWith subtract ts' (tail ts')) = ts
139 | otherwise = error "odeSolve requires increasing times"
140 where ts' = toList ts