Backward Differentiation Formula

Starting from an unknown continuous function of time \(\mathbf q(T)\), the solution can be found numerically at discretized time points \(T_n\) with \(n\in \mathbb N\) where, to simplify the notations \(\mathbf q_n = \mathbf q(T_n)\). The goal is to be able to compute \(\mathbf q_{n+1}\) from previous known discrete values. The BDF way to proceed is to assume that the real solution can be interpolated using Lagrange polynomials.

Let be \(\mathbf{q}\) be an unknown continuous function of time \(T\) and \(\mathbf{\dot q}\) be its first derivative \[ \begin{aligned} \mathbf{q} &= \mathbf{q}(T)\\ \mathbf{\dot q} &= \frac{\mathrm{d} \mathbf{q}}{\mathrm{d}T}\\ \end{aligned} \] The solution can be found numerically at discretized time points \(T_n\) with \(n\in \mathbb N\) where, to simplify the notations, \(\mathbf q_n = \mathbf q(T_n)\). The continuous solution can be interpolated from previous discrete values, e.g. using Lagrange polynomials. \[ \mathbf{q}^h(T) = \mathbb L^p(T,\mathbf q_{n-p},\ldots,\mathbf q_n) \] Which have the property that \(\mathbf q^h(T_n)=\mathbf q_n\). Using the interpolation function instead of the unknown field \(\mathbf q\) yields an ordinary differential equation where \[ \mathbf{\dot q} = \frac{\mathrm{d} \mathbf{q}^h}{\mathrm{d}T} = \frac{\mathrm{d} \mathbb L^{p}}{\mathrm{d}T} \] Evaluating the derivative at a new location \(T_{n+1}\) yields an implicit equation \[ \mathbf{\dot q}_{n+1} = \frac{\mathrm{d}\mathbb L^p}{\mathrm dT}(T,\mathbf q_{n+1-p},\ldots,\mathbf q_{n+1}) \]

Backward Differentiation Formula is a family of implicit multistep methods for numerical integration of ODE. Implicit implies that the new value can only be found by an iterative process, like Newton’s method, which is computationally costly. Multistep implies that the previous values are needed to compute the new one, which may be costly in memory. BDF is particularly suited for stiff smooth equations, yet the scheme is only A-stable up to order 2. Note that the BDF formula of degree 1 yields the famous backward Euler method.

BDF2

The second order Lagrange polynomials can be written as the following \[ \mathbf{q}^h(T) = \frac{\left(T-T_n\right)\left(T-T_{n-1}\right)}{\left(T_{n+1}-T_n\right)\left(T_{n+1}-T_{n-1}\right)}\mathbf{q}_{n+1} +\frac{\left(T-T_{n+1}\right)\left(T-T_{n-1}\right)}{\left(T_n-T_{n+1}\right)\left(T_n-T_{n-1}\right)}\mathbf{q}_{n} +\frac{\left(T-T_{n+1}\right)\left(T-T_{n}\right)}{\left(T_{n-1}-T_{n+1}\right)\left(T_{n-1}-T_{n}\right)}\mathbf{q}_{n-1} \]

The derivative yields \[ \mathbf{\dot{q}}^h(T) = \frac{2T-T_{n}-T_{n-1}}{(T_{n+1}-T_{n})(T_{n+1}-T_{n-1})}\mathbf{q}_{n+1} + \frac{2T-T_{n+1}-T_{n-1}}{(T_{n}-T_{n+1})(T_{n}-T_{n-1})}\mathbf{q}_{n} + \frac{2T-T_{n+1}-T_{n}}{(T_{n-1}-T_{n+1})(T_{n-1}-T_{n})}\mathbf{q}_{n-1} \] For the generic case \(T=T_{n+1}=t+\tau_1,\;T_n=t,\;T_{n-1}=t-\tau_2\), we have \[ \mathbf{\dot{q}}^h(t+\tau_1) = \frac{2\tau_1+\tau_2}{\tau_1(\tau_1+\tau_2)}\mathbf{q}_{n+1} - \frac{\tau_1+\tau_2}{\tau_1\tau_2}\mathbf{q}_{n} + \frac{\tau_1}{\tau_2(\tau_1+\tau_2)}\mathbf{q}_{n-1} \]

For another representation where \(T=T_{n+1}=t+\tau,\;T_{n}=t+\gamma\tau,\;T_{n-1}=t\), we have \[ \mathbf{\dot{q}}^h(t+\tau) = \frac{1}{\tau}\left(\frac{2-\gamma}{1-\gamma}\mathbf{q}_{n+1} - \frac{1}{\gamma(1-\gamma)}\mathbf{q}_{n} + \frac{1-\gamma}{\gamma}\mathbf{q}_{n-1}\right) \]

For the case where the timestep is fixed \(\tau_1=\tau_2\) or \(\gamma=\tfrac{1}{2}\) , then \[ \mathbf{\dot{q}}^h(t+\tau) = \frac{1}{\tau}\left(\frac{3}{2}\mathbf{q}_{n+1} - 2\mathbf{q}_{n} + \frac{1}{2}\mathbf{q}_{n-1}\right) \]

Derivation using Taylor’s theorem

A simpler way to approach it, at least for fixed timestep, is to use a Taylor expansion of \(\mathbf q\) \[ \mathbf{q}(t+\tau)\approx \mathbf{q}(t) + \tau\mathbf{\dot q}(t)+\frac{\tau^2}{2}\mathbf{\ddot q}(t) \] Replacing the derivatives in the formula by their discrete counterpart yields \[ \mathbf{q}(t+\tau) \approx \mathbf{q}(t) + \tau\frac{\mathbf{q}(t+\tau)-\mathbf{q}(t)}{\tau} +\frac{\tau^2}{2}\frac{\mathbf{q}(t+\tau)-2\mathbf{q}(t)+\mathbf{q}(t-\tau)}{\tau^2} \] where, by simplifying and reorganizing, \[ \mathbf{q}(t+\tau) - \mathbf{q}(t)\approx \frac{3}{2}\mathbf{q}(t+\tau)-2\mathbf{q}(t) +\frac{1}{2}\mathbf{q}(t-\tau) \]