OrdinaryDiffEqBDF

BDF (Backward Differentiation Formula) methods for mass matrix differential-algebraic equations (DAEs) and stiff ODEs with singular mass matrices. These methods provide robust, high-order integration for systems with algebraic constraints and mixed differential-algebraic structure.

Key Properties

Mass matrix BDF methods provide:

DAE capability for index-1 differential-algebraic equations
Mass matrix support for singular and non-diagonal mass matrices
High-order accuracy up to 5th order with good stability
L-stable behavior for stiff problems with excellent damping
Automatic differentiation for efficient Jacobian computation
Variable order and stepsize adaptation for efficiency

When to Use Mass Matrix BDF Methods

These methods are recommended for:

Differential-algebraic equations (DAEs) with index-1 structure
Constrained mechanical systems with holonomic constraints
Electrical circuit simulation with algebraic loop equations
Chemical reaction networks with conservation constraints
Multibody dynamics with kinematic constraints
Semi-explicit DAEs arising from spatial discretizations

Mathematical Background

Mass matrix DAEs have the form: M du/dt = f(u,t)

where M is a potentially singular mass matrix. When M is singular, some equations become algebraic constraints rather than differential equations, leading to a DAE system.

Problem Formulation

Use ODEFunction with a mass_matrix:

using LinearAlgebra: Diagonal
function rober(du, u, p, t)
    y₁, y₂, y₃ = u
    k₁, k₂, k₃ = p
    du[1] = -k₁ * y₁ + k₃ * y₂ * y₃
    du[2] = k₁ * y₁ - k₃ * y₂ * y₃ - k₂ * y₂^2
    du[3] = y₁ + y₂ + y₃ - 1
    nothing
end
M = Diagonal([1.0, 1.0, 0])  # Singular mass matrix
f = ODEFunction(rober, mass_matrix = M)
prob_mm = ODEProblem(f, [1.0, 0.0, 0.0], (0.0, 1e5), (0.04, 3e7, 1e4))
sol = solve(prob_mm, FBDF(), reltol = 1e-8, abstol = 1e-8)

Solver Selection Guide

Recommended Methods

FBDF: Recommended - Fixed leading coefficient BDF with excellent stability
QNDF: Quasi-constant stepsize Nordsieck BDF with good efficiency
QBDF: Alternative quasi-constant stepsize BDF formulation

Specific order methods

QNDF1: First-order method for simple problems
QNDF2: Second-order method balancing accuracy and stability
QBDF1, QBDF2: Alternative second-order formulations
ABDF2: Adams-type BDF for specific applications
MEBDF2: Modified extended BDF for enhanced stability

Performance Guidelines

When mass matrix BDF methods excel

Index-1 DAE systems with well-separated differential and algebraic variables
Large stiff systems with algebraic constraints
Problems with conservation laws naturally expressed as constraints
Multiphysics simulations combining differential and algebraic equations
Systems where constraints are essential to the physics

Mass matrix considerations

Singular mass matrices require consistent initial conditions
Index determination affects solver performance and stability
Constraint violations may accumulate and require projection
Well-conditioned problems generally perform better

Important Considerations

Initial conditions

Must be consistent with algebraic constraints
Use initialization procedures if constraints are not satisfied initially
Index-1 assumption requires that constraints uniquely determine algebraic variables

Numerical challenges

Constraint drift may occur over long integrations
Index higher than 1 not directly supported
Ill-conditioned mass matrices can cause numerical difficulties
Discontinuities in constraints require special handling

Alternative Approaches

Consider these alternatives:

Implicit Runge-Kutta methods for higher accuracy requirements
Rosenbrock methods for moderately stiff DAEs
Projection methods for constraint preservation
Index reduction techniques for higher-index 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, FBDF())

Full list of solvers

OrdinaryDiffEqBDF.ABDF2 — Type

