| Copyright | Dominic Steinitz 2018 Novadiscovery 2018 |
|---|---|
| License | BSD |
| Maintainer | Dominic Steinitz |
| Stability | provisional |
| Safe Haskell | None |
| Language | Haskell2010 |
Numeric.Sundials.ARKode.ODE
Description
Solution of ordinary differential equation (ODE) initial value problems.
https://computation.llnl.gov/projects/sundials/sundials-software
A simple example:
import Numeric.Sundials.ARKode.ODE
import Numeric.LinearAlgebra
import Plots as P
import qualified Diagrams.Prelude as D
import Diagrams.Backend.Rasterific
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
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
main = do
let res1 = 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))
KVAERNO_4_2_3
\[ \begin{array}{c|cccc} 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\ 0.871733043 & 0.4358665215 & 0.4358665215 & 0.0 & 0.0 \\ 1.0 & 0.490563388419108 & 7.3570090080892e-2 & 0.4358665215 & 0.0 \\ 1.0 & 0.308809969973036 & 1.490563388254106 & -1.235239879727145 & 0.4358665215 \\ \hline & 0.308809969973036 & 1.490563388254106 & -1.235239879727145 & 0.4358665215 \\ & 0.490563388419108 & 7.3570090080892e-2 & 0.4358665215 & 0.0 \\ \end{array} \]
SDIRK_2_1_2
\[ \begin{array}{c|cc} 1.0 & 1.0 & 0.0 \\ 0.0 & -1.0 & 1.0 \\ \hline & 0.5 & 0.5 \\ & 1.0 & 0.0 \\ \end{array} \]
SDIRK_5_3_4
\[ \begin{array}{c|ccccc} 0.25 & 0.25 & 0.0 & 0.0 & 0.0 & 0.0 \\ 0.75 & 0.5 & 0.25 & 0.0 & 0.0 & 0.0 \\ 0.55 & 0.34 & -4.0e-2 & 0.25 & 0.0 & 0.0 \\ 0.5 & 0.2727941176470588 & -5.036764705882353e-2 & 2.7573529411764705e-2 & 0.25 & 0.0 \\ 1.0 & 1.0416666666666667 & -1.0208333333333333 & 7.8125 & -7.083333333333333 & 0.25 \\ \hline & 1.0416666666666667 & -1.0208333333333333 & 7.8125 & -7.083333333333333 & 0.25 \\ & 1.2291666666666667 & -0.17708333333333334 & 7.03125 & -7.083333333333333 & 0.0 \\ \end{array} \]
- odeSolve :: (Double -> [Double] -> [Double]) -> [Double] -> Vector Double -> Matrix Double
- odeSolveV :: ODEMethod -> Maybe Double -> Double -> Double -> (Double -> Vector Double -> Vector Double) -> Vector Double -> Vector Double -> Matrix Double
- odeSolveVWith :: ODEMethod -> StepControl -> Maybe Double -> (Double -> Vector Double -> Vector Double) -> Vector Double -> Vector Double -> Either Int (Vector Double, SundialsDiagnostics)
- odeSolveVWith' :: ODEMethod -> StepControl -> Maybe Double -> (Double -> Vector Double -> Vector Double) -> Vector Double -> Vector Double -> Matrix Double
- data ButcherTable = ButcherTable {}
- butcherTable :: ODEMethod -> ButcherTable
- data ODEMethod
- = SDIRK_2_1_2 Jacobian
- | SDIRK_2_1_2'
- | BILLINGTON_3_3_2 Jacobian
- | BILLINGTON_3_3_2'
- | TRBDF2_3_3_2 Jacobian
- | TRBDF2_3_3_2'
- | KVAERNO_4_2_3 Jacobian
- | KVAERNO_4_2_3'
- | ARK324L2SA_DIRK_4_2_3 Jacobian
- | ARK324L2SA_DIRK_4_2_3'
- | CASH_5_2_4 Jacobian
- | CASH_5_2_4'
- | CASH_5_3_4 Jacobian
- | CASH_5_3_4'
- | SDIRK_5_3_4 Jacobian
- | SDIRK_5_3_4'
- | KVAERNO_5_3_4 Jacobian
- | KVAERNO_5_3_4'
- | ARK436L2SA_DIRK_6_3_4 Jacobian
- | ARK436L2SA_DIRK_6_3_4'
- | KVAERNO_7_4_5 Jacobian
- | KVAERNO_7_4_5'
- | ARK548L2SA_DIRK_8_4_5 Jacobian
- | ARK548L2SA_DIRK_8_4_5'
- | HEUN_EULER_2_1_2 Jacobian
- | HEUN_EULER_2_1_2'
- | BOGACKI_SHAMPINE_4_2_3 Jacobian
- | BOGACKI_SHAMPINE_4_2_3'
- | ARK324L2SA_ERK_4_2_3 Jacobian
- | ARK324L2SA_ERK_4_2_3'
- | ZONNEVELD_5_3_4 Jacobian
- | ZONNEVELD_5_3_4'
- | ARK436L2SA_ERK_6_3_4 Jacobian
- | ARK436L2SA_ERK_6_3_4'
- | SAYFY_ABURUB_6_3_4 Jacobian
- | SAYFY_ABURUB_6_3_4'
- | CASH_KARP_6_4_5 Jacobian
- | CASH_KARP_6_4_5'
- | FEHLBERG_6_4_5 Jacobian
- | FEHLBERG_6_4_5'
- | DORMAND_PRINCE_7_4_5 Jacobian
- | DORMAND_PRINCE_7_4_5'
- | ARK548L2SA_ERK_8_4_5 Jacobian
- | ARK548L2SA_ERK_8_4_5'
- | VERNER_8_5_6 Jacobian
- | VERNER_8_5_6'
- | FEHLBERG_13_7_8 Jacobian
- | FEHLBERG_13_7_8'
- data StepControl
- type Jacobian = Double -> Vector Double -> Matrix Double
- data SundialsDiagnostics = SundialsDiagnostics {
- aRKodeGetNumSteps :: Int
- aRKodeGetNumStepAttempts :: Int
- aRKodeGetNumRhsEvals_fe :: Int
- aRKodeGetNumRhsEvals_fi :: Int
- aRKodeGetNumLinSolvSetups :: Int
- aRKodeGetNumErrTestFails :: Int
- aRKodeGetNumNonlinSolvIters :: Int
- aRKodeGetNumNonlinSolvConvFails :: Int
- aRKDlsGetNumJacEvals :: Int
- aRKDlsGetNumRhsEvals :: Int
Documentation
Arguments
| :: (Double -> [Double] -> [Double]) | The RHS of the system \(\dot{y} = f(t,y)\) |
| -> [Double] | initial conditions |
| -> Vector Double | desired solution times |
| -> Matrix Double | solution |
A version of odeSolveV with reasonable default parameters and
system of equations defined using lists. FIXME: we should say
something about the fact we could use the Jacobian but don't for
compatibility with hmatrix-gsl.
Arguments
| :: ODEMethod | |
| -> Maybe Double | initial step size - by default, ARKode estimates the initial step size to be the solution \(h\) of the equation \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where \(\ddot{y}\) is an estimated value of the second derivative of the solution at \(t_0\) |
| -> Double | absolute tolerance for the state vector |
| -> Double | relative tolerance for the state vector |
| -> (Double -> Vector Double -> Vector Double) | The RHS of the system \(\dot{y} = f(t,y)\) |
| -> Vector Double | initial conditions |
| -> Vector Double | desired solution times |
| -> Matrix Double | solution |
A version of odeSolveVWith with reasonable default step control.
Arguments
| :: ODEMethod | |
| -> StepControl | |
| -> Maybe Double | initial step size - by default, ARKode estimates the initial step size to be the solution \(h\) of the equation \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where \(\ddot{y}\) is an estimated value of the second derivative of the solution at \(t_0\) |
| -> (Double -> Vector Double -> Vector Double) | The RHS of the system \(\dot{y} = f(t,y)\) |
| -> Vector Double | Initial conditions |
| -> Vector Double | Desired solution times |
| -> Either Int (Vector Double, SundialsDiagnostics) | Error code or solution |
Arguments
| :: ODEMethod | |
| -> StepControl | |
| -> Maybe Double | initial step size - by default, ARKode estimates the initial step size to be the solution \(h\) of the equation \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where \(\ddot{y}\) is an estimated value of the second derivative of the solution at \(t_0\) |
| -> (Double -> Vector Double -> Vector Double) | The RHS of the system \(\dot{y} = f(t,y)\) |
| -> Vector Double | Initial conditions |
| -> Vector Double | Desired solution times |
| -> Matrix Double | Error code or solution |
butcherTable :: ODEMethod -> ButcherTable Source #
Stepping functions
Constructors
data StepControl Source #
Adaptive step-size control functions.
GSL allows the user to control the step size adjustment using \(D_i = \epsilon^{abs}s_i + \epsilon^{rel}(a_{y} |y_i| + a_{dy/dt} h |\dot{y}_i|)\) where \(\epsilon^{abs}\) is the required absolute error, \(\epsilon^{rel}\) is the required relative error, \(s_i\) is a vector of scaling factors, \(a_{y}\) is a scaling factor for the solution \(y\) and \(a_{dydt}\) is a scaling factor for the derivative of the solution \(dy/dt\).
ARKode allows the user to control the step size adjustment using \(\eta^{rel}|y_i| + \eta^{abs}_i\). For compatibility with hmatrix-gsl, tolerances for \(y\) and \(\dot{y}\) can be specified but the latter have no effect.
Constructors
| X Double Double | absolute and relative tolerance for \(y\); in GSL terms, \(a_{y} = 1\) and \(a_{dy/dt} = 0\); in ARKode terms, the \(\eta^{abs}_i\) are identical |
| X' Double Double | absolute and relative tolerance for \(\dot{y}\); in GSL terms, \(a_{y} = 0\) and \(a_{dy/dt} = 1\); in ARKode terms, the latter is treated as the relative tolerance for \(y\) so this is the same as specifying |
| XX' Double Double Double Double | include both via relative tolerance scaling factors \(a_y\), \(a_{{dy}/{dt}}\); in ARKode terms, the latter is ignored and \(\eta^{rel} = a_{y}\epsilon^{rel}\) |
| ScXX' Double Double Double Double (Vector Double) | scale absolute tolerance of \(y_i\); in ARKode terms, \(a_{{dy}/{dt}}\) is ignored, \(\eta^{abs}_i = s_i \epsilon^{abs}\) and \(\eta^{rel} = a_{y}\epsilon^{rel}\) |
data SundialsDiagnostics Source #
Constructors
Instances