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