summaryrefslogtreecommitdiff
path: root/packages/sundials
diff options
context:
space:
mode:
authorDominic Steinitz <dominic@steinitz.org>2018-04-17 13:50:56 +0100
committerDominic Steinitz <dominic@steinitz.org>2018-04-17 13:50:56 +0100
commit4cfa8c60c7652302904e33dd7cb8120b4596e1e3 (patch)
tree61a1b0436d44a8878b7a111df11964d30caedfd4 /packages/sundials
parent09f2d4cea649be1185183164fc65789c5944a2c8 (diff)
Support RKF45
Diffstat (limited to 'packages/sundials')
-rw-r--r--packages/sundials/src/Arkode.hsc25
-rw-r--r--packages/sundials/src/Main.hs4
-rw-r--r--packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs24
3 files changed, 44 insertions, 9 deletions
diff --git a/packages/sundials/src/Arkode.hsc b/packages/sundials/src/Arkode.hsc
index 2904c36..4ef95e2 100644
--- a/packages/sundials/src/Arkode.hsc
+++ b/packages/sundials/src/Arkode.hsc
@@ -46,6 +46,8 @@ arkSMax :: Int
46arkSMax = #const ARK_S_MAX 46arkSMax = #const ARK_S_MAX
47 47
48-- FIXME: We could just use inline-c instead 48-- FIXME: We could just use inline-c instead
49
50-- /* Butcher table accessors -- implicit */
49sDIRK_2_1_2 :: Int 51sDIRK_2_1_2 :: Int
50sDIRK_2_1_2 = #const SDIRK_2_1_2 52sDIRK_2_1_2 = #const SDIRK_2_1_2
51-- #define BILLINGTON_3_3_2 13 53-- #define BILLINGTON_3_3_2 13
@@ -67,3 +69,26 @@ sDIRK_5_3_4 = #const SDIRK_5_3_4
67-- #define DEFAULT_DIRK_3 ARK324L2SA_DIRK_4_2_3 69-- #define DEFAULT_DIRK_3 ARK324L2SA_DIRK_4_2_3
68-- #define DEFAULT_DIRK_4 SDIRK_5_3_4 70-- #define DEFAULT_DIRK_4 SDIRK_5_3_4
69-- #define DEFAULT_DIRK_5 ARK548L2SA_DIRK_8_4_5 71-- #define DEFAULT_DIRK_5 ARK548L2SA_DIRK_8_4_5
72
73-- /* Butcher table accessors -- explicit */
74-- #define HEUN_EULER_2_1_2 0
75-- #define BOGACKI_SHAMPINE_4_2_3 1
76-- #define ARK324L2SA_ERK_4_2_3 2
77-- #define ZONNEVELD_5_3_4 3
78-- #define ARK436L2SA_ERK_6_3_4 4
79-- #define SAYFY_ABURUB_6_3_4 5
80-- #define CASH_KARP_6_4_5 6
81fEHLBERG_6_4_5 :: Int
82fEHLBERG_6_4_5 = #const FEHLBERG_6_4_5
83-- #define FEHLBERG_6_4_5 7
84-- #define DORMAND_PRINCE_7_4_5 8
85-- #define ARK548L2SA_ERK_8_4_5 9
86-- #define VERNER_8_5_6 10
87-- #define FEHLBERG_13_7_8 11
88
89-- #define DEFAULT_ERK_2 HEUN_EULER_2_1_2
90-- #define DEFAULT_ERK_3 BOGACKI_SHAMPINE_4_2_3
91-- #define DEFAULT_ERK_4 ZONNEVELD_5_3_4
92-- #define DEFAULT_ERK_5 CASH_KARP_6_4_5
93-- #define DEFAULT_ERK_6 VERNER_8_5_6
94-- #define DEFAULT_ERK_8 FEHLBERG_13_7_8
diff --git a/packages/sundials/src/Main.hs b/packages/sundials/src/Main.hs
index e96522f..d8f3f3d 100644
--- a/packages/sundials/src/Main.hs
+++ b/packages/sundials/src/Main.hs
@@ -128,6 +128,10 @@ main = do
128 putStrLn $ show resB 128 putStrLn $ show resB
129 putStrLn $ butcherTableauTex resB 129 putStrLn $ butcherTableauTex resB
130 130
131 -- let resC = butcherTable (FEHLBERG_6_4_5 undefined)
132 -- putStrLn $ show resC
133 -- putStrLn $ butcherTableauTex resC
134
131 let res1 = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0]) 135 let res1 = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
132 renderRasterific "diagrams/brusselator.png" 136 renderRasterific "diagrams/brusselator.png"
133 (D.dims2D 500.0 500.0) 137 (D.dims2D 500.0 500.0)
diff --git a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
index 444138d..380b1d6 100644
--- a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
+++ b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
@@ -137,7 +137,7 @@ import Numeric.LinearAlgebra.HMatrix (Vector, Matrix, toList, (><),
137 size, subVector) 137 size, subVector)
138 138
139import qualified Types as T 139import qualified Types as T
140import Arkode (sDIRK_2_1_2, kVAERNO_4_2_3, sDIRK_5_3_4) 140import Arkode (sDIRK_2_1_2, kVAERNO_4_2_3, sDIRK_5_3_4, fEHLBERG_6_4_5)
141import qualified Arkode as B 141import qualified Arkode as B
142 142
143import Debug.Trace 143import Debug.Trace
@@ -225,18 +225,24 @@ data ODEMethod = SDIRK_2_1_2 Jacobian
225 | KVAERNO_4_2_3 Jacobian 225 | KVAERNO_4_2_3 Jacobian
226 | SDIRK_5_3_4 Jacobian 226 | SDIRK_5_3_4 Jacobian
227 | SDIRK_5_3_4' 227 | SDIRK_5_3_4'
228 | FEHLBERG_6_4_5 Jacobian
229 | FEHLBERG_6_4_5'
228 230
229getMethod :: ODEMethod -> Int 231getMethod :: ODEMethod -> Int
230getMethod (SDIRK_2_1_2 _) = sDIRK_2_1_2 232getMethod (SDIRK_2_1_2 _) = sDIRK_2_1_2
231getMethod (KVAERNO_4_2_3 _) = kVAERNO_4_2_3 233getMethod (KVAERNO_4_2_3 _) = kVAERNO_4_2_3
232getMethod (SDIRK_5_3_4 _) = sDIRK_5_3_4 234getMethod (SDIRK_5_3_4 _) = sDIRK_5_3_4
233getMethod (SDIRK_5_3_4' ) = sDIRK_5_3_4 235getMethod (SDIRK_5_3_4' ) = sDIRK_5_3_4
236getMethod (FEHLBERG_6_4_5 _) = fEHLBERG_6_4_5
237getMethod (FEHLBERG_6_4_5' ) = fEHLBERG_6_4_5
234 238
235getJacobian :: ODEMethod -> Maybe Jacobian 239getJacobian :: ODEMethod -> Maybe Jacobian
236getJacobian (SDIRK_2_1_2 j) = Just j 240getJacobian (SDIRK_2_1_2 j) = Just j
237getJacobian (KVAERNO_4_2_3 j) = Just j 241getJacobian (KVAERNO_4_2_3 j) = Just j
238getJacobian (SDIRK_5_3_4 j) = Just j 242getJacobian (SDIRK_5_3_4 j) = Just j
239getJacobian (SDIRK_5_3_4' ) = Nothing 243getJacobian (SDIRK_5_3_4' ) = Nothing
244getJacobian (FEHLBERG_6_4_5 j) = Just j
245getJacobian (FEHLBERG_6_4_5' ) = Nothing
240 246
241-- | A version of 'odeSolveVWith' with reasonable default step control. 247-- | A version of 'odeSolveVWith' with reasonable default step control.
242odeSolveV 248odeSolveV