OrdinaryDiffEqBDF

Backward Differentiation Formula (BDF) methods are multistep implicit methods specifically designed for solving large stiff systems of differential equations. They are the preferred choice for very large systems (>1000 equations) where other implicit methods become computationally expensive.

Key Properties

BDF methods offer:

Excellent efficiency for large systems (>1000 ODEs)
L-stable behavior for orders 1 and 2 only
Adaptive order and stepsize control for optimal performance
Alpha-stability for higher orders (but less stable than L-stable methods for problems with large complex eigenvalues)

When to Use BDF Methods

BDF methods are recommended for:

Large stiff systems with more than 1000 equations
Very stiff problems where other implicit methods struggle
Long-time integration of stiff systems
Parabolic PDEs after spatial discretization
Reaction-diffusion systems and chemical kinetics
Circuit simulation and other engineering applications with large stiff systems

Solver Selection Guide

Recommended methods

QNDF: Adaptive order quasi-constant timestep BDF, best general choice for large systems
FBDF: Fixed-leading coefficient BDF, often more efficient than QNDF

Performance Characteristics

Most efficient for systems with >1000 equations
Outperform Runge-Kutta methods on very large stiff systems
Memory efficient due to multistep structure
Excel at very low tolerances (1e-9 and below)
Particularly effective for problems arising from PDE discretizations

Comparison with Other Methods

Choose BDF methods over:

Rosenbrock methods: When system size > 1000 equations
SDIRK methods: For very large stiff systems where RK methods become expensive
Explicit methods: For any stiff problem

Choose other methods over BDF when:

System size < 100: Rosenbrock or SDIRK methods often more efficient
Problems with large complex eigenvalues: Rosenbrock and L-stable SDIRK methods are more stable due to BDF methods only being alpha-stable
Moderate stiffness: SDIRK methods may be more robust
Non-stiff problems: Use explicit methods like Tsit5

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

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

OrdinaryDiffEqBDF.NordsieckBDF — Type

NordsieckBDF(; nlsolve = NLNewton(max_iter = 3),
               extrapolant = :linear,
               max_order::Val{MO} = Val{5}(),
               step_limiter! = trivial_limiter!)

Multistep Method An adaptive-order, adaptive-time BDF method on a propagated Nordsieck history array zn[j] = h^j/j! * y^(j)(t_n), following SUNDIALS CVODE.

Where FBDF stores raw (t_i, u_i) history and rebuilds the predictor and the error/order estimates from it every step, this method propagates one array: predicting is a Pascal-triangle shift, changing the step size is zn[j] *= eta^j, and accepting a step is a rank-1 update. Because nothing is reconstructed from stored points, a loose nonlinear solve cannot be amplified into a step-size collapse, so the corrector can be solved to a fraction of the local error budget (nlsolve = NLNewton(κ = …), CVODE's NLSCOEF) instead of to full accuracy. On stiff benchmarks that is worth roughly a factor of two in f-evaluations relative to FBDF.

Dense output is the Nordsieck polynomial itself and is free.

Keyword Arguments

nlsolve: nonlinear solver for the implicit stage. Its κ acts as CVODE's NLSCOEF, i.e. the fraction of the local error budget the corrector is allowed to consume, because the increment norm is scaled by the test quantity tq[2]. The default caps the corrector at 3 iterations, matching CVODE's MAXCOR: past that it is cheaper to give up, refresh W, and re-converge than to keep iterating with a stale Jacobian — a trade this method can take because refactorizing is comparatively rare for it.
max_order: maximum BDF order (1–5).
step_limiter!: function of the form limiter!(u, integrator, p, t).

References

@article{byrne1975polyalgorithm, title={A polyalgorithm for the numerical solution of ordinary differential equations}, author={Byrne, George D and Hindmarsh, Alan C}, journal={ACM Transactions on Mathematical Software}, volume={1}, number={1}, pages={71–96}, year={1975}} @article{hindmarsh2005sundials, title={{SUNDIALS}: Suite of nonlinear and differential/algebraic equation solvers}, author={Hindmarsh, Alan C and Brown, Peter N and Grant, Keith E and Lee, Steven L and Serban, Radu and Shumaker, Dan E and Woodward, Carol S}, journal={ACM Transactions on Mathematical Software}, volume={31}, number={3}, pages={363–396}, year={2005}}

source

OrdinaryDiffEqBDF.DNordsieckBDF — Type

DNordsieckBDF(; nlsolve = NLNewton(max_iter = 3),
                extrapolant = :linear,
                max_order::Val{MO} = Val{5}())

Fully Implicit Multistep Method Fully implicit DAE solver: the NordsieckBDF method applied to f(du, u, p, t) = 0. The Nordsieck array supplies both the state predictor and the derivative du = zn[1]/h, so the corrector solves f((zn[1] + l[1]*acor)/h, ypred + acor, p, t) = 0 with leading coefficient cj = l[1]/h — the same role IDA's cj plays.

Keyword Arguments

nlsolve: nonlinear solver for the implicit stage; its κ acts as NLSCOEF.
max_order: maximum BDF order (1–5).

References

@article{hindmarsh2005sundials, title={{SUNDIALS}: Suite of nonlinear and differential/algebraic equation solvers}, author={Hindmarsh, Alan C and Brown, Peter N and Grant, Keith E and Lee, Steven L and Serban, Radu and Shumaker, Dan E and Woodward, Carol S}, journal={ACM Transactions on Mathematical Software}, volume={31}, number={3}, pages={363–396}, year={2005}}

source