diff options
Diffstat (limited to 'packages/gsl/src/Numeric/GSL/ODE.hs')
-rw-r--r-- | packages/gsl/src/Numeric/GSL/ODE.hs | 16 |
1 files changed, 10 insertions, 6 deletions
diff --git a/packages/gsl/src/Numeric/GSL/ODE.hs b/packages/gsl/src/Numeric/GSL/ODE.hs index 7549a65..3258b83 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 |
@@ -32,7 +35,7 @@ module Numeric.GSL.ODE ( | |||
32 | odeSolve, odeSolveV, ODEMethod(..), Jacobian | 35 | odeSolve, odeSolveV, ODEMethod(..), Jacobian |
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) |
@@ -68,7 +71,7 @@ odeSolve | |||
68 | -> Vector Double -- ^ desired solution times | 71 | -> Vector Double -- ^ desired solution times |
69 | -> Matrix Double -- ^ solution | 72 | -> Matrix Double -- ^ solution |
70 | odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) (fromList xi) ts | 73 | odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) (fromList xi) ts |
71 | where hi = (ts@>1 - ts@>0)/100 | 74 | where hi = (ts!1 - ts!0)/100 |
72 | epsAbs = 1.49012e-08 | 75 | epsAbs = 1.49012e-08 |
73 | epsRel = 1.49012e-08 | 76 | epsRel = 1.49012e-08 |
74 | l2v f = \t -> fromList . f t . toList | 77 | l2v f = \t -> fromList . f t . toList |
@@ -107,14 +110,14 @@ odeSolveV' | |||
107 | -> Vector Double -- ^ desired solution times | 110 | -> Vector Double -- ^ desired solution times |
108 | -> Matrix Double -- ^ solution | 111 | -> Matrix Double -- ^ solution |
109 | odeSolveV' method mbjac h epsAbs epsRel f xiv ts = unsafePerformIO $ do | 112 | odeSolveV' method mbjac h epsAbs epsRel f xiv ts = unsafePerformIO $ do |
110 | let n = dim xiv | 113 | let n = size xiv |
111 | fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t)) | 114 | fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t)) |
112 | jp <- case mbjac of | 115 | jp <- case mbjac of |
113 | Just jac -> mkDoubleVecMatfun (\t -> aux_vTom (checkdim2 n . jac t)) | 116 | Just jac -> mkDoubleVecMatfun (\t -> aux_vTom (checkdim2 n . jac t)) |
114 | Nothing -> return nullFunPtr | 117 | Nothing -> return nullFunPtr |
115 | sol <- vec xiv $ \xiv' -> | 118 | sol <- vec xiv $ \xiv' -> |
116 | vec (checkTimes ts) $ \ts' -> | 119 | vec (checkTimes ts) $ \ts' -> |
117 | createMIO (dim ts) n | 120 | createMIO (size ts) n |
118 | (ode_c (method) h epsAbs epsRel fp jp // xiv' // ts' ) | 121 | (ode_c (method) h epsAbs epsRel fp jp // xiv' // ts' ) |
119 | "ode" | 122 | "ode" |
120 | freeHaskellFunPtr fp | 123 | freeHaskellFunPtr fp |
@@ -126,7 +129,7 @@ foreign import ccall safe "ode" | |||
126 | ------------------------------------------------------- | 129 | ------------------------------------------------------- |
127 | 130 | ||
128 | checkdim1 n v | 131 | checkdim1 n v |
129 | | dim v == n = v | 132 | | size v == n = v |
130 | | otherwise = error $ "Error: "++ show n | 133 | | otherwise = error $ "Error: "++ show n |
131 | ++ " components expected in the result of the function supplied to odeSolve" | 134 | ++ " components expected in the result of the function supplied to odeSolve" |
132 | 135 | ||
@@ -135,6 +138,7 @@ checkdim2 n m | |||
135 | | otherwise = error $ "Error: "++ show n ++ "x" ++ show n | 138 | | otherwise = error $ "Error: "++ show n ++ "x" ++ show n |
136 | ++ " Jacobian expected in odeSolve" | 139 | ++ " Jacobian expected in odeSolve" |
137 | 140 | ||
138 | checkTimes ts | dim ts > 1 && all (>0) (zipWith subtract ts' (tail ts')) = ts | 141 | checkTimes ts | size ts > 1 && all (>0) (zipWith subtract ts' (tail ts')) = ts |
139 | | otherwise = error "odeSolve requires increasing times" | 142 | | otherwise = error "odeSolve requires increasing times" |
140 | where ts' = toList ts | 143 | where ts' = toList ts |
144 | |||