Backward differentiation formula
The backward differentiation formula (BDF) is a family of implicit multistep numerical methods designed for solving initial value problems in ordinary differential equations (ODEs), where the approximation of the derivative at the current time step relies on a backward difference polynomial fitted to past solution values.[1] These methods generate a sequence of approximations y_n to the exact solution y(t_n) at discrete times t_n = t_0 + n h, using the general form h y_n' = \sum_{j=0}^k \alpha_{k j} y_{n-j}, where h is the step size, k is the method order, and the coefficients \alpha_{k j} are determined by the requirement that the formula is exact for polynomials up to degree k.[1] For example, the first-order BDF (Backward Euler method) is y_n = y_{n-1} + h f(t_n, y_n), while higher orders incorporate more previous points for increased accuracy.[1] Introduced in the early 1950s as part of efforts to address stiff ODEs—systems where components evolve on vastly different time scales—BDFs were first proposed by Charles F. Curtiss and Joseph O. Hirschfelder in their work on integrating stiff equations, recognizing the need for stable implicit schemes to avoid restrictive step sizes imposed by explicit methods. The methods gained prominence through the contributions of C. William Gear, who formalized and analyzed them in the late 1960s, including variable-order and variable-step implementations that adapt to the problem's stiffness.[1] Gear's seminal book further established BDFs as a cornerstone for stiff solvers, emphasizing their stability properties for orders 1 through 6, which allow large step sizes without oscillations in stiff components.[1] A key advantage of BDFs lies in their stability region, which includes the entire negative real axis in the complex plane for orders up to 6 (A-stable for orders 1-2 and A(α)-stable for 3-6), making them particularly effective for stiff systems arising in chemical kinetics, circuit simulation, and mechanical vibrations.[1][2] However, orders 7 and higher exhibit instability along the negative real axis, limiting practical use to lower orders in most implementations.[1] BDFs are also applicable to differential-algebraic equations (DAEs) of index at most 1, where they handle algebraic constraints alongside differential components.[1] Modern solvers, such as those in scientific software, often employ variable-order BDFs with error control to balance accuracy and efficiency.[1]Overview
Definition and Purpose
The backward differentiation formula (BDF) is a family of implicit linear multistep methods designed for the numerical integration of initial value problems in ordinary differential equations (ODEs) of the form y' = f(t, y). These methods generate approximations to the solution by leveraging information from multiple previous time steps to advance the solution at the current step.[3] BDFs approximate the derivative y'(t_n) through backward differences, which arise from fitting a polynomial to the solution values at the current point t_n and the preceding k-1 points, then differentiating this interpolant and equating it to f(t_n, y_n). This interpolation-based approach ensures the method's consistency with the underlying ODE while incorporating implicitness to handle challenging dynamics.[3] The primary purpose of BDFs is to offer robust stability for stiff ODEs, where the system's eigenvalues have widely varying magnitudes, leading to rapid transients that explicit methods cannot capture without impractically small step sizes. In stiff systems, explicit Runge-Kutta or Adams methods become unstable unless the step size is restricted to h < O(1 / |\lambda_{\max}|), where \lambda_{\max} is the eigenvalue with the largest magnitude, often resulting in inefficient computation. BDFs, being implicit and A-stable (or L-stable for lower orders), permit larger steps by damping these fast components effectively, making them suitable for applications like chemical kinetics or circuit simulation.[4] To illustrate the need for such stable implicit methods, consider the prototypical stiff problem y' = -1000 y, y(0) = 1, with exact solution y(t) = e^{-1000 t}, which exhibits extremely rapid decay. Explicit methods fail to integrate this stably beyond step sizes around h \approx 0.002, whereas BDFs maintain accuracy and stability with steps up to h \approx 0.1 or larger, highlighting their utility in avoiding the stiffness-induced step-size bottleneck.[4]Historical Development
The backward differentiation formulas (BDFs) emerged as a specialized class of linear multistep methods for solving ordinary differential equations (ODEs), building on foundational work in multistep integration techniques developed in the late 19th and early 20th centuries. Early multistep methods, such as the Adams-Bashforth explicit formulas introduced by John Couch Adams in 1883 for predictor steps and extended by Forest Ray Moulton in the 1920s into implicit Adams-Moulton corrector methods, provided efficient ways to approximate solutions by leveraging multiple past points, particularly for non-stiff problems. These approaches laid the groundwork for implicit methods like BDFs, which shift focus from quadrature to direct differentiation to handle the growing recognition of stiffness in chemical kinetics and other applications. The specific introduction of BDFs is credited to Charles F. Curtiss and Joseph O. Hirschfelder in their 1952 paper, where they proposed implicit multistep formulas based on backward differences to address the challenges of stiff ODEs arising in reaction kinetics simulations. This work marked an early formal recognition of stiffness issues and advocated for methods with improved stability over explicit Runge-Kutta approaches, though the formulas were not yet termed "backward differentiation" at the time. Their contribution highlighted the need for A-stable methods capable of resolving disparate timescales in physical systems without severe step-size restrictions. BDFs gained widespread adoption through the efforts of C. William Gear, who in 1967 formalized and popularized these methods in his paper on the automatic integration of stiff ODEs, emphasizing their suitability for variable-step implementations in computational chemistry and engineering. Gear's subsequent 1971 book further detailed the algorithms, including predictor-corrector variants and stability analysis, solidifying BDFs as a cornerstone for stiff system solvers and influencing early software like DIFSUB. In the 1980s, advancements focused on enhancing BDF flexibility, with J.R. Cash introducing extended backward differentiation formulas (EBDFs) and modified extended variants (MEBDFs) to support variable-order and variable-step computations, improving efficiency for complex stiff problems without sacrificing stability up to order six.[5] These developments addressed limitations in fixed-order BDFs and paved the way for robust implementations in modern numerical libraries.Mathematical Formulation
General Formula
The general k-step backward differentiation formula (BDF) for numerically integrating the ordinary differential equation y' = f(t, y) is expressed as \sum_{j=0}^{k} \alpha_j y_{n+j} = h f(t_{n+k}, y_{n+k}), where h > 0 is the uniform step size, t_{n+j} = t_n + j h, the coefficients \{\alpha_j\}_{j=0}^k satisfy \alpha_k = 1 and are determined to ensure consistency and order k. This implicit linear multistep method is particularly suited for stiff systems due to its favorable stability properties.[1] The backward difference operator provides a foundational basis for the BDF coefficients and is defined recursively as \nabla y_n = y_n - y_{n-1} for the first-order difference, with higher-order differences given by \nabla^m y_n = \nabla (\nabla^{m-1} y_n) = \nabla^{m-1} y_n - \nabla^{m-1} y_{n-1} for m = 2, \dots, k. The k-th backward difference \nabla^k y_{n+k} approximates h^k y^{(k)}(\xi) for some \xi in the interval, linking the discrete differences to the continuous derivatives underlying the method. The BDF formula arises from polynomial interpolation: consider the unique polynomial \pi_k(t) of degree at most k that interpolates the points (t_{n+j}, y_{n+j}) for j = 0, \dots, k; the method then approximates y'(t_{n+k}) \approx \pi_k'(t_{n+k}), yielding the implicit equation above after substitution and normalization. This interpolation perspective ensures the method's order by matching the Taylor expansion up to order k. In variable step-size formulations, the BDF relates to divided differences, where the coefficients are computed using generalized divided differences f[t_{n}, t_{n-1}, \dots, t_{n-k}] \approx y^{(k)}(\xi)/k! to maintain accuracy without uniform spacing. This connection facilitates extensions to non-uniform grids while preserving the method's core approximation properties.[1]Derivation of Coefficients
The derivation of the coefficients in the backward differentiation formula (BDF) relies on polynomial interpolation to approximate the solution of the ordinary differential equation y' = f(t, y). Consider k+1 equally spaced points (t_{n+j}, y_{n+j}) for j = 0, 1, ..., k, where the step size is constant h, so t_{n+j} = t_n + j h. The interpolating polynomial π_k(t) of degree at most k is constructed to pass through these points, satisfying π_k(t_{n+j}) = y_{n+j} for each j. The BDF equation is then formed by imposing the condition that the derivative of this polynomial at the most recent point equals the right-hand side of the differential equation: π_k'(t_{n+k}) = f(t_{n+k}, y_{n+k}). This condition yields a linear relation among the y_{n+j} values and f_{n+k}, with coefficients determined by the interpolation properties. To obtain the explicit form, express the interpolating polynomial using the Lagrange basis: \pi_k(t) = \sum_{j=0}^k y_{n+j} \, l_j(t), where the basis polynomials are l_j(t) = \prod_{\substack{m=0 \\ m \neq j}}^k \frac{t - t_{n+m}}{t_{n+j} - t_{n+m}} = \prod_{\substack{m=0 \\ m \neq j}}^k \frac{t - t_n - m h}{j h - m h}. Differentiating gives \pi_k'(t) = \sum_{j=0}^k y_{n+j} \, l_j'(t). Evaluating at t = t_{n+k} produces the characteristic equation f_{n+k} = \sum_{j=0}^k y_{n+j} \, l_j'(t_{n+k}), where the coefficients are the evaluated basis derivatives l_j'(t_{n+k}). Since the points are equally spaced, these derivatives scale with 1/h, reflecting the geometric progression in the denominators. The general BDF form is thus \sum_{j=0}^k \alpha_j y_{n+j} = h f_{n+k}, with α_j = h , l_j'(t_{n+k}) for normalization such that the coefficient of f is h. This interpolation ensures the method is exact for polynomials of degree at most k.[1] The coefficients α_j can be derived using properties of the Lagrange basis under equal spacing. Standard values for low orders are as follows (with indexing shifted to match the common backward form where the sum is over y_n to y_{n-k}, but equivalent here):| Order k | α_k (current) | α_{k-1} | α_{k-2} | α_{k-3} | α_{k-4} | α_{k-5} | α_{k-6} |
|---|---|---|---|---|---|---|---|
| 1 | 1 | -1 | |||||
| 2 | 1 | -4/3 | 1/3 | ||||
| 3 | 1 | -18/11 | 9/11 | -2/11 | |||
| 4 | 1 | -48/25 | 36/25 | -16/25 | 3/25 | ||
| 5 | 1 | -300/137 | 300/137 | -200/137 | 75/137 | -12/137 | |
| 6 | 1 | -360/49 | 450/49 | -400/49 | 225/49 | -72/49 | 10/49 |