, 94 min read

Stability Mountains for Fredebeul's Formulas

Beware, this three-dimensional graphic needs some time to load. You can rotate the graphic around any axis.

Here we analyze the three methods from Christoph Fredebeul, which were mentioned in Albrecht (1989). These are three cyclic linear multistep methods of order $p=3$ and $p=4$. They should not to be confused with the A-BDF from Fredebeul.

Bibliography.

  1. Peter Albrecht: “Elements of a General Theory of Composite Integration Methods”, Applied Mathematics and Computation, Vol. 31, May 1989, pp. 1–17
  2. Christoph Fredebeul: “A-BDF: A Generalization of the Backward Differentiation Formulae”, SIAM Journal on Numerical Analysis, Vol. 35, Issue 5, pp.1917–1938

Order 3. The right one is insofar special as the three implicit formulas can be computed in parallel. This is due to the diagonal structure of the $A_1$ and $B_1$ matrix, see below. That makes this cycle particular attractive for GPUs. However, the predictor will likely give worse estimates for the first iterate in Newton iteration.

$$ \begin{array}{r|rrr} p=3 & 1 & 2 & 3 \\ \hline -2 & 12 & 12 & 12 \\ -1 & 0 & 0 & 0 \\ 0 & -72 & -72 & -72 \\ 1 & 60 & 0 & 0 \\ 2 & 0 & 60 & 0 \\ 3 & 0 & 0 & 60 \\ \hline -2 & -9 & -14 & -19 \\ -1 & -6 & 4 & 14 \\ 0 & 21 & 16 & 11 \\ 1 & 30 & 60 & 60 \\ 2 & 0 & 30 & 60 \\ 3 & 0 & 0 & 30 \\ \end{array} \qquad\qquad \begin{array}{r|rrr} p=3 & 1 & 2 & 3 \\ \hline -2 & 30 & 30 & 30 \\ -1 & 0 & 0 & 0 \\ 0 & -90 & -90 & -90 \\ 1 & 60 & 0 & 0 \\ 2 & 0 & 60 & 0 \\ 3 & 0 & 0 & 60 \\ \hline -2 & -3 & 40 & 71 \\ -1 & -66 & -200 & -310 \\ 0 & 51 & 190 & 305 \\ 1 & 18 & 0 & 0 \\ 2 & 0 & 30 & 0 \\ 3 & 0 & 0 & 54 \\ \end{array} $$

Order 4.

$$ \begin{array}{r|rrr} p=4 & 1 & 2 & 3 \\ \hline -2 & 300 & 300 & 300 \\ -1 & 0 & 0 & 0 \\ 0 & -900 & -900 & -900 \\ 1 & 600 & 0 & 0 \\ 2 & 0 & 600 & 0 \\ 3 & 0 & 0 & 600 \\ \hline -2 & -75 & 0 & 113 \\ -1 & -525 & -800 & -1132 \\ 0 & 375 & 700 & 923 \\ 1 & 225 & 400 & 543 \\ 2 & 0 & 300 & 408 \\ 3 & 0 & 0 & 345 \\ \end{array} $$

The method, written as

$$ A_1 y_{n+1} + A_0 y_n = h\cdot \left(B_1 f(y_{n+1} + B_0 f(y_n)\right) $$

has the stability polynomial

$$ Q(\mu,H) = \det\left( A_1\mu + A_0 - H\cdot(B_1\mu + B_0) \right). $$

For $\mathop{\rm Re} H\to -\infty$ the matrix $B_0$ will dominate.

The methods are all $A[\alpha]$-stable but not $A_\infty^0[\alpha]$-stable. See Tendler-like Formulas for Stiff ODEs for the definitions of all these stability properties.

Widlund-wedge angle and distance are as below. For comparison the Widlund-wedge and distance from the corresponding Tendler formulas.

Method α δ r(0) r(∞) Method α δ r(0) r(∞)
Fred3 88.165° 0.0282 0 0.85539679 Tendler3 89.427° 0.0048 0.55371901 0
Fred3p 61.526° 1.0474 0 0.94061062 eTendler3 89.724° 0.0016 0.70756795 0
Fred4 77.345° 0.3766 0 0.88440508 Tendler4 80.882° 0.2441 0.35406989 0

The very high absolute radius at infinity makes the Fredebeul formulas less suitable for stiff system. However, the Fredebeul formulas have two highly sought properties:

  1. Zero parasitic roots
  2. After $k=p$ steps the formulas are stable in step-changing!

The comparison with the Tischer formulas:

  1. All Tischer formulas up to $p\le4$ are A-stable with $r_\infty=0$.
  2. All parasitic roots are zero.

The error constant is

$$ c_{p+1} = \frac{1}{\alpha_{i}\,(p+1)!} \sum_{i=0}^k\bigl(\alpha_ii^{p+1}-(p+1)\beta_ii^p\bigr). $$

The order 3 method.

Fredebeul3, p=3, k=3, l=3
                12.0000        12.0000        12.0000
                 0.0000         0.0000         0.0000
               -72.0000       -72.0000       -72.0000
                60.0000         0.0000         0.0000
                 0.0000        60.0000         0.0000
                 0.0000         0.0000        60.0000
                -9.0000       -14.0000       -19.0000
                -6.0000         4.0000        14.0000
                21.0000        16.0000        11.0000
                30.0000        60.0000        60.0000
                 0.0000        30.0000        60.0000
                 0.0000         0.0000        30.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.000000000            0.000000000            0.000000000
rho_4              -0.125000000           -0.333333333           -0.625000000   <-----

The parallel order 3 method.

Fredebeul3p, p=3, k=3, l=3
                30.0000        30.0000        30.0000
                 0.0000         0.0000         0.0000
               -90.0000       -90.0000       -90.0000
                60.0000         0.0000         0.0000
                 0.0000        60.0000         0.0000
                 0.0000         0.0000        60.0000
                -3.0000        40.0000        71.0000
               -66.0000      -200.0000      -310.0000
                51.0000       190.0000       305.0000
                18.0000         0.0000         0.0000
                 0.0000        30.0000         0.0000
                 0.0000         0.0000        54.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.000000000            0.000000000            0.000000000
rho_4               0.075000000            0.666666667            0.375000000   <-----

The order 4 method.

Fredebeul4, p=4, k=3, l=3
               300.0000       300.0000       300.0000
                 0.0000         0.0000         0.0000
              -900.0000      -900.0000      -900.0000
               600.0000         0.0000         0.0000
                 0.0000       600.0000         0.0000
                 0.0000         0.0000       600.0000
               -75.0000         0.0000       113.0000
              -525.0000      -800.0000     -1132.0000
               375.0000       700.0000       923.0000
               225.0000       400.0000       543.0000
                 0.0000       300.0000       408.0000
                 0.0000         0.0000       345.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.000000000            0.000000000            0.000000000
rho_4               0.000000000            0.000000000            0.000000000
rho_5              -0.020833333           -0.172222222           -0.586944444   <-----

1. Stability regions

Using stabregion3 to graph the stability regions.

stabregion3 -f Fredebeul3 -oj -r1800

stabregion3 -f Fredebeul3p -oj -r1800

stabregion3 -f Fredebeul4 -oj -r1800

2. Stability mountain

Below is the output of:

stabregion3 -f Fredebeul3 -o3

Fredebeul3 stability mountain. So visually it is obvious that the method is not $A_\infty^0[\alpha]$-stable.

Output for the "parallel" order 3 method.

stabregion3 -f Fredebeul3p -o3 -r100 -L4:-15:10:25

Again, this method is not $A_\infty^0[\alpha]$-stable.

Output for the order 4 method.

stabregion3 -f Fredebeul4 -o3 -r100 -L5:-20:20:20

This method is not $A_\infty^0[\alpha]$-stable.