BDF Methods

CVODE, CVODES, and IDA all use Backward Differentiation Formulas methods (CVODE and CVODES also provide the option of using Adams-Moulton methods for the solution of nonstiff problems).

Basic formulas

Denote by $t_n$ and $y_n = y(t_n)$ the discrete values of the independent variable $t$ and solution $y(t)$ and let $h = t_n - t_{n-1}$ be the stepsize.

BDF is one of many Linear Multistep Methods (LMM). It is called BDF because it differentiates the solution $y$ using past (backward) values; i.e.,

(1)

\begin{align} \dot y_n = (\text{linear combination of } y_n , y_{n-1}, y_{n-2}, \ldots ) \end{align}

The simplest BDF method (BDF1) is Backward Euler (BE):

(2)

\begin{align} y_n = y_{n-1} + h \dot y_n \end{align}

The fixed-step BDF of order $q$ (BDFq) is

(3)

\begin{align} y_n = \sum_{i=1}^q \alpha_i y_{n-i} + h \beta_0 \dot y_n \end{align}

where the coefficients $\alpha_i$ and $\beta_0$ depend only on the method order $q$.

For an ODE in explicit form; i.e. $\dot y = f(t,y)$, the BDFq method gives $y_n$ as the solution of

(4)

\begin{align} y_n = \sum_{i=1}^q \alpha_i y_{n-i} + h \beta_0 f(t_n , y_n) \end{align}

For an ODE in implicit form or a DAE; i.e. $F(t,y,\dot y) = 0$, the BDFq method gives $y_n$ as the solution of

(5)

\begin{align} F \left( t_n , y_n , \frac{y_n - \sum_{i=1}^q \alpha_i y_{n-i} }{h\beta_0} \right) = 0 \end{align}

Symbolic operator derivation

For any sequence $\{x_n\}$, define the operators

(6)

\begin{eqnarray} 1 x_n = x_n && \text{(identity)} \\ E x_n = x_{n+1} && \text{(increment)} \\ \Delta x_n = x_{n+1} - x_n && \text{(forward difference)} \\ \nabla x_n = x_{n} - x_{n-1} && \text{(backward difference)} \\ \end{eqnarray}

It is easy to develop an algebra of operators by observing that the following relationships hold:

(7)

\begin{eqnarray} \Delta &=& E - 1 \\ \nabla &=& 1 - E^{-1} \end{eqnarray}

where the inverse increment operator is defined by

(8)

\begin{equation} E^{-1} x_n = x_{n-1} \end{equation}

A given sequence may or may not be the set of discrete values of a smooth functino of $t$. When it is, we can define a differentiation operator, $Df_n = \dot f(t_n)$ and write the infinite Taylor series expansion of $f(t)$ at $t = t_n$ as

(9)

\begin{eqnarray} Ef_n &=& f_{n+1} = f(t_{n+1}) = f_n + h \dot f_n + \frac{h^2}{2} \ddot f + \ldots \\ &=& f_n + hDf_n + \frac{h^2D^2}{2} f_n + \ldots \\ &=& \left( 1+ hD + \frac{h^2D^2}{2} + \ldots \right) f_n \\ &=& e^{hD} f_n \end{eqnarray}

In other words, the increment and differentiation operators are related through

(10)

\begin{equation} E = e^{hD} \end{equation}

Next, using (7), we get $E = (1-\nabla)^{-1}$ and therefore

(11)

\begin{align} hD = \ln (E) = \ln \left[ (1-\nabla)^{-1} \right] = \nabla + \frac{1}{2} \nabla^2 + \frac{1}{3} \nabla^3 + \ldots \end{align}

which, applied to $y(t)$ gives:

(12)

\begin{align} h \dot y_n = \nabla y_n + \frac{1}{2} \nabla^2 y_n + \frac{1}{3} \nabla^3 y_n + \ldots \end{align}

to obtain the BDFq method, simply truncate the series (12) at the $\nabla^q$ term. Through simple manipulations, it is easy to obtain

