OrdinaryDiffEqSDIRK

Singly Diagonally Implicit Runge-Kutta (SDIRK) methods are a family of implicit Runge-Kutta methods designed for solving stiff ordinary differential equations. These methods are particularly effective for problems where explicit methods become unstable due to stiffness.

Key Properties

SDIRK methods have several important characteristics:

A-stable and L-stable: Can handle highly stiff problems without numerical instability
Stiffly accurate: Many SDIRK methods provide additional numerical stability for stiff problems
Diagonally implicit structure: The implicit system only requires solving a sequence of nonlinear equations rather than a large coupled system
Good for moderate to large systems: More efficient than fully implicit RK methods for many problems

When to Use SDIRK Methods

SDIRK methods are recommended for:

Stiff differential equations where explicit methods fail or require very small timesteps
Problems requiring good stability properties at moderate to high tolerances
Systems where Rosenbrock methods (which require Jacobians) are not suitable or available
IMEX problems using the KenCarp family, which can split stiff and non-stiff terms

Solver Selection Guide

High tolerance (>1e-2)

TRBDF2: Second-order A-B-L-S-stable method, good for oscillatory problems

Medium tolerance (1e-8 to 1e-2)

KenCarp4: Fourth-order method with excellent stability, good all-around choice
KenCarp47: Seventh-stage fourth-order method, enhanced stability
Kvaerno4 or Kvaerno5: High-order stiffly accurate methods

Low tolerance (<1e-8)

Kvaerno5: Fifth-order stiffly accurate method for high accuracy
KenCarp5: Fifth-order method with splitting capabilities

Special Cases

ImplicitEuler: First-order method, only recommended for problems with discontinuities or when f is not differentiable
Trapezoid: Second-order symmetric method, reversible but not symplectic. Good for eliminating damping often seen with L-stable methods
ImplicitMidpoint: Second-order A-stable symplectic method for energy-preserving systems
SSPSDIRK2: Strong stability preserving variant for problems requiring monotonicity preservation

Installation

To be able to access the solvers in OrdinaryDiffEqSDIRK, you must first install them use the Julia package manager:

using Pkg
Pkg.add("OrdinaryDiffEqSDIRK")

This will only install the solvers listed at the bottom of this page. If you want to explore other solvers for your problem, you will need to install some of the other libraries listed in the navigation bar on the left.

Example usage

using OrdinaryDiffEqSDIRK

function lorenz!(du, u, p, t)
    du[1] = 10.0 * (u[2] - u[1])
    du[2] = u[1] * (28.0 - u[3]) - u[2]
    du[3] = u[1] * u[2] - (8 / 3) * u[3]
end
u0 = [1.0; 0.0; 0.0]
tspan = (0.0, 100.0)
prob = ODEProblem(lorenz!, u0, tspan)
sol = solve(prob, KenCarp4())

Full list of solvers

OrdinaryDiffEqSDIRK.ImplicitEuler — Type

ImplicitEuler(; chunk_size = Val{0}(),
                autodiff = AutoForwardDiff(),
                standardtag = Val{true}(),
                concrete_jac = nothing,
                diff_type = Val{:forward}(),
                linsolve = nothing,
                precs = DEFAULT_PRECS,
                nlsolve = NLNewton(),
                extrapolant = :constant,
                controller = :PI,
                step_limiter! = trivial_limiter!)

SDIRK Method. A 1st order implicit solver. A-B-L-stable. Adaptive timestepping through a divided differences estimate. Strong-stability preserving (SSP). Good for highly stiff equations.

Keyword Arguments

autodiff: Uses ADTypes.jl to specify whether to use automatic differentiation via ForwardDiff.jl or finite differencing via FiniteDiff.jl. Defaults to AutoForwardDiff() for automatic differentiation, which by default uses chunksize = 0, and thus uses the internal ForwardDiff.jl algorithm for the choice. To use FiniteDiff.jl, the AutoFiniteDiff() ADType can be used, which has a keyword argument fdtype with default value Val{:forward}(), and alternatives Val{:central}() and Val{:complex}().
standardtag: Specifies whether to use package-specific tags instead of the ForwardDiff default function-specific tags. For more information, see this blog post. Defaults to Val{true}().
concrete_jac: Specifies whether a Jacobian should be constructed. Defaults to nothing, which means it will be chosen true/false depending on circumstances of the solver, such as whether a Krylov subspace method is used for linsolve.
linsolve: Any LinearSolve.jl compatible linear solver. For example, to use KLU.jl, specify ImplicitEuler(linsolve = KLUFactorization()). When nothing is passed, uses DefaultLinearSolver.
precs: Any LinearSolve.jl-compatible preconditioner can be used as a left or right preconditioner. Preconditioners are specified by the Pl,Pr = precs(W,du,u,p,t,newW,Plprev,Prprev,solverdata) function where the arguments are defined as:
- W: the current Jacobian of the nonlinear system. Specified as either $I - \gamma J$ or $I/\gamma - J$ depending on the algorithm. This will commonly be a WOperator type defined by OrdinaryDiffEq.jl. It is a lazy representation of the operator. Users can construct the W-matrix on demand by calling convert(AbstractMatrix,W) to receive an AbstractMatrix matching the jac_prototype.
- du: the current ODE derivative
- u: the current ODE state
- p: the ODE parameters
- t: the current ODE time
- newW: a Bool which specifies whether the W matrix has been updated since the last call to precs. It is recommended that this is checked to only update the preconditioner when newW == true.
- Plprev: the previous Pl.
- Prprev: the previous Pr.
- solverdata: Optional extra data the solvers can give to the precs function. Solver-dependent and subject to change.
The return is a tuple (Pl,Pr) of the LinearSolve.jl-compatible preconditioners. To specify one-sided preconditioning, simply return nothing for the preconditioner which is not used. Additionally, precs must supply the dispatch:
```
Pl, Pr = precs(W, du, u, p, t, ::Nothing, ::Nothing, ::Nothing, solverdata)
```
which is used in the solver setup phase to construct the integrator type with the preconditioners (Pl,Pr). The default is precs=DEFAULT_PRECS where the default preconditioner function is defined as:
```
DEFAULT_PRECS(W, du, u, p, t, newW, Plprev, Prprev, solverdata) = nothing, nothing
```
nlsolve: TBD
- extrapolant: TBD
controller: TBD
step_limiter!: function of the form limiter!(u, integrator, p, t)

References

@book{wanner1996solving, title={Solving ordinary differential equations II}, author={Wanner, Gerhard and Hairer, Ernst}, volume={375}, year={1996}, publisher={Springer Berlin Heidelberg New York}}