GitHub

OrdinaryDiffEqBDF

BDF (Backward Differentiation Formula) methods for fully implicit differential-algebraic equations (DAEs) in the form F(du/dt, u, t) = 0. These methods provide robust integration for index-1 DAE systems with fully implicit formulations.

Performance Consideration

DFBDF and family have not been made fully efficient yet, and thus Sundials.jl IDA is recommended for production use.

Key Properties

Fully implicit DAE BDF methods provide:

  • General DAE capability for F(du/dt, u, t) = 0 formulations
  • Index-1 DAE support for properly formulated DAE systems
  • Robust nonlinear solver integration for implicit equation systems
  • High-order accuracy with excellent stability properties
  • Large stiff system capability with efficient linear algebra

When to Use Fully Implicit DAE BDF Methods

These methods are recommended for:

  • Fully implicit DAE systems where F(du/dt, u, t) = 0 cannot be easily rearranged
  • Index-1 DAE problems that cannot be easily rearranged to semi-explicit form
  • Multibody dynamics with complex kinematic constraints
  • Electrical circuits with ideal components and algebraic loops
  • Chemical engineering with equilibrium and conservation constraints
  • Large-scale DAE systems requiring robust implicit integration

Mathematical Background

Fully implicit DAEs have the general form: F(du/dt, u, t) = 0

Unlike semi-explicit forms, these cannot be written as du/dt = f(u,t) even after constraint elimination. BDF methods discretize the time derivative using backward differences and solve the resulting nonlinear system at each timestep.

Problem Formulation

Use DAEProblem with implicit function specification:

function f2(out, du, u, p, t)
    out[1] = -0.04u[1] + 1e4 * u[2] * u[3] - du[1]
    out[2] = +0.04u[1] - 3e7 * u[2]^2 - 1e4 * u[2] * u[3] - du[2]
    out[3] = u[1] + u[2] + u[3] - 1.0
end
u₀ = [1.0, 0, 0]
du₀ = [-0.04, 0.04, 0.0]
tspan = (0.0, 100000.0)
differential_vars = [true, true, false]
prob = DAEProblem(f2, du₀, u₀, tspan, differential_vars = differential_vars)
sol = solve(prob, DFBDF())

Solver Selection Guide

  • DFBDF: Recommended - Variable-order BDF for general DAE systems
  • DImplicitEuler: For non-smooth problems with discontinuities

Method characteristics

  • DFBDF: Most robust and efficient for general smooth DAE problems
  • DImplicitEuler: Best choice for problems with discontinuities, events, or non-smooth behavior

Performance Guidelines

When fully implicit DAE BDF methods excel

  • Index-1 DAE systems with complex implicit structure
  • Complex constraint structures with multiple algebraic relationships
  • Large-scale problems where specialized DAE methods are essential
  • Multiphysics simulations with mixed differential-algebraic structure
  • Problems where semi-explicit formulation is impractical

Index considerations

  • Index-1 formulation required: Problems should be written in index-1 form
  • Compare with mass matrix methods: For some index-1 problems, mass matrix formulation may be more efficient
  • Higher-index problems: Should be reduced to index-1 form before using these methods

Important DAE Requirements

Initial conditions

  • Both u₀ and du₀ must be provided and consistent with constraints
  • differential_vars specification helps identify algebraic variables
  • Consistent initialization is crucial for index-1 DAE problems

Function specification

  • Residual form: F(du/dt, u, t) = 0 with F returning zero for satisfied equations
  • Proper scaling: Ensure equations are well-conditioned numerically
  • Jacobian availability: Analytical Jacobians improve performance when available

Alternative Approaches

Consider these alternatives:

  • Mass matrix DAE methods for index-1 problems with M du/dt = f(u,t) structure
  • Index reduction techniques using ModelingToolkit.jl to convert problems to index-1 form if needed
  • Constraint stabilization methods for drift control
  • Projection methods for manifold preservation

For more details on DAE formulations and alternative approaches, see this blog post on Neural DAEs.

Installation

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

using Pkg
Pkg.add("OrdinaryDiffEqBDF")

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 OrdinaryDiffEqBDF

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, DFBDF())

Full list of solvers

DAE

OrdinaryDiffEqBDF.DImplicitEuler — Type
DImplicitEuler(; chunk_size = Val{0}(),
                 autodiff = AutoForwardDiff(),
                 standardtag = Val{true}(),
                 concrete_jac = nothing,
                 linsolve = nothing,
                 precs = DEFAULT_PRECS,
                 nlsolve = NLNewton(),
                 extrapolant = :constant,
                 controller = :Standard)

Multistep Method. 1st order A-L and stiffly stable adaptive implicit Euler. Implicit Euler for implicit DAE form. It uses an apriori error estimator for adaptivity based on a finite differencing approximation from SPICE.

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 DImplicitEuler(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
    /n- nlsolve: TBD
  • extrapolant: TBD
  • controller: TBD

References

source
OrdinaryDiffEqBDF.DABDF2 — Type
DABDF2(; chunk_size = Val{0}(),
         autodiff = AutoForwardDiff(),
         standardtag = Val{true}(),
         concrete_jac = nothing,
         linsolve = nothing,
         precs = DEFAULT_PRECS,
         nlsolve = NLNewton(),
         extrapolant = :constant,
         controller = :Standard)

Multistep Method. 2nd order A-L stable adaptive BDF method. Fully implicit implementation of BDF2.

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 DABDF2(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
    /n- nlsolve: TBD
  • extrapolant: TBD
  • controller: TBD

References

@article{celaya2014implementation, title={Implementation of an Adaptive BDF2 Formula and Comparison with the MATLAB Ode15s}, author={Celaya, E Alberdi and Aguirrezabala, JJ Anza and Chatzipantelidis, Panagiotis}, journal={Procedia Computer Science}, volume={29}, pages={1014–1026}, year={2014}, publisher={Elsevier}}

source
OrdinaryDiffEqBDF.DFBDF — Type
DFBDF(; chunk_size = Val{0}(),
        autodiff = AutoForwardDiff(),
        standardtag = Val{true}(),
        concrete_jac = nothing,
        linsolve = nothing,
        precs = DEFAULT_PRECS,
        κ = nothing,
        tol = nothing,
        nlsolve = NLNewton(),
        extrapolant = :linear,
        controller = :Standard,
        max_order::Val{MO} = Val{5}())

Multistep Method. Fixed-leading coefficient adaptive-order adaptive-time BDF method. Fully implicit implementation of FBDF based on Shampine's

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 DFBDF(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
    /n- κ: TBD
  • tol: TBD
  • nlsolve: TBD
  • extrapolant: TBD
  • controller: TBD
  • max_order: TBD

References

@article{shampine2002solving, title={Solving 0= F (t, y (t), y′(t)) in Matlab}, author={Shampine, Lawrence F}, year={2002}, publisher={Walter de Gruyter GmbH and Co. KG} }

source