diff options
Diffstat (limited to 'packages/sundials/src/Main.hs')
-rw-r--r-- | packages/sundials/src/Main.hs | 72 |
1 files changed, 60 insertions, 12 deletions
diff --git a/packages/sundials/src/Main.hs b/packages/sundials/src/Main.hs index 729d35a..16c21c5 100644 --- a/packages/sundials/src/Main.hs +++ b/packages/sundials/src/Main.hs | |||
@@ -1,6 +1,7 @@ | |||
1 | {-# OPTIONS_GHC -Wall #-} | 1 | {-# OPTIONS_GHC -Wall #-} |
2 | 2 | ||
3 | import Numeric.Sundials.ARKode.ODE | 3 | import qualified Numeric.Sundials.ARKode.ODE as ARK |
4 | import qualified Numeric.Sundials.CVode.ODE as CV | ||
4 | import Numeric.LinearAlgebra | 5 | import Numeric.LinearAlgebra |
5 | 6 | ||
6 | import Plots as P | 7 | import Plots as P |
@@ -80,6 +81,23 @@ _stiffJac _t _v = (1><1) [ lamda ] | |||
80 | where | 81 | where |
81 | lamda = -100.0 | 82 | lamda = -100.0 |
82 | 83 | ||
84 | predatorPrey :: Double -> [Double] -> [Double] | ||
85 | predatorPrey _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 | |||
83 | lSaxis :: [[Double]] -> P.Axis B D.V2 Double | 101 | lSaxis :: [[Double]] -> P.Axis B D.V2 Double |
84 | lSaxis xs = P.r2Axis &~ do | 102 | lSaxis xs = P.r2Axis &~ do |
85 | let ts = xs!!0 | 103 | let ts = xs!!0 |
@@ -97,33 +115,37 @@ kSaxis xs = P.r2Axis &~ do | |||
97 | main :: IO () | 115 | main :: IO () |
98 | main = do | 116 | main = do |
99 | 117 | ||
100 | let res1 = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0]) | 118 | let res1 = ARK.odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0]) |
101 | renderRasterific "diagrams/brusselator.png" | 119 | renderRasterific "diagrams/brusselator.png" |
102 | (D.dims2D 500.0 500.0) | 120 | (D.dims2D 500.0 500.0) |
103 | (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1)) | 121 | (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1)) |
104 | 122 | ||
105 | let res1a = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0]) | 123 | let res1a = ARK.odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0]) |
106 | renderRasterific "diagrams/brusselatorA.png" | 124 | renderRasterific "diagrams/brusselatorA.png" |
107 | (D.dims2D 500.0 500.0) | 125 | (D.dims2D 500.0 500.0) |
108 | (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1a)) | 126 | (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1a)) |
109 | 127 | ||
110 | let res2 = odeSolve stiffish [0.0] (fromList [0.0, 0.1 .. 10.0]) | 128 | let res2 = ARK.odeSolve stiffish [0.0] (fromList [0.0, 0.1 .. 10.0]) |
111 | renderRasterific "diagrams/stiffish.png" | 129 | renderRasterific "diagrams/stiffish.png" |
112 | (D.dims2D 500.0 500.0) | 130 | (D.dims2D 500.0 500.0) |
113 | (renderAxis $ kSaxis $ zip [0.0, 0.1 .. 10.0] (concat $ toLists res2)) | 131 | (renderAxis $ kSaxis $ zip [0.0, 0.1 .. 10.0] (concat $ toLists res2)) |
114 | 132 | ||
115 | let res2a = odeSolveV (SDIRK_5_3_4') Nothing 1e-6 1e-10 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0]) | 133 | let res2a = ARK.odeSolveV (ARK.SDIRK_5_3_4') Nothing 1e-6 1e-10 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0]) |
116 | 134 | ||
117 | let res2b = odeSolveV (TRBDF2_3_3_2') Nothing 1e-6 1e-10 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0]) | 135 | let res2b = ARK.odeSolveV (ARK.TRBDF2_3_3_2') Nothing 1e-6 1e-10 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0]) |
118 | 136 | ||
119 | let maxDiff = maximum $ map abs $ | 137 | let maxDiffA = maximum $ map abs $ |
120 | zipWith (-) ((toLists $ tr res2a)!!0) ((toLists $ tr res2b)!!0) | 138 | zipWith (-) ((toLists $ tr res2a)!!0) ((toLists $ tr res2b)!!0) |
121 | 139 | ||
122 | hspec $ describe "Compare results" $ do | 140 | let res2c = CV.odeSolveV (CV.BDF) Nothing 1e-6 1e-10 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0]) |
123 | it "for two different RK methods" $ | 141 | |
124 | maxDiff < 1.0e-6 | 142 | let maxDiffB = maximum $ map abs $ |
143 | zipWith (-) ((toLists $ tr res2a)!!0) ((toLists $ tr res2c)!!0) | ||
125 | 144 | ||
126 | let res3 = odeSolve lorenz [-5.0, -5.0, 1.0] (fromList [0.0, 0.01 .. 10.0]) | 145 | let maxDiffC = maximum $ map abs $ |
146 | zipWith (-) ((toLists $ tr res2b)!!0) ((toLists $ tr res2c)!!0) | ||
147 | |||
148 | let res3 = ARK.odeSolve lorenz [-5.0, -5.0, 1.0] (fromList [0.0, 0.01 .. 10.0]) | ||
127 | 149 | ||
128 | renderRasterific "diagrams/lorenz.png" | 150 | renderRasterific "diagrams/lorenz.png" |
129 | (D.dims2D 500.0 500.0) | 151 | (D.dims2D 500.0 500.0) |
@@ -136,3 +158,29 @@ main = do | |||
136 | renderRasterific "diagrams/lorenz2.png" | 158 | renderRasterific "diagrams/lorenz2.png" |
137 | (D.dims2D 500.0 500.0) | 159 | (D.dims2D 500.0 500.0) |
138 | (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 | |||