diff options
Diffstat (limited to 'packages/sundials/src')
-rw-r--r-- | packages/sundials/src/Main.hs | 68 |
1 files changed, 15 insertions, 53 deletions
diff --git a/packages/sundials/src/Main.hs b/packages/sundials/src/Main.hs index 78921a7..729d35a 100644 --- a/packages/sundials/src/Main.hs +++ b/packages/sundials/src/Main.hs | |||
@@ -8,9 +8,9 @@ import qualified Diagrams.Prelude as D | |||
8 | import Diagrams.Backend.Rasterific | 8 | import Diagrams.Backend.Rasterific |
9 | 9 | ||
10 | import Control.Lens | 10 | import Control.Lens |
11 | import Data.List (intercalate) | ||
12 | 11 | ||
13 | import Text.PrettyPrint.HughesPJClass | 12 | import Test.Hspec |
13 | |||
14 | 14 | ||
15 | lorenz :: Double -> [Double] -> [Double] | 15 | lorenz :: Double -> [Double] -> [Double] |
16 | lorenz _t u = [ sigma * (y - x) | 16 | lorenz _t u = [ sigma * (y - x) |
@@ -75,8 +75,8 @@ stiffishV t v = fromList [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ] | |||
75 | lamda = -100.0 | 75 | lamda = -100.0 |
76 | u = v ! 0 | 76 | u = v ! 0 |
77 | 77 | ||
78 | stiffJac :: Double -> Vector Double -> Matrix Double | 78 | _stiffJac :: Double -> Vector Double -> Matrix Double |
79 | stiffJac _t _v = (1><1) [ lamda ] | 79 | _stiffJac _t _v = (1><1) [ lamda ] |
80 | where | 80 | where |
81 | lamda = -100.0 | 81 | lamda = -100.0 |
82 | 82 | ||
@@ -94,44 +94,9 @@ kSaxis :: [(Double, Double)] -> P.Axis B D.V2 Double | |||
94 | kSaxis xs = P.r2Axis &~ do | 94 | kSaxis xs = P.r2Axis &~ do |
95 | P.linePlot' xs | 95 | P.linePlot' xs |
96 | 96 | ||
97 | butcherTableauTex :: ButcherTable -> String | ||
98 | butcherTableauTex (ButcherTable m c b b2) = | ||
99 | render $ | ||
100 | vcat [ text ("\n\\begin{array}{c|" ++ (concat $ replicate n "c") ++ "}") | ||
101 | , us | ||
102 | , text "\\hline" | ||
103 | , text bs <+> text "\\\\" | ||
104 | , text b2s <+> text "\\\\" | ||
105 | , text "\\end{array}" | ||
106 | ] | ||
107 | where | ||
108 | n = rows m | ||
109 | rs = toLists m | ||
110 | ss = map (\r -> intercalate " & " $ map show r) rs | ||
111 | ts = zipWith (\i r -> show i ++ " & " ++ r) (toList c) ss | ||
112 | us = vcat $ map (\r -> text r <+> text "\\\\") ts | ||
113 | bs = " & " ++ (intercalate " & " $ map show $ toList b) | ||
114 | b2s = " & " ++ (intercalate " & " $ map show $ toList b2) | ||
115 | |||
116 | main :: IO () | 97 | main :: IO () |
117 | main = do | 98 | main = do |
118 | 99 | ||
119 | let res = butcherTable (SDIRK_2_1_2 undefined) | ||
120 | putStrLn $ show res | ||
121 | putStrLn $ butcherTableauTex res | ||
122 | |||
123 | let resA = butcherTable (KVAERNO_4_2_3 undefined) | ||
124 | putStrLn $ show resA | ||
125 | putStrLn $ butcherTableauTex resA | ||
126 | |||
127 | let resB = butcherTable (SDIRK_5_3_4 undefined) | ||
128 | putStrLn $ show resB | ||
129 | putStrLn $ butcherTableauTex resB | ||
130 | |||
131 | let resC = butcherTable (FEHLBERG_6_4_5 undefined) | ||
132 | putStrLn $ show resC | ||
133 | putStrLn $ butcherTableauTex resC | ||
134 | |||
135 | let res1 = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0]) | 100 | let res1 = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0]) |
136 | renderRasterific "diagrams/brusselator.png" | 101 | renderRasterific "diagrams/brusselator.png" |
137 | (D.dims2D 500.0 500.0) | 102 | (D.dims2D 500.0 500.0) |
@@ -143,29 +108,26 @@ main = do | |||
143 | (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1a)) | 108 | (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1a)) |
144 | 109 | ||
145 | let res2 = odeSolve stiffish [0.0] (fromList [0.0, 0.1 .. 10.0]) | 110 | let res2 = odeSolve stiffish [0.0] (fromList [0.0, 0.1 .. 10.0]) |
146 | putStrLn $ show res2 | ||
147 | renderRasterific "diagrams/stiffish.png" | 111 | renderRasterific "diagrams/stiffish.png" |
148 | (D.dims2D 500.0 500.0) | 112 | (D.dims2D 500.0 500.0) |
149 | (renderAxis $ kSaxis $ zip [0.0, 0.1 .. 10.0] (concat $ toLists res2)) | 113 | (renderAxis $ kSaxis $ zip [0.0, 0.1 .. 10.0] (concat $ toLists res2)) |
150 | 114 | ||
151 | let res2a = odeSolveV (SDIRK_5_3_4 stiffJac) Nothing 1e-3 1e-6 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0]) | 115 | let res2a = odeSolveV (SDIRK_5_3_4') Nothing 1e-6 1e-10 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0]) |
152 | putStrLn "Lower tolerances" | ||
153 | putStrLn $ show res2a | ||
154 | 116 | ||
155 | let res3 = odeSolve lorenz [-5.0, -5.0, 1.0] (fromList [0.0, 0.01 .. 10.0]) | 117 | let res2b = odeSolveV (TRBDF2_3_3_2') Nothing 1e-6 1e-10 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0]) |
156 | putStrLn $ show $ last ((toLists $ tr res3)!!0) | ||
157 | 118 | ||
158 | let res3A = odeSolve lorenz [-5.0, -5.0, 1.0] (fromList [0.0, 0.01 .. 10.0]) | 119 | let maxDiff = maximum $ map abs $ |
159 | putStrLn $ show $ last ((toLists $ tr res3A)!!0) | 120 | zipWith (-) ((toLists $ tr res2a)!!0) ((toLists $ tr res2b)!!0) |
160 | 121 | ||
161 | renderRasterific "diagrams/lorenz.png" | 122 | hspec $ describe "Compare results" $ do |
162 | (D.dims2D 500.0 500.0) | 123 | it "for two different RK methods" $ |
163 | (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!0) ((toLists $ tr res3A)!!1)) | 124 | maxDiff < 1.0e-6 |
164 | 125 | ||
165 | let res3a = odeSolve lorenz [-5.0, -5.0, 1.0] (fromList [0.0, 0.01 .. 10.0]) | 126 | let res3 = odeSolve lorenz [-5.0, -5.0, 1.0] (fromList [0.0, 0.01 .. 10.0]) |
166 | renderRasterific "diagrams/lorenzA.png" | 127 | |
128 | renderRasterific "diagrams/lorenz.png" | ||
167 | (D.dims2D 500.0 500.0) | 129 | (D.dims2D 500.0 500.0) |
168 | (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!0) ((toLists $ tr res3a)!!1)) | 130 | (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!0) ((toLists $ tr res3)!!1)) |
169 | 131 | ||
170 | renderRasterific "diagrams/lorenz1.png" | 132 | renderRasterific "diagrams/lorenz1.png" |
171 | (D.dims2D 500.0 500.0) | 133 | (D.dims2D 500.0 500.0) |