, 9 min read

Stability Mountain for Majid/Suleiman/Omar

Original post is here eklausmeier.goip.de/blog/2026/02-18-stability-mountain-for-majid-suleiman-omar.


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

Here we analyze the method from Majid/Suleiman/Omar. This is a block implicit method of order 4.

$$ \begin{array}{r|rrr} p=4 & 1 & 2 & 3\cr \hline -1 & -24 & 0 & 0 \cr 0 & 24 & -24 & 0 \cr 1 & 0 & 24 & -24 \cr 2 & 0 & 0 & 24 \cr \hline -1 & 9 & -1 & 1 \cr 0 & 19 & 13 & -5 \cr 1 & -5 & 13 & 19 \cr 2 & 1 & -1 & 9 \cr \hline c_{5i} & -0.026 & 0.015 & -0.026\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). $$ 
Majid, p=4, k=1, l=3
       -24.0000         0.0000         0.0000
        24.0000       -24.0000         0.0000
         0.0000        24.0000       -24.0000
         0.0000         0.0000        24.0000
         9.0000        -1.0000         1.0000
        19.0000        13.0000        -5.0000
        -5.0000        13.0000        19.0000
         1.0000        -1.0000         9.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.026388889           0.015277778          -0.026388889   <-----

Parasitic roots of Majid
        nr      real                    imag                    abs                     3-th root
         0      1.00000000         0.00000000              1.00000000                1.00000000
         1      0.00000000         0.00000000              0.00000000                0.00000000
         2      0.00000000         0.00000000              0.00000000                0.00000000

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). $$

The matrices are

$$ A_1 = \pmatrix{ 24 & 0 & 0\cr -24 & 24 & 0\cr 0 & -24 & 24 }, \quad A_0 = \pmatrix{ 0 & 0 & -24\cr 0&0&0\cr 0&0&0 }, \quad B_1 = \pmatrix{ 19 & -5 & 1\cr 13 & 13 & -1\cr -5 & 19 & 9 }, \quad B_0 = \pmatrix{ 0&0&9\cr 0&0&-1\cr 0&0&1 }. $$

For $\mathop{\rm Re} H\to -\infty$ the matrix $B_0$ will dominate. Matrix $B_0$ has eigenvalues 1, 0, 0.

We assume that the implicitness is resolved by a full Newton method. I.e., for a scaler differential equation we would have to solve a three-dimensional equation at each step. The method is $A$-stable but not $A_\infty^0$-stable (also called $L$-stable). See Tendler-like Formulas for Stiff ODEs for the definitions of all these stability properties.

Majid/Suleiman/Omar, however, treat the implicitness via fixed point iteration. Also, they do not solve a three-dimensional system but solve the block by either Jacobi or Gauß-Seidel iteration. In that case the method is no longer $A$-stable. The stability diagrams are given in their paper.

It would be interesting to investigate the effect of Gauß-Seidel combined with a Newton iteration per stage.

1. Stability region.

Below is the output of:

stabregion3 -f MajidOmar -oj -r120

2. Stability mountain.

Below is the output of:

stabregion3 -fMajidOmar -o3 -L5:-20:5:20

Majid/Suleiman/Omar stability mountain. One can clearly see a trough around zero. Also, clearly the roots approach 1 at infinity. So visually it is obvious that the method is not $A_\infty^0$-stable.