Derivation of BDF in "fixed-leading coefficient" form

The fixed-leading coefficient form of BDF, the choice in all SUNDIALS integrators can be derived (and implemented) most conveniently via interpolating polynomials. The derivation for BDFq can be obtained in the following 10 steps.

Suppose we have computed $y_j \approx y(t_j), j = 1,2,....,n-1$. Let the past step sizes be $h_j = t_j - t_{j-1}, j=1,2,...,n-1$ and $h_n = t_n - t_{n-1}$ be the current step size. We wish to compute $y_n \approx y(t_n)$.

Step 1
Build a predictor polynomial $\omega_p$ of order $q$ that satisfies the following $(q+1)$ conditions:

(1)

\begin{eqnarray} \omega_p(t_{n-i}) &=& y_{n-i}, \quad i=1,...,q \\ \dot \omega_p(t_{n-1}) &=& \dot y_{n-1} \end{eqnarray}

In other words, $\omega_p$ interpolates the past $q$ computed values plus the derivative at the last computed step.

Step 2
Build a corrector polynomial $\omega_c$ of order $q$ that satisfies the following $(q+1)$ conditions:

(2)

\begin{eqnarray} \omega_c(t_n) &=& y_n \\ \omega_c(t_n - i h_n) &=& \omega_p(t_n - i h_n), \quad i=1,...,q \end{eqnarray}

In other words, $\omega_c$ passes through the (unknown at this point) solution at $t_n$ and coincides with the predictor polynomial at $q$ past equidistant points (separated by the current step size $h_n$).

Step 3
It is easy to see that the $q$-th order FLC BDF formula can be written simply as

(3)

\begin{align} \dot y_n = \omega_c(t_n) \end{align}

Step 4
Since both $\omega_p$ and $\omega_c$ are $q$-th order polynomial, let's relate them through a new polynomial $c(t)$ (which we will define later) as

(4)

\begin{align} \omega_c(t) = \omega_p(t) + c(t) [ \omega_c(t_n) - \omega_p(t_n) ] \end{align}

It's obvious that, once plugged into the ODE, the above BDF formula will lead to a nonlinear equation in the desired solution $y_n$. We will use some form of iteration to solve this nonlinear equation for $y$ at $t_n$ (recall that this is, by construction, $\omega_c(t_n)$ ) and we will use the predicted value at $t_n$ as our initial guess.

Step 5
Denote then

(5)

\begin{eqnarray} y_{n(0)} = \omega_p(t_n) \\ \dot y_{n(0)} = \dot \omega_p(t_n) \end{eqnarray}

Note that $y_{n(m)}, m=1.2,...$ will be successive iterates of the nonlinear solver (and $y_n$ will be $y_{n(M)}$ where $M$ is the iteration at which we declare convergence).

Step 6
With the above definitions we have

(6)

\begin{align} \omega_c(t) = \omega_p(t) + c(t) [ y_n - y_{n(0)} ] \end{align}

Step 7
Now use the relations (1) and (2) to derive conditions on $c(t)$. In other words, evaluate (6) at $t_n - i h_n, i=0,1,...,q$ and use (1) and (2) to get:

(7)

\begin{eqnarray} &&c(t_n) = 1 \\ &&c(t_n - h_n) = c(t_n - 2 h_n) = ... = c(t_n - q h_n) = 0 \end{eqnarray}

which are $(q+1)$ conditions uniquely defining the $q$-th order polynomial $c(t)$.

Step 8
Using the relation (6), the BDF formula in (3) can now be written as:

(8)

\begin{align} \dot y_n = \dot \omega_p(t_n) + \dot c(t_n) [ y_n - y_{n(0)} ] \end{align}

or, using (5)

(9)

\begin{align} \dot y_n = \dot y_{n(0)} + \dot c(t_n) [ y_n - y_{n(0)} ] \end{align}

All we need to do now to complete the derivation of the BDF formula is find the derivative of $c(t)$ at $t_n$.

Step 9
Given the conditions (7), i.e. the roots of the $q$-th order polynomial $c(t)$, we have:

(10)

\begin{align} c(t) = K \Pi_{i=1}^q { [ t - ( t_n - i h_n ) ] } \end{align}

And at $t=t_n$, we get

(11)

\begin{align} 1 = c(t_n) = K \cdot h_n \cdot (2 h_n) \cdot ... \cdot (q h_n) \end{align}

Therefore $K = {1}/{q!h_n^q}$ and $c(t)$ is completely determined.

Step 10
Differentiate (10) and evaluate at $t=t_n$ to get

(12)

\begin{align} \dot c(t_n) = \frac{1}{h_n} [ 1/1 + 1/2 + ... + 1/q ] \end{align}

or

(13)

\begin{align} h_n \dot c(t_n) = \frac{1}{\beta_0} \end{align}

where we define $\beta_0 = 1/(1+1/2+...+1/q)$.

All done. Multiply (9) by $h_n$ and use (13) to get the BDF formula as:

(14)

\begin{align} h_n \dot y_n = h_n \dot y_{n(0)} + \frac{1}{\beta_0} [ y_n - y_{n(0)} ] \end{align}

or

(15)

\begin{align} y_n = h_n \beta_0 \dot y_n + [ y_{n(0)} - h_n \beta_0 \dot y_{n(0)} ] \end{align}

or

(16)

\begin{align} y_n - y_{n(0)} = h_n \beta_0 [\dot y_n - \dot y_{n(0)} ] \end{align}