One of the key features of DifferentialEquations.jl is its flexibility. Keeping with this trend, many of the native Julia solvers provided by DifferentialEquations.jl allow you to choose the method for linear and nonlinear solving. This section details how to make that choice.
We highly recommend looking at the Solving Large Stiff Equations tutorial which goes through these options in a real-world example.
For linear solvers, DifferentialEquations.jl uses LinearSolve.jl. Any LinearSolve.jl algorithm can be used as the linear solver simply by passing the algorithm choice to linsolve. For example, the following tells
TRBDF2 to use KLU.jl
TRBDF2(linsolve = KLUFactorization())
Many choices exist, including GPU offloading, so consult the LinearSolve.jl documentation for more details on the choices.
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
WOperatortype 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
du: the current ODE derivative
u: the current ODE state
p: the ODE parameters
t: the current ODE time
Boolwhich specifies whether the
Wmatrix 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
Prprev: the previous
solverdata: Optional extra data the solvers can give to the
precsfunction. 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.
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 in order to construct the integrator type with the preconditioners
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
All of the Julia-based implicit solvers (OrdinaryDiffEq.jl, StochasticDiffEq.jl, etc.) allow for choosing the nonlinear solver that is used to handle the implicit system. While fully modifiable and customizable, most users should stick to the pre-defined nonlinear solver choices. These are:
NLNewton(; κ=1//100, max_iter=10, fast_convergence_cutoff=1//5, new_W_dt_cutoff=1//5): A quasi-Newton method. The default.
NLAnderson(; κ=1//100, max_iter=10, max_history::Int=5, aa_start::Int=1, droptol=nothing, fast_convergence_cutoff=1//5): Anderson acceleration. While more stable than functional iteration, this method is less stable than Newton's method but does not require a Jacobian.
NLFunctional(; κ=1//100, max_iter=10, fast_convergence_cutoff=1//5): This method is the least stable but does not require Jacobians. Should only be used for non-stiff ODEs.