1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
|
{-# OPTIONS_GHC -Wall #-}
import qualified Numeric.Sundials.ARKode.ODE as ARK
import qualified Numeric.Sundials.CVode.ODE as CV
import Numeric.LinearAlgebra
import Plots as P
import qualified Diagrams.Prelude as D
import Diagrams.Backend.Rasterific
import Control.Lens
import Test.Hspec
lorenz :: Double -> [Double] -> [Double]
lorenz _t u = [ sigma * (y - x)
, x * (rho - z) - y
, x * y - beta * z
]
where
rho = 28.0
sigma = 10.0
beta = 8.0 / 3.0
x = u !! 0
y = u !! 1
z = u !! 2
_lorenzJac :: Double -> Vector Double -> Matrix Double
_lorenzJac _t u = (3><3) [ (-sigma), rho - z, y
, sigma , -1.0 , x
, 0.0 , (-x) , (-beta)
]
where
rho = 28.0
sigma = 10.0
beta = 8.0 / 3.0
x = u ! 0
y = u ! 1
z = u ! 2
brusselator :: Double -> [Double] -> [Double]
brusselator _t x = [ a - (w + 1) * u + v * u * u
, w * u - v * u * u
, (b - w) / eps - w * u
]
where
a = 1.0
b = 3.5
eps = 5.0e-6
u = x !! 0
v = x !! 1
w = x !! 2
_brussJac :: Double -> Vector Double -> Matrix Double
_brussJac _t x = (3><3) [ (-(w + 1.0)) + 2.0 * u * v, w - 2.0 * u * v, (-w)
, u * u , (-(u * u)) , 0.0
, (-u) , u , (-1.0) / eps - u
]
where
y = toList x
u = y !! 0
v = y !! 1
w = y !! 2
eps = 5.0e-6
stiffish :: Double -> [Double] -> [Double]
stiffish t v = [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ]
where
lamda = -100.0
u = v !! 0
stiffishV :: Double -> Vector Double -> Vector Double
stiffishV t v = fromList [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ]
where
lamda = -100.0
u = v ! 0
_stiffJac :: Double -> Vector Double -> Matrix Double
_stiffJac _t _v = (1><1) [ lamda ]
where
lamda = -100.0
predatorPrey :: Double -> [Double] -> [Double]
predatorPrey _t v = [ x * a - b * x * y
, d * x * y - c * y - e * y * z
, (-f) * z + g * y * z
]
where
x = v!!0
y = v!!1
z = v!!2
a = 1.0
b = 1.0
c = 1.0
d = 1.0
e = 1.0
f = 1.0
g = 1.0
lSaxis :: [[Double]] -> P.Axis B D.V2 Double
lSaxis xs = P.r2Axis &~ do
let ts = xs!!0
us = xs!!1
vs = xs!!2
ws = xs!!3
P.linePlot' $ zip ts us
P.linePlot' $ zip ts vs
P.linePlot' $ zip ts ws
kSaxis :: [(Double, Double)] -> P.Axis B D.V2 Double
kSaxis xs = P.r2Axis &~ do
P.linePlot' xs
main :: IO ()
main = do
let res1 = ARK.odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
renderRasterific "diagrams/brusselator.png"
(D.dims2D 500.0 500.0)
(renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1))
let res1a = ARK.odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
renderRasterific "diagrams/brusselatorA.png"
(D.dims2D 500.0 500.0)
(renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1a))
let res2 = ARK.odeSolve stiffish [0.0] (fromList [0.0, 0.1 .. 10.0])
renderRasterific "diagrams/stiffish.png"
(D.dims2D 500.0 500.0)
(renderAxis $ kSaxis $ zip [0.0, 0.1 .. 10.0] (concat $ toLists res2))
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])
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])
let maxDiffA = maximum $ map abs $
zipWith (-) ((toLists $ tr res2a)!!0) ((toLists $ tr res2b)!!0)
let res2c = CV.odeSolveV (CV.BDF) Nothing 1e-6 1e-10 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0])
let maxDiffB = maximum $ map abs $
zipWith (-) ((toLists $ tr res2a)!!0) ((toLists $ tr res2c)!!0)
let maxDiffC = maximum $ map abs $
zipWith (-) ((toLists $ tr res2b)!!0) ((toLists $ tr res2c)!!0)
let res3 = ARK.odeSolve lorenz [-5.0, -5.0, 1.0] (fromList [0.0, 0.01 .. 10.0])
renderRasterific "diagrams/lorenz.png"
(D.dims2D 500.0 500.0)
(renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!0) ((toLists $ tr res3)!!1))
renderRasterific "diagrams/lorenz1.png"
(D.dims2D 500.0 500.0)
(renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!0) ((toLists $ tr res3)!!2))
renderRasterific "diagrams/lorenz2.png"
(D.dims2D 500.0 500.0)
(renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!1) ((toLists $ tr res3)!!2))
let res4 = CV.odeSolve predatorPrey [0.5, 1.0, 2.0] (fromList [0.0, 0.01 .. 10.0])
renderRasterific "diagrams/predatorPrey.png"
(D.dims2D 500.0 500.0)
(renderAxis $ kSaxis $ zip ((toLists $ tr res4)!!0) ((toLists $ tr res4)!!1))
renderRasterific "diagrams/predatorPrey1.png"
(D.dims2D 500.0 500.0)
(renderAxis $ kSaxis $ zip ((toLists $ tr res4)!!0) ((toLists $ tr res4)!!2))
renderRasterific "diagrams/predatorPrey2.png"
(D.dims2D 500.0 500.0)
(renderAxis $ kSaxis $ zip ((toLists $ tr res4)!!1) ((toLists $ tr res4)!!2))
let res4a = ARK.odeSolve predatorPrey [0.5, 1.0, 2.0] (fromList [0.0, 0.01 .. 10.0])
let maxDiffPpA = maximum $ map abs $
zipWith (-) ((toLists $ tr res4)!!0) ((toLists $ tr res4a)!!0)
hspec $ describe "Compare results" $ do
it "for SDIRK_5_3_4' and TRBDF2_3_3_2'" $ maxDiffA < 1.0e-6
it "for SDIRK_5_3_4' and BDF" $ maxDiffB < 1.0e-6
it "for TRBDF2_3_3_2' and BDF" $ maxDiffC < 1.0e-6
it "for CV and ARK for the Predator Prey model" $ maxDiffPpA < 1.0e-3
|