(13)

\begin{align} \beta_0 = \frac{1}{\sum_{j=1}^q j^{-1}} \end{align}

Varying the stepsize

There are three possible ways of allowing a variable stepsize: interpolated fixed-step BDF, fully variable-step BDF, and fixed-leading coefficient BDF.

(1) In the interpolated fixed-step approach, one uses the fixed-step BDF method with stepsize $h$, but interpolates at a spacing equal to the new stepsize $h'$. However, the behavior of the resulting method can be unstable if the stepsize changes too frequently.

(2) A fully variable-step BDF can be obtained, by allowing the method coefficients to also depend on the stepsize $h$, as

(14)

\begin{align} y_n = \sum_{i=1}^q \alpha_{n,i} y_{n-i} + h_n \beta_{n,0} \dot y_n \end{align}

Unlike the method (3), the coefficients $\alpha_{n,i}$ and $\beta_{n,0}$ depend on the method order $q$ and the last $q$ stepsizes $h_n , h_{n-1} , \ldots , h_{n-q+1}$. Denoting $h = \max \{ h_n , h_{n-1} , \ldots , h_{n-q+1} \}$, the method coefficients are uniquely determined by

(15)

\begin{align} \frac{1}{\beta_{n,0}} \left[ y(t_n) - \sum_{i=1}^q \alpha_{n,i} y(t_{n-i}) \right] = h_n \dot y (t_n) + O(h^{q+1}) \end{align}

The disadvantage of this approach is that the variation of the quantity $h_n \beta_{n,0}$ makes it harder to keep the Newton iteration matrix (for the solution of the resulting nonlinear system) current.

(3) The fixed-leading coefficient form of BDF is a compromise between the above two choices in that it has the fixed-step value (13) of $\beta_{n,0} \equiv \beta_0 = 1 / \sum_{j=1}^q j^{-1}$ but at the cost of introducing an additional term $\beta_{n,1}} \dot y_{n-1}$:

(16)

\begin{align} y_n = \sum_{i=1}^q \alpha_{n,i} y_{n-i} + h_n \beta_0 \dot y_n + h_n \beta_{n,1} \dot y_{n-1} \end{align}

As expected, the coefficient $\beta_{n,1}$ vanishes when the stepsizes are equal. The fixed-leading coefficient form of BDF is most conveniently derived via interpolating polynomials. This derivation is explained in details in a separate article.
The main advantages of FLC BDF is that it does not suffer from the unstable behavior of the interpolated fixed-step method and, since the only variation in the quantity $h_n \beta_0$ is now through the stepsize $h_n$ only, the Newton iteration matrix can be reused for more steps than in a fully variable-step approach.

The fixed-leading coefficient form of BDF is the choice in all SUNDIALS integrators (CVODE, CVODES, IDA).

Example: BDF2

Fixed-step

(17)

\begin{eqnarray} && h \dot y_n = \nabla y_n + \frac{1}{2} \nabla^2 y_n = y_n - y_{n-1} + \frac{1}{2} (y_n -2 y_{n-1} + y_{n-2}) \\ \text{or}&&\\ && y_n = \frac{4}{3} y_{n-1} - \frac{1}{3} y_{n-2} + \frac{2}{3} h \dot y_n \end{eqnarray}

Variable-step

(18)

\begin{eqnarray} &&y_n = \frac{1}{3\rho + 1} \left[ (\rho+1)^2 y_{n-1} - \rho^2 y_{n-2} + (\rho+1) h \dot y_n \right] \\ \text{where } && \rho = h_n / h_{n-1} \end{eqnarray}

Fixed-leading coefficient

(19)

\begin{eqnarray} &&y_n = \left( 1 + \frac{\rho^2}{3} \right) y_{n-1} - \frac{\rho^2}{3} y_{n-2} + \frac{2}{3} h \dot y_n - \left( \frac{\rho-1}{3} \right) h_n \dot y_{n-1} \\ \text{where } && \rho = h_n / h_{n-1} \end{eqnarray}