diff options
Diffstat (limited to 'packages/sundials/src/Main.hs')
-rw-r--r-- | packages/sundials/src/Main.hs | 186 |
1 files changed, 0 insertions, 186 deletions
diff --git a/packages/sundials/src/Main.hs b/packages/sundials/src/Main.hs deleted file mode 100644 index 16c21c5..0000000 --- a/packages/sundials/src/Main.hs +++ /dev/null | |||
@@ -1,186 +0,0 @@ | |||
1 | {-# OPTIONS_GHC -Wall #-} | ||
2 | |||
3 | import qualified Numeric.Sundials.ARKode.ODE as ARK | ||
4 | import qualified Numeric.Sundials.CVode.ODE as CV | ||
5 | import Numeric.LinearAlgebra | ||
6 | |||
7 | import Plots as P | ||
8 | import qualified Diagrams.Prelude as D | ||
9 | import Diagrams.Backend.Rasterific | ||
10 | |||
11 | import Control.Lens | ||
12 | |||
13 | import Test.Hspec | ||
14 | |||
15 | |||
16 | lorenz :: Double -> [Double] -> [Double] | ||
17 | lorenz _t u = [ sigma * (y - x) | ||
18 | , x * (rho - z) - y | ||
19 | , x * y - beta * z | ||
20 | ] | ||
21 | where | ||
22 | rho = 28.0 | ||
23 | sigma = 10.0 | ||
24 | beta = 8.0 / 3.0 | ||
25 | x = u !! 0 | ||
26 | y = u !! 1 | ||
27 | z = u !! 2 | ||
28 | |||
29 | _lorenzJac :: Double -> Vector Double -> Matrix Double | ||
30 | _lorenzJac _t u = (3><3) [ (-sigma), rho - z, y | ||
31 | , sigma , -1.0 , x | ||
32 | , 0.0 , (-x) , (-beta) | ||
33 | ] | ||
34 | where | ||
35 | rho = 28.0 | ||
36 | sigma = 10.0 | ||
37 | beta = 8.0 / 3.0 | ||
38 | x = u ! 0 | ||
39 | y = u ! 1 | ||
40 | z = u ! 2 | ||
41 | |||
42 | brusselator :: Double -> [Double] -> [Double] | ||
43 | brusselator _t x = [ a - (w + 1) * u + v * u * u | ||
44 | , w * u - v * u * u | ||
45 | , (b - w) / eps - w * u | ||
46 | ] | ||
47 | where | ||
48 | a = 1.0 | ||
49 | b = 3.5 | ||
50 | eps = 5.0e-6 | ||
51 | u = x !! 0 | ||
52 | v = x !! 1 | ||
53 | w = x !! 2 | ||
54 | |||
55 | _brussJac :: Double -> Vector Double -> Matrix Double | ||
56 | _brussJac _t x = (3><3) [ (-(w + 1.0)) + 2.0 * u * v, w - 2.0 * u * v, (-w) | ||
57 | , u * u , (-(u * u)) , 0.0 | ||
58 | , (-u) , u , (-1.0) / eps - u | ||
59 | ] | ||
60 | where | ||
61 | y = toList x | ||
62 | u = y !! 0 | ||
63 | v = y !! 1 | ||
64 | w = y !! 2 | ||
65 | eps = 5.0e-6 | ||
66 | |||
67 | stiffish :: Double -> [Double] -> [Double] | ||
68 | stiffish t v = [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ] | ||
69 | where | ||
70 | lamda = -100.0 | ||
71 | u = v !! 0 | ||
72 | |||
73 | stiffishV :: Double -> Vector Double -> Vector Double | ||
74 | stiffishV t v = fromList [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ] | ||
75 | where | ||
76 | lamda = -100.0 | ||
77 | u = v ! 0 | ||
78 | |||
79 | _stiffJac :: Double -> Vector Double -> Matrix Double | ||
80 | _stiffJac _t _v = (1><1) [ lamda ] | ||
81 | where | ||
82 | lamda = -100.0 | ||
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 | |||
101 | lSaxis :: [[Double]] -> P.Axis B D.V2 Double | ||
102 | lSaxis xs = P.r2Axis &~ do | ||
103 | let ts = xs!!0 | ||
104 | us = xs!!1 | ||
105 | vs = xs!!2 | ||
106 | ws = xs!!3 | ||
107 | P.linePlot' $ zip ts us | ||
108 | P.linePlot' $ zip ts vs | ||
109 | P.linePlot' $ zip ts ws | ||
110 | |||
111 | kSaxis :: [(Double, Double)] -> P.Axis B D.V2 Double | ||
112 | kSaxis xs = P.r2Axis &~ do | ||
113 | P.linePlot' xs | ||
114 | |||
115 | main :: IO () | ||
116 | main = do | ||
117 | |||
118 | let res1 = ARK.odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0]) | ||
119 | renderRasterific "diagrams/brusselator.png" | ||
120 | (D.dims2D 500.0 500.0) | ||
121 | (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1)) | ||
122 | |||
123 | let res1a = ARK.odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0]) | ||
124 | renderRasterific "diagrams/brusselatorA.png" | ||
125 | (D.dims2D 500.0 500.0) | ||
126 | (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1a)) | ||
127 | |||
128 | let res2 = ARK.odeSolve stiffish [0.0] (fromList [0.0, 0.1 .. 10.0]) | ||
129 | renderRasterific "diagrams/stiffish.png" | ||
130 | (D.dims2D 500.0 500.0) | ||
131 | (renderAxis $ kSaxis $ zip [0.0, 0.1 .. 10.0] (concat $ toLists res2)) | ||
132 | |||
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]) | ||
134 | |||
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]) | ||
136 | |||
137 | let maxDiffA = maximum $ map abs $ | ||
138 | zipWith (-) ((toLists $ tr res2a)!!0) ((toLists $ tr res2b)!!0) | ||
139 | |||
140 | let res2c = CV.odeSolveV (CV.BDF) Nothing 1e-6 1e-10 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0]) | ||
141 | |||
142 | let maxDiffB = maximum $ map abs $ | ||
143 | zipWith (-) ((toLists $ tr res2a)!!0) ((toLists $ tr res2c)!!0) | ||
144 | |||
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]) | ||
149 | |||
150 | renderRasterific "diagrams/lorenz.png" | ||
151 | (D.dims2D 500.0 500.0) | ||
152 | (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!0) ((toLists $ tr res3)!!1)) | ||
153 | |||
154 | renderRasterific "diagrams/lorenz1.png" | ||
155 | (D.dims2D 500.0 500.0) | ||
156 | (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!0) ((toLists $ tr res3)!!2)) | ||
157 | |||
158 | renderRasterific "diagrams/lorenz2.png" | ||
159 | (D.dims2D 500.0 500.0) | ||
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 | |||