ABDF2(; autodiff = AutoForwardDiff(),
        concrete_jac = nothing,
        linsolve = nothing,
        κ = nothing,
        tol = nothing,
        nlsolve = NLNewton(),
        smooth_est = true,
        extrapolant = :linear,
        step_limiter! = trivial_limiter!)

Multistep Method. An adaptive order 2 L-stable fixed leading coefficient multistep BDF method.

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}().
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 ABDF2(linsolve = KLUFactorization()). When nothing is passed, uses DefaultLinearSolver. /n- κ: coefficient for the order and stability control of the BDF method. When nothing, the default value is used.
tol: tolerance for the nonlinear solver. When nothing, uses the default tolerance.
nlsolve: nonlinear solver algorithm used for solving the implicit system.
smooth_est: whether to use a smoothed estimate for error control.
extrapolant: extrapolation method used for the initial guess in the nonlinear solve.
step_limiter!: function of the form limiter!(u, integrator, p, t)

References

E. Alberdi Celayaa, J. J. Anza Aguirrezabalab, P. Chatzipantelidisc. Implementation of an Adaptive BDF2 Formula and Comparison with The MATLAB Ode15s. Procedia Computer Science, 29, pp 1014-1026, 2014. doi: https://doi.org/10.1016/j.procs.2014.05.091

source

OrdinaryDiffEqBDF.QNDF — Type

QNDF(; autodiff = AutoForwardDiff(),
       concrete_jac = nothing,
       linsolve = nothing,
       κ = nothing,
       tol = nothing,
       nlsolve = NLNewton(),
       extrapolant = :linear,
       kappa =  promote(-0.1850, -1 // 9, -0.0823, -0.0415, 0),
       step_limiter! = trivial_limiter!)

Multistep Method. An adaptive order quasi-constant timestep NDF method. Similar to MATLAB's ode15s. Uses Shampine's accuracy-optimal coefficients. Performance improves with larger, more complex ODEs. Good for medium to highly stiff problems. Recommended for large systems (>1000 ODEs).

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}().
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 QNDF(linsolve = KLUFactorization()). When nothing is passed, uses DefaultLinearSolver. /n- κ: coefficient for the order and stability control of the BDF method. When nothing, the default value is used.
tol: tolerance for the nonlinear solver. When nothing, uses the default tolerance.
nlsolve: nonlinear solver algorithm used for solving the implicit system.
extrapolant: extrapolation method used for the initial guess in the nonlinear solve.
kappa: coefficient for the BDF error estimator.
step_limiter!: function of the form limiter!(u, integrator, p, t)

References

@article{shampine1997matlab, title={The matlab ode suite}, author={Shampine, Lawrence F and Reichelt, Mark W}, journal={SIAM journal on scientific computing}, volume={18}, number={1}, pages={1–22}, year={1997}, publisher={SIAM} }

source

OrdinaryDiffEqBDF.QNDF1 — Type

QNDF1(; autodiff = AutoForwardDiff(),
        concrete_jac = nothing,
        linsolve = nothing,
        nlsolve = NLNewton(),
        extrapolant = :linear,
        kappa = -0.1850,
        step_limiter! = trivial_limiter!)

Multistep Method. An adaptive order 1 quasi-constant timestep L-stable numerical differentiation function method.

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}().
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 QNDF1(linsolve = KLUFactorization()). When nothing is passed, uses DefaultLinearSolver. /n- nlsolve: nonlinear solver algorithm used for solving the implicit system.
extrapolant: extrapolation method used for the initial guess in the nonlinear solve.
kappa: coefficient for the BDF error estimator.
step_limiter!: function of the form limiter!(u, integrator, p, t)

References

source

OrdinaryDiffEqBDF.QNDF2 — Type

QNDF2(; autodiff = AutoForwardDiff(),
        concrete_jac = nothing,
        linsolve = nothing,
        nlsolve = NLNewton(),
        extrapolant = :linear,
        kappa =  -1 // 9,
        step_limiter! = trivial_limiter!)

