, 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:

Command is

$ stabregion3 -r400 -df BDF234
BDF234, p=2, k=2, l=3
                 1.0000        -2.0000         3.0000
                -4.0000         9.0000       -16.0000
                 3.0000       -18.0000        36.0000
                 0.0000        11.0000       -48.0000
                 0.0000         0.0000        25.0000
                 0.0000         0.0000         0.0000
                 0.0000         0.0000         0.0000
                 2.0000         0.0000         0.0000
                 0.0000         6.0000         0.0000
                 0.0000         0.0000        12.0000
rho_0              0.000000000           0.000000000           0.000000000
rho_1               0.000000000            0.000000000            0.000000000
rho_2               0.000000000            0.000000000            0.000000000
rho_3              -0.222222222            0.000000000            0.000000000   <-----
k=2, l=3, rest=1, n=3, nrest=3, nsq=9, colLen=10
BDF234
A1
              3.00000      0.00000      0.00000
            -18.00000     11.00000      0.00000
             36.00000    -48.00000     25.00000
A0
              0.00000      1.00000     -4.00000
              0.00000     -2.00000      9.00000
              0.00000      3.00000    -16.00000
B1
              2.00000      0.00000      0.00000
              0.00000      6.00000      0.00000
              0.00000      0.00000     12.00000
B0
              0.00000      0.00000      0.00000
              0.00000      0.00000      0.00000
              0.00000      0.00000      0.00000
Parasitic roots of BDF234
        nr      real                    imag                    abs                     3-th root
            0      1.00000000         0.00000000              1.00000000                1.00000000
            1     -0.02545455         0.00000000              0.02545455                0.29416327
            2      0.00000000         0.00000000              0.00000000                0.00000000
            0         1.87500000              3.75624480,         0.00000        90.00000
            0         0.00000000              0.00000000,         0.00000        90.00000
            0         1.87500000             -3.75624480,         0.00000        90.00000
            1         1.89967097              3.76664079,         0.00000        90.00000
. . .
          398        -0.00000000             -0.01047212,        -0.00108        89.90379
          398         1.82627526              3.73500348,        -0.00108        89.90379
          399         1.89967097             -3.76664079,        -0.00108        89.90379
          399        -0.00000000             -0.00523601,        -0.00108        89.90379
          399         1.85053464              3.74569794,        -0.00108        89.90379

3. Instability for variable stepsize. 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)