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.hs186
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
3import qualified Numeric.Sundials.ARKode.ODE as ARK
4import qualified Numeric.Sundials.CVode.ODE as CV
5import Numeric.LinearAlgebra
6
7import Plots as P
8import qualified Diagrams.Prelude as D
9import Diagrams.Backend.Rasterific
10
11import Control.Lens
12
13import Test.Hspec
14
15
16lorenz :: Double -> [Double] -> [Double]
17lorenz _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
42brusselator :: Double -> [Double] -> [Double]
43brusselator _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
67stiffish :: Double -> [Double] -> [Double]
68stiffish t v = [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ]
69 where
70 lamda = -100.0
71 u = v !! 0
72
73stiffishV :: Double -> Vector Double -> Vector Double
74stiffishV 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
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
101lSaxis :: [[Double]] -> P.Axis B D.V2 Double
102lSaxis 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
111kSaxis :: [(Double, Double)] -> P.Axis B D.V2 Double
112kSaxis xs = P.r2Axis &~ do
113 P.linePlot' xs
114
115main :: IO ()
116main = 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
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-02 12:53:28 +0000