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

Recommended DAE Methods

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(; autodiff = AutoForwardDiff(),
                 concrete_jac = nothing,
                 linsolve = nothing,
                 nlsolve = NLNewton(),
                 extrapolant = :constant)

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}().
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. /n- nlsolve: nonlinear solver algorithm used for solving the implicit system.
extrapolant: extrapolation method used for the initial guess in the nonlinear solve.

References

source

OrdinaryDiffEqBDF.DABDF2 — Type

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

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}().
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. /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{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(; autodiff = AutoForwardDiff(),
        concrete_jac = nothing,
        linsolve = nothing,
        κ = nothing,
        tol = nothing,
        nlsolve = NLNewton(),
        extrapolant = :linear,
        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}().
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. /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.
max_order: maximum order of the adaptive-order BDF method.

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