Backward Differentiation Formulas (BDF) are linear recurrence relations that approximate the derivative of a function at a time $$t_n$$ in terms of its function values $$y(t_n)$$ and k earlier values. The $$k^{th}$$ order BDF is derived by performing a Taylor series expansion of the function at k previous time steps:

$$y(t_n-h) = y(t_n) - h\left.\frac{dy}{dt}\right|_{t=t_n} + \frac{h^2}{2}\left.\frac{d^2y}{dt^2}\right|_{t=t_n} - \frac{h^3}{6}\left.\frac{d^3y}{dt^3}\right|_{t=t_n} + \frac{h^4}{24}\left.\frac{d^4y}{dt^4}\right|_{t=t_n} - \dots \\ \label{eq:bdf1}$$ $$y(t_n-2h) = y(t_n) - 2h\left.\frac{dy}{dt}\right|_{t=t_n} + \frac{4h^2}{2}\left.\frac{d^2y}{dt^2}\right|_{t=t_n} - \frac{8h^3}{6}\left.\frac{d^3y}{dt^3}\right|_{t=t_n} + \frac{16h^4}{24}\left.\frac{d^4y}{dt^4}\right|_{t=t_n} - \dots \\ \vdots \\ \label{eq:bdf2}$$ $$y(t_n-kh) = y(t_n) - kh\left.\frac{dy}{dt}\right|_{t=t_n} + \frac{(kh)^2}{2}\left.\frac{d^2y}{dt^2}\right|_{t=t_n} - \frac{(kh)^3}{6}\left.\frac{d^3y}{dt^3}\right|_{t=t_n} + \frac{(kh)^4}{24}\left.\frac{d^4y}{dt^4}\right|_{t=t_n} - \dots \\ \label{eq:bdf4}$$

The generation of a particular BDF recurrence relation involves adding and subtracting multiples of the above formulas such that derivatives of order one and greater are eliminated from the resulting equation. For example, for the $$2^{nd}$$ order BDF, multiplying equation \eqref{eq:bdf1} by 4 and subtracting equation \eqref{eq:bdf2} from the result to eliminate the $$2^{nd}$$ order derivates, and ignoring terms of order greater than $$h^3$$ (since they are negligible for small h):

$$\left.\frac{dy}{dt}\right|_{t=t_n} = \frac{1}{2h}\left[y(t_n-2h) - 4y(t_n-h) + 3y(t_n)\right] + \frac{h^2}{3}\left.\frac{d^3y}{dt^3}\right|_{t=t_n} \\ \label{eq:bdf5}$$ $$\therefore \left.\frac{dy}{dt}\right|_{t=t_n} = \frac{1}{2h}\left[y(t_n-2h) - 4y(t_n-h) + 3y(t_n)\right] + O(h^2) \\ \label{eq:bdf6}$$

In this example, the dominant error term is proportional to $$O(h^2)$$. In general, the dominant error term is proportional to $$O(h^k)$$. Therefore, as the order of the BDF method is increased, the approximation error decreases.

In general, an order k bdf method may be summerised as:

$$\left.\frac{dy}{dt}\right|_{t=t_n} = \frac{1}{h}\sum_{j=0}^{k}a_{kj}y(t_n - jh) \\ \label{eq:bdf7}$$

The constants $$a_{kj}$$, for BDF orders 1 to 6, are given in the Table below

Order k $$a_{k0}$$ $$a_{k1}$$ $$a_{k2}$$ $$a_{k3}$$ $$a_{k4}$$ $$a_{k5}$$ $$a_{k6}$$
1 1 -1
2 $$\frac{3}{2}$$ -2 $$\frac{1}{2}$$
3 $$\frac{11}{6}$$ -3 $$\frac{3}{2}$$ $$-\frac{1}{3}$$
4 $$\frac{25}{12}$$ -4 3 $$-\frac{4}{3}$$ $$\frac{1}{4}$$
5 $$\frac{137}{60}$$ -5 5 $$-\frac{10}{3}$$ $$\frac{5}{4}$$ $$-\frac{1}{5}$$
6 $$\frac{147}{60}$$ -6 $$\frac{15}{2}$$ $$-\frac{20}{3}$$ $$\frac{15}{4}$$ $$-\frac{6}{5}$$ $$-\frac{1}{6}$$