summaryrefslogtreecommitdiff
path: root/packages/sundials/src/Main.hs
diff options
context:
space:
mode:
Diffstat (limited to 'packages/sundials/src/Main.hs')
-rw-r--r--packages/sundials/src/Main.hs14
1 files changed, 8 insertions, 6 deletions
diff --git a/packages/sundials/src/Main.hs b/packages/sundials/src/Main.hs
index 01d3595..4e0cf4b 100644
--- a/packages/sundials/src/Main.hs
+++ b/packages/sundials/src/Main.hs
@@ -1,6 +1,5 @@
1{-# OPTIONS_GHC -Wall #-} 1{-# OPTIONS_GHC -Wall #-}
2 2
3import qualified Data.Vector.Storable as V
4import Numeric.Sundials.Arkode.ODE 3import Numeric.Sundials.Arkode.ODE
5import Numeric.LinearAlgebra 4import Numeric.LinearAlgebra
6 5
@@ -9,14 +8,14 @@ import qualified Diagrams.Prelude as D
9import Diagrams.Backend.Rasterific 8import Diagrams.Backend.Rasterific
10 9
11import Control.Lens 10import Control.Lens
12import Data.List (zip4) 11import Data.List (intercalate)
13 12
14import Text.PrettyPrint.HughesPJClass 13import Text.PrettyPrint.HughesPJClass
15import Data.List (intercalate)
16 14
17 15
18brusselator _t x = [ a - (w + 1) * u + v * u^2 16brusselator :: Double -> [Double] -> [Double]
19 , w * u - v * u^2 17brusselator _t x = [ a - (w + 1) * u + v * u * u
18 , w * u - v * u * u
20 , (b - w) / eps - w * u 19 , (b - w) / eps - w * u
21 ] 20 ]
22 where 21 where
@@ -27,6 +26,7 @@ brusselator _t x = [ a - (w + 1) * u + v * u^2
27 v = x !! 1 26 v = x !! 1
28 w = x !! 2 27 w = x !! 2
29 28
29brussJac :: Double -> Vector Double -> Matrix Double
30brussJac _t x = (3><3) [ (-(w + 1.0)) + 2.0 * u * v, w - 2.0 * u * v, (-w) 30brussJac _t x = (3><3) [ (-(w + 1.0)) + 2.0 * u * v, w - 2.0 * u * v, (-w)
31 , u * u , (-(u * u)) , 0.0 31 , u * u , (-(u * u)) , 0.0
32 , (-u) , u , (-1.0) / eps - u 32 , (-u) , u , (-1.0) / eps - u
@@ -38,11 +38,13 @@ brussJac _t x = (3><3) [ (-(w + 1.0)) + 2.0 * u * v, w - 2.0 * u * v, (-w)
38 w = y !! 2 38 w = y !! 2
39 eps = 5.0e-6 39 eps = 5.0e-6
40 40
41stiffish :: Double -> [Double] -> [Double]
41stiffish t v = [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ] 42stiffish t v = [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ]
42 where 43 where
43 lamda = -100.0 44 lamda = -100.0
44 u = v !! 0 45 u = v !! 0
45 46
47stiffJac :: Double -> Vector Double -> Matrix Double
46stiffJac _t _v = (1><1) [ lamda ] 48stiffJac _t _v = (1><1) [ lamda ]
47 where 49 where
48 lamda = -100.0 50 lamda = -100.0
@@ -71,7 +73,7 @@ butcherTableauTex m = render $
71 n = rows m 73 n = rows m
72 rs = toLists m 74 rs = toLists m
73 ss = map (\r -> intercalate " & " $ map show r) rs 75 ss = map (\r -> intercalate " & " $ map show r) rs
74 ts = zipWith (\n r -> "c_" ++ show n ++ " & " ++ r) [1..n] ss 76 ts = zipWith (\i r -> "c_" ++ show i ++ " & " ++ r) [1..n] ss
75 us = vcat $ map (\r -> text r <+> text "\\\\") ts 77 us = vcat $ map (\r -> text r <+> text "\\\\") ts
76 78
77main :: IO () 79main :: IO ()