Backward differentiation formulas
| Bill Gear (2007), Scholarpedia, 2(8):3162. | doi:10.4249/scholarpedia.3162 | revision #91024 [link to/cite this article] |
Contents |
Backward Differentiation Methods
These are numerical integration methods based on Backward Differentiation Formulas (BDFs). They are particularly useful for stiff differential equations and Differential-Algebraic Equations (DAEs). BDFs are formulas that give an approximation to a derivative of a variable at a time \(t_n\) in terms of its function values \(y(t) \) at \(t_n\) and earlier times (hence the "backward" in the name). They are derived by forming the \(k\)-th degree interpolating polynomial approximating the function \(y(t) \) using \(y(t_n), y(t_{n-1}), \cdots, y(t_{n-k})\ ,\) differentiating it, and evaluating it at \(t_n\ .\)
For example, the linear interpolating polynomial through \(y(t_n) \) and \(y(t_{n-1})\) is
\[ y(t) \approx y(t_n) + (t - t_n)\frac{y_n - y_{n-1}}{t_n - t_{n-1}} \]
so the approximation to the derivative is the familiar
\[ y^\prime_n \approx \frac{y_n - y_{n-1}}{t_n - t_{n-1}} \]
If this is used to obtain a numerical approximation to the ordinary differential equation
\[\tag{1} y^\prime = f(y,t) \]
by replacing the derivative on the left hand side of equation
(1), one obtains the Backward Euler method
\[\tag{2} y_n = y_{n-1} + (t_n - t_{n-1})f(y_n,t_n) \]
If \(y_{n-1}\) is known, then equation (2)
is implicit in \(y_n\) --- it occurs on both sides of
the equation. (Implicitness is essential for arbitrarily
Stiff Systems.) Because equation (2) is based on a linear
approximation to \(y(t)\ ,\) it is a first-order method.
BDFs can also be used directly to solve some DAEs of the form \[\tag{3} F(y^\prime,y,t) = 0 \]
by replacing the derivative term by an expression such as in equation (2) to get \[\tag{4} F(\frac{y_n - y_{n-1}}{t_n - t_{n-1}},y_n,t_n) = 0 \]
The reader should be aware that this method will work only if the DAE has index no greater than 1 or has other special properties (see DAEs).
When higher-order methods are discussed, it is easier to discuss constant step size methods. The step size is the gap between adjacent time points, \(h_n = t_n - t_{n-1}\ ,\) and they are constant if \(h_n = h\) independent of \(n\ .\) The \(k\)-th degree polynomial interpolating \(y_n, y_{n-1}, \cdots, y_{n-k}\) can be written using backward differences as
\(\tag{5} y(t) \approx y_n + \frac{1}{h}(t-t_n)\nabla y_n + \frac{1}{2 h^2}(t-t_n)(t-t_{n-1})\nabla^2y_n + \cdots +\frac{1}{h^k k!}(t-t_n)\cdots(t-t_{n-k+1})\nabla^k\)
where \(\nabla\) is the backward difference operator
\(\nabla y_n = y_n - y_{n-1}\) and \(\nabla^p y_n =
\nabla^{p-1}y_n - \nabla^{p-1}y_{n-1}\ .\) If equation
(5) is differentiated and \(t\) is set to
\(t_n\ ,\) the \(k\)-th order BDF is obtained. The
formula has the form
\[\tag{6} hy^\prime_n = \sum_{j = 0}^k \alpha_{kj} y_{n-j}\]
where the coefficients \(\{\alpha_{kj}\}\) are
| Order k | \(\alpha_{k0}\ !\)! \(\alpha_{k1}\ !\)
\(\alpha_{k2}\ !\)! \(\alpha_{k3}\ !\)
\(\alpha_{k4}\ !\)! \(\alpha_{k5}\ !\)
\(\alpha_{k6}\)
| ||||||
|---|---|---|---|---|---|---|---|
| 1 | 1 | -1 | |||||
| 2 | 3/2 | -2 | 1/2 | ||||
| 3 | 11/6 | -3 | 3/2 | -1/3 | |||
| 4 | 25/12 | -4 | 3 | -4/3 | 1/4 | ||
| 5 | 137/60 | -5 | 5 | -10/3 | 5/4 | -1/5 | |
| 6 | 49/20 | -6 | 15/2 | -20/3 | 15/4 | -6/5 | 1/6 |
While equation (6) provides an easy way to discuss BDF's, quality codes implement a variable step size (and variable order) version of these methods, often using modified divided differences, which are an unequal step size version of the backwards differences used above. Another representation often used (with either equal or unequal step sizes) is the Nordsieck array, in which the derivatives of \(y\) at \(t_n\) are approximated by difference quotients.
The importance of the BDF-based methods is their stability. The region of absolute stability of a method is that set of \(h\lambda\) such that when the method is applied to the test equation \(y^\prime = \lambda y\) with complex \(\lambda\ ,\) the numerical solution is non-increasing in magnitude. The stability regions for the methods of order 1 through 6 are shown inFigure 1 and Figure 2. The stable regions are the exterior of the contours indicated (the full contour is plotted only for orders 1 and 2). The other contours close in the right-half plane.) Note that these six methods are stable along the whole of the negative real axis, the fact that makes them suitable for stiff equations. The contour of the BDF method of order 7 crosses the negative real axis, making it and any higher order BDF methods of no value.
The first use of BDF methods appears to date back to Curtiss and Hirschfelder (1952), although they were not given that name at the time. Later, Henrici (1962), in Section 5.1-4 of his book, discussed "methods based on differentiation." The methods were dismissed by Henrici because they are "less accurate that the corresponding Adams-Moulton formula." For non-stiff equations this is a valid point. For stiff systems, the value of BDF methods lies in their superior stability properties which allow them to take much larger stepsizes than would be possible with explicit methods.
References
- C. F. Curtiss and J. O. Hirschfelder, Integration of Stiff Equations, Proc US Nat. Acad. Sci, 38, #3 pp235-243, March, 1952
- P. Henrici, Discrete variable methods in ordinary differential equations, New York, Wiley (1962)
Internal references
- Lawrence F. Shampine and Skip Thompson (2007) Initial value problems. Scholarpedia, 2(3):2861.
- Philip Holmes and Eric T. Shea-Brown (2006) Stability. Scholarpedia, 1(10):1838.
- Lawrence F. Shampine and Skip Thompson (2007) Stiff systems. Scholarpedia, 2(3):2855.
