blob: 5e513727d0e8244163a0becc8cefc40e38f8fb0e (
plain)
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
|
{-# OPTIONS_GHC -Wall #-}
import qualified Data.Vector.Storable as V
import Numeric.Sundials.Arkode.ODE
import Numeric.LinearAlgebra
import Plots as P
import qualified Diagrams.Prelude as D
import Diagrams.Backend.Rasterific
import Control.Lens
import Data.List (zip4)
brusselator _t x = [ a - (w + 1) * u + v * u^2
, w * u - v * u^2
, (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
stiffish t v = [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ]
where
lamda = -100.0
u = v !! 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 res = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
putStrLn $ show res
renderRasterific "diagrams/brusselator.png"
(D.dims2D 500.0 500.0)
(renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res))
let res = odeSolve stiffish [0.0] (fromList [0.0, 0.1 .. 10.0])
putStrLn $ show res
renderRasterific "diagrams/stiffish.png"
(D.dims2D 500.0 500.0)
(renderAxis $ kSaxis $ zip [0.0, 0.1 .. 10.0] (concat $ toLists res))
|
