, 2 min read

Cyclic Combination of BDF

1. Backward differentiation formulas. The BDF are as follows.

Name Formula
BDF1 $y_{n+1} - y_n = h f_{n+1}$
BDF2 $3 y_{n+2} - 4 y_{n+1} + y_n = 2 h f_{n+2}$
BDF3 $11 y_{n+3} - 18 y_{n+2} + 9 y_{n+1} - 2 y_n = 6 h f_{n+3}$
BDF4 $25 y_{n+4} - 48 y_{n+3} + 36 y_{n+2} - 16 y_{n+1} + 3 y_n = 12 h f_{n+4}$
BDF5 $137 y_{n+5} - 300 y_{n+4} + 300 y_{n+3} - 200 y_{n+2} + 75 y_{n+1} -12 y_n = 60 h f_{n+5}$
BDF6 $147 y_{n+6} - 360 y_{n+5} + 450 y_{n+4} - 400 y_{n+3} + 225 y_{n+2} - 72 y_{n+1} + 10 y_n = 60 h f_{n+6}$

2. Cycle. Combining BDF2, BDF3, BDF4 in a cyclic manner leads to the two matrix polynomials

$$ \rho(\mu) = A_1\mu + A_0, \quad \sigma(\mu) = B_1\mu + B_0 $$

with coefficients

$$ A_0 = \begin{pmatrix} 0 & 1 & -4\\ 0 & -2 & 9\\ 0 & 3 & -16 \end{pmatrix}, \, A_1 = \begin{pmatrix} 3 & 0 & 0\\ -18 & 11 & 0\\ 36 & -48 & 25 \end{pmatrix}, \quad B_0 = 0, \, B_1 = \begin{pmatrix} 2 & 0 & 0\\ 0 & 6 & 0\\ 0 & 0 & 12 \end{pmatrix}. $$

Using stabregion3.c one computes:

3. Instability. Nishikawa (2019) showed that a certain combination of BDF1 and BDF2 can lead to instabilities. However, the combination of BDF1 and BDF2 for fixed stepsizes will remain L-stable. The same holds true for any combination:

4. Cycle as product of multistep methods. Instead of using the matrices $A_1, A_0, B_1, B_0$ one can write the above 3-stage cycle as

$$ Y_{n+1} = K Y_n + h F(Y_n, Y_{n+1}), $$

where

$$ K = K_4 \, K_3 \, K_2 $$

with

$$ K_2 = \begin{pmatrix} 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1\\ 0 & 0 & -\frac{1}{3} & \frac{4}{3} \end{pmatrix}, \, K_3 = \begin{pmatrix} 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1\\ 0 & -\frac{2}{11} & -\frac{9}{11} & \frac{11}{11} \end{pmatrix}, \, K_4 = \begin{pmatrix} 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1\\ -\frac{3}{25} & \frac{16}{25} & -\frac{36}{25} & \frac{48}{25} \end{pmatrix} $$

The eigenvalues of $\det(A_1\mu+A_0)$ and $\det(I\mu-K)$ are the same. They are

$$ \begin{pmatrix} 0\\ 0\\ -0.02545455\\ 1 \end{pmatrix} $$

For example, you can use Octave:

k2=[0, 1, 0, 0; 0, 0, 1, 0; 0, 0, 0, 1; 0, 0, -1/3, 4/3]
k3=[0, 1, 0, 0; 0, 0, 1, 0; 0, 0, 0, 1; 0, 2/11, -9/11, 18/11]
k4=[0, 1, 0, 0; 0, 0, 1, 0; 0, 0, 0, 1; -3/25, 16/25, -36/25, 48/25]
eig(k4 * k3 * k2)

Similarly

a0=[0,1,-4; 0,-2,9; 0,3,-16]
a1=[3,0,0; -18,11,0; 36,-48,25]
eig(a0,-a1)