Multistep Method. An adaptive order 2 quasi-constant timestep L-stable numerical differentiation function (NDF) method.

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}().
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 QNDF2(linsolve = KLUFactorization()). When nothing is passed, uses DefaultLinearSolver. /n- nlsolve: nonlinear solver algorithm used for solving the implicit system.
extrapolant: extrapolation method used for the initial guess in the nonlinear solve.
kappa: coefficient for the BDF error estimator.
step_limiter!: function of the form limiter!(u, integrator, p, t)

References

source

OrdinaryDiffEqBDF.QBDF — Function

QBDF: Multistep Method

An alias of QNDF with κ=0.

source

OrdinaryDiffEqBDF.QBDF1 — Function

QBDF1: Multistep Method

An alias of QNDF1 with κ=0.

source

OrdinaryDiffEqBDF.QBDF2 — Function

QBDF2: Multistep Method

An alias of QNDF2 with κ=0.

source

OrdinaryDiffEqBDF.MEBDF2 — Type

MEBDF2(; autodiff = AutoForwardDiff(),
         concrete_jac = nothing,
         linsolve = nothing,
         nlsolve = NLNewton(),
         extrapolant = :constant)

Multistep Method. The second order Modified Extended BDF method, which has improved stability properties over the standard BDF. Fixed timestep only.

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}().
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 MEBDF2(linsolve = KLUFactorization()). When nothing is passed, uses DefaultLinearSolver. /n- nlsolve: nonlinear solver algorithm used for solving the implicit system.
extrapolant: extrapolation method used for the initial guess in the nonlinear solve.

References

@article{cash2000modified, title={Modified extended backward differentiation formulae for the numerical solution of stiff initial value problems in ODEs and DAEs}, author={Cash, JR}, journal={Journal of Computational and Applied Mathematics}, volume={125}, number={1-2}, pages={117–130}, year={2000}, publisher={Elsevier}}

source

OrdinaryDiffEqBDF.FBDF — Type

FBDF(; autodiff = AutoForwardDiff(),
       concrete_jac = nothing,
       linsolve = nothing,
       κ = nothing,
       tol = nothing,
       nlsolve = NLNewton(),
       extrapolant = :linear,
       step_limiter! = trivial_limiter!,
       max_order::Val{MO} = Val{5}(),
       stald = true,
       stald_rrcut = 0.98,
       stald_vrrtol = 1e-4,
       stald_vrrt2 = 5e-4,
       stald_sqtol = 1e-3,
       stald_rrtol = 1e-2,
       stald_tiny = 1e-90)

Multistep Method. An adaptive order quasi-constant timestep NDF method. Fixed leading coefficient BDF. Utilizes Shampine's accuracy-optimal kappa values as defaults (has a keyword argument for a tuple of kappa coefficients).

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}().
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 FBDF(linsolve = KLUFactorization()). When nothing is passed, uses DefaultLinearSolver. /n- κ: coefficient for the order and stability control of the BDF method. When nothing, the default value is used.
tol: tolerance for the nonlinear solver. When nothing, uses the default tolerance.
nlsolve: nonlinear solver algorithm used for solving the implicit system.
extrapolant: extrapolation method used for the initial guess in the nonlinear solve.
step_limiter!: function of the form limiter!(u, integrator, p, t)
max_order: maximum order of the adaptive-order BDF method.
stald: Enable Stability Limit Detection (STALD) for BDF orders 3-5. Default: true.
stald_rrcut: STALD cutoff for characteristic root magnitude. Default: 0.98.
stald_vrrtol: STALD tolerance for variance of ratios. Default: 1e-4.
stald_vrrt2: STALD secondary variance tolerance. Default: 5e-4.
stald_sqtol: STALD tolerance for quartic residual. Default: 1e-3.
stald_rrtol: STALD tolerance for rr cross-verification. Default: 1e-2.
stald_tiny: STALD tiny value to avoid division by zero. Default: 1e-90.

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 \& Co. KG}}

source