summaryrefslogtreecommitdiff
path: root/packages
diff options
context:
space:
mode:
authorDominic Steinitz <dominic@steinitz.org>2018-04-08 10:38:59 +0100
committerDominic Steinitz <dominic@steinitz.org>2018-04-08 10:38:59 +0100
commit56c51c7476beeab5f1b701da8339f9bc1422cd52 (patch)
tree709c2f103097fbe0e35b45ccbf9f73cca34a24b1 /packages
parentc15409b09cb5133bf7f226aa8c8de06c9dbfdbd2 (diff)
Tidy and add picture
Diffstat (limited to 'packages')
-rw-r--r--packages/sundials/hmatrix-sundials.cabal9
-rw-r--r--packages/sundials/src/Main.hs12
-rw-r--r--packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs18
3 files changed, 13 insertions, 26 deletions
diff --git a/packages/sundials/hmatrix-sundials.cabal b/packages/sundials/hmatrix-sundials.cabal
index 53afa58..5aa4373 100644
--- a/packages/sundials/hmatrix-sundials.cabal
+++ b/packages/sundials/hmatrix-sundials.cabal
@@ -1,16 +1,17 @@
1name: hmatrix-sundials 1name: hmatrix-sundials
2version: 0.1.0.0 2version: 0.1.0.0
3-- synopsis: 3synopsis: hmatrix interface to sundials
4-- description: 4description: Foo bar baz
5homepage: https://github.com/idontgetoutmuch/hmatrix/tree/sundials 5homepage: https://github.com/idontgetoutmuch/hmatrix/tree/sundials
6license: BSD3 6license: BSD3
7license-file: LICENSE 7license-file: LICENSE
8author: Dominic Steinitz 8author: Dominic Steinitz
9maintainer: dominic@steinitz.org 9maintainer: dominic@steinitz.org
10-- copyright: 10copyright: Dominic Steinitz 2018, Novadiscovery 2018
11category: Math 11category: Math
12build-type: Simple 12build-type: Simple
13extra-source-files: ChangeLog.md 13extra-source-files: ChangeLog.md, README.md, diagrams/*.png
14extra-doc-files: diagrams/*.png
14cabal-version: >=1.10 15cabal-version: >=1.10
15 16
16 17
diff --git a/packages/sundials/src/Main.hs b/packages/sundials/src/Main.hs
index 6d6a397..ac19e7f 100644
--- a/packages/sundials/src/Main.hs
+++ b/packages/sundials/src/Main.hs
@@ -133,17 +133,17 @@ main = do
133 -- putStrLn $ show res 133 -- putStrLn $ show res
134 -- putStrLn $ butcherTableauTex res 134 -- putStrLn $ butcherTableauTex res
135 135
136 let res1 = odeSolve' (SDIRK_5_3_4 brussJac) brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0]) 136 let res1 = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
137 renderRasterific "diagrams/brusselator.png" 137 renderRasterific "diagrams/brusselator.png"
138 (D.dims2D 500.0 500.0) 138 (D.dims2D 500.0 500.0)
139 (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1)) 139 (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1))
140 140
141 let res1a = odeSolve' (SDIRK_5_3_4') brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0]) 141 let res1a = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
142 renderRasterific "diagrams/brusselatorA.png" 142 renderRasterific "diagrams/brusselatorA.png"
143 (D.dims2D 500.0 500.0) 143 (D.dims2D 500.0 500.0)
144 (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1a)) 144 (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1a))
145 145
146 let res2 = odeSolve' (SDIRK_5_3_4 stiffJac) stiffish [0.0] (fromList [0.0, 0.1 .. 10.0]) 146 let res2 = odeSolve stiffish [0.0] (fromList [0.0, 0.1 .. 10.0])
147 putStrLn $ show res2 147 putStrLn $ show res2
148 renderRasterific "diagrams/stiffish.png" 148 renderRasterific "diagrams/stiffish.png"
149 (D.dims2D 500.0 500.0) 149 (D.dims2D 500.0 500.0)
@@ -153,17 +153,17 @@ main = do
153 putStrLn "Lower tolerances" 153 putStrLn "Lower tolerances"
154 putStrLn $ show res2a 154 putStrLn $ show res2a
155 155
156 let res3 = odeSolve' (SDIRK_5_3_4 lorenzJac) lorenz [-5.0, -5.0, 1.0] (fromList [0.0, 0.01 .. 10.0]) 156 let res3 = odeSolve lorenz [-5.0, -5.0, 1.0] (fromList [0.0, 0.01 .. 10.0])
157 putStrLn $ show $ last ((toLists $ tr res3)!!0) 157 putStrLn $ show $ last ((toLists $ tr res3)!!0)
158 158
159 let res3 = odeSolve' (SDIRK_5_3_4') lorenz [-5.0, -5.0, 1.0] (fromList [0.0, 0.01 .. 10.0]) 159 let res3 = odeSolve lorenz [-5.0, -5.0, 1.0] (fromList [0.0, 0.01 .. 10.0])
160 putStrLn $ show $ last ((toLists $ tr res3)!!0) 160 putStrLn $ show $ last ((toLists $ tr res3)!!0)
161 161
162 renderRasterific "diagrams/lorenz.png" 162 renderRasterific "diagrams/lorenz.png"
163 (D.dims2D 500.0 500.0) 163 (D.dims2D 500.0 500.0)
164 (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!0) ((toLists $ tr res3)!!1)) 164 (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!0) ((toLists $ tr res3)!!1))
165 165
166 let res3a = odeSolve' (SDIRK_5_3_4') lorenz [-5.0, -5.0, 1.0] (fromList [0.0, 0.01 .. 10.0]) 166 let res3a = odeSolve lorenz [-5.0, -5.0, 1.0] (fromList [0.0, 0.01 .. 10.0])
167 renderRasterific "diagrams/lorenzA.png" 167 renderRasterific "diagrams/lorenzA.png"
168 (D.dims2D 500.0 500.0) 168 (D.dims2D 500.0 500.0)
169 (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!0) ((toLists $ tr res3a)!!1)) 169 (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!0) ((toLists $ tr res3a)!!1))
diff --git a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
index 8f83fe7..2577b8e 100644
--- a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
+++ b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
@@ -35,6 +35,8 @@
35-- main = mplot (ts : toColumns sol) 35-- main = mplot (ts : toColumns sol)
36-- @ 36-- @
37-- 37--
38-- <<diagrams/brusselator.png#diagram=brusselator&height=400&width=500>>
39--
38-- KVAERNO_4_2_3 40-- KVAERNO_4_2_3
39-- 41--
40-- \[ 42-- \[
@@ -70,7 +72,6 @@
70module Numeric.Sundials.ARKode.ODE ( odeSolve 72module Numeric.Sundials.ARKode.ODE ( odeSolve
71 , odeSolveV 73 , odeSolveV
72 , odeSolveVWith 74 , odeSolveVWith
73 , odeSolve'
74 , getButcherTable 75 , getButcherTable
75 , getBT 76 , getBT
76 , btGet 77 , btGet
@@ -244,21 +245,6 @@ odeSolve f y0 ts =
244 nC = length y0 245 nC = length y0
245 g t x0 = V.fromList $ f t (V.toList x0) 246 g t x0 = V.fromList $ f t (V.toList x0)
246 247
247odeSolve' :: ODEMethod
248 -> (Double -> [Double] -> [Double]) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
249 -> [Double] -- ^ initial conditions
250 -> Vector Double -- ^ desired solution times
251 -> Matrix Double -- ^ solution
252odeSolve' method f y0 ts =
253 case odeSolveVWith method (XX' 1.0e-6 1.0e-10 1 1) g (V.fromList y0) (V.fromList $ toList ts) of
254 Left c -> error $ show c -- FIXME
255 Right (v, d) -> trace (show d) $ (nR >< nC) (V.toList v)
256 where
257 us = toList ts
258 nR = length us
259 nC = length y0
260 g t x0 = V.fromList $ f t (V.toList x0)
261
262odeSolveVWith :: 248odeSolveVWith ::
263 ODEMethod 249 ODEMethod
264 -> StepControl 250 -> StepControl