summaryrefslogtreecommitdiff
path: root/packages/sundials/src/Main.hs
diff options
context:
space:
mode:
authorDominic Steinitz <dominic@steinitz.org>2018-05-02 14:42:43 +0100
committerDominic Steinitz <dominic@steinitz.org>2018-05-02 14:42:43 +0100
commit4ba859636396d211637b5507f19722b6953656a5 (patch)
tree9493c4851e6141a400e6345efe59a07197709f63 /packages/sundials/src/Main.hs
parent149dedfc6ec8dea039a4df7ad1d31880820c52eb (diff)
Add more options
Diffstat (limited to 'packages/sundials/src/Main.hs')
-rw-r--r--packages/sundials/src/Main.hs48
1 files changed, 43 insertions, 5 deletions
diff --git a/packages/sundials/src/Main.hs b/packages/sundials/src/Main.hs
index 85928e2..16c21c5 100644
--- a/packages/sundials/src/Main.hs
+++ b/packages/sundials/src/Main.hs
@@ -81,6 +81,23 @@ _stiffJac _t _v = (1><1) [ lamda ]
81 where 81 where
82 lamda = -100.0 82 lamda = -100.0
83 83
84predatorPrey :: Double -> [Double] -> [Double]
85predatorPrey _t v = [ x * a - b * x * y
86 , d * x * y - c * y - e * y * z
87 , (-f) * z + g * y * z
88 ]
89 where
90 x = v!!0
91 y = v!!1
92 z = v!!2
93 a = 1.0
94 b = 1.0
95 c = 1.0
96 d = 1.0
97 e = 1.0
98 f = 1.0
99 g = 1.0
100
84lSaxis :: [[Double]] -> P.Axis B D.V2 Double 101lSaxis :: [[Double]] -> P.Axis B D.V2 Double
85lSaxis xs = P.r2Axis &~ do 102lSaxis xs = P.r2Axis &~ do
86 let ts = xs!!0 103 let ts = xs!!0
@@ -128,11 +145,6 @@ main = do
128 let maxDiffC = maximum $ map abs $ 145 let maxDiffC = maximum $ map abs $
129 zipWith (-) ((toLists $ tr res2b)!!0) ((toLists $ tr res2c)!!0) 146 zipWith (-) ((toLists $ tr res2b)!!0) ((toLists $ tr res2c)!!0)
130 147
131 hspec $ describe "Compare results" $ do
132 it "for SDIRK_5_3_4' and TRBDF2_3_3_2'" $ maxDiffA < 1.0e-6
133 it "for SDIRK_5_3_4' and BDF" $ maxDiffB < 1.0e-6
134 it "for TRBDF2_3_3_2' and BDF" $ maxDiffC < 1.0e-6
135
136 let res3 = ARK.odeSolve lorenz [-5.0, -5.0, 1.0] (fromList [0.0, 0.01 .. 10.0]) 148 let res3 = ARK.odeSolve lorenz [-5.0, -5.0, 1.0] (fromList [0.0, 0.01 .. 10.0])
137 149
138 renderRasterific "diagrams/lorenz.png" 150 renderRasterific "diagrams/lorenz.png"
@@ -146,3 +158,29 @@ main = do
146 renderRasterific "diagrams/lorenz2.png" 158 renderRasterific "diagrams/lorenz2.png"
147 (D.dims2D 500.0 500.0) 159 (D.dims2D 500.0 500.0)
148 (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!1) ((toLists $ tr res3)!!2)) 160 (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!1) ((toLists $ tr res3)!!2))
161
162 let res4 = CV.odeSolve predatorPrey [0.5, 1.0, 2.0] (fromList [0.0, 0.01 .. 10.0])
163
164 renderRasterific "diagrams/predatorPrey.png"
165 (D.dims2D 500.0 500.0)
166 (renderAxis $ kSaxis $ zip ((toLists $ tr res4)!!0) ((toLists $ tr res4)!!1))
167
168 renderRasterific "diagrams/predatorPrey1.png"
169 (D.dims2D 500.0 500.0)
170 (renderAxis $ kSaxis $ zip ((toLists $ tr res4)!!0) ((toLists $ tr res4)!!2))
171
172 renderRasterific "diagrams/predatorPrey2.png"
173 (D.dims2D 500.0 500.0)
174 (renderAxis $ kSaxis $ zip ((toLists $ tr res4)!!1) ((toLists $ tr res4)!!2))
175
176 let res4a = ARK.odeSolve predatorPrey [0.5, 1.0, 2.0] (fromList [0.0, 0.01 .. 10.0])
177
178 let maxDiffPpA = maximum $ map abs $
179 zipWith (-) ((toLists $ tr res4)!!0) ((toLists $ tr res4a)!!0)
180
181 hspec $ describe "Compare results" $ do
182 it "for SDIRK_5_3_4' and TRBDF2_3_3_2'" $ maxDiffA < 1.0e-6
183 it "for SDIRK_5_3_4' and BDF" $ maxDiffB < 1.0e-6
184 it "for TRBDF2_3_3_2' and BDF" $ maxDiffC < 1.0e-6
185 it "for CV and ARK for the Predator Prey model" $ maxDiffPpA < 1.0e-3
186