1st February 2026, 103 min read

Rubin's 4-th Order Method is Neither A-stable Nor D-stable

Original post is here eklausmeier.goip.de/blog/2026/02-01-rubin-4-th-order-method-is-neither-a-stable-nor-d-stable.

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

This is in continuation of:

We analyze below method from Rubin, see his Fig. 4.2. Also see Weitere zyklische und blockimplizite lineare Mehrschritt-Verfahren.

$$ \begin{array}{r|rr} p=4 && 1 & 2\cr%R4A \hline -1 && 0 & 56\cr 0 && 24 & -72\cr 1 && -24 & 0\cr 2 && 0 & 16\cr \hline -1 && 1 & -21\cr 0 && -13 & -39\cr 1 && -13 & 33\cr 2 && 1 & 3\cr \hline c_{5i} && 0.0153 & 0.08125\cr \end{array} $$

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 stability polynomial has roots at 1 and 3.5, and therefore by definition is not D-stable. By a continuity argument it can hence not be A-stable. I.e., in the vicinity of zero its roots are less than one in magnitude, by definition of A-stability. But at exactly zero the zero jumps to 3.5. That cannot be.

Rubin1, p=4, k=2, l=2
            0.0000        56.0000
           24.0000       -72.0000
          -24.0000         0.0000
            0.0000        16.0000
            1.0000       -21.0000
          -13.0000       -39.0000
          -13.0000        33.0000
            1.0000         3.0000
rho_0       0.000000000           0.000000000
rho_1      -0.000000000           0.000000000
rho_2      -0.000000000           0.000000000
rho_3      -0.000000000           0.000000000
rho_4      -0.000000000           0.000000000
rho_5       0.015277778           0.081250000  <-----

1. Stability region.

Below is the output of:

stabregion2 -f Rubin1 -oj -r600

Just looking at this stability region one could assume that the method is A-stable.

Truth is, it is not.

2. Stability mountain.

Below is the output of:

stabregion2 -f Rubin1 -o3 -r600 -L29

Rubin1 stability mountain.