1st June 2025, 68 min read

Stability Regions for BDF and Tendler's Formulas

In Die verwendeten zyklischen Formeln im Programm TENDLER we listed the Widlund wedge α and Widlund distance δ for the seven formulas of Joel Tendler, which he developed during his PhD thesis in 1973.

Here we compute these values by using the QZ algorithm. We list the values for the BDF (backward differentiation formulas) and the cyclic linear multistep formulas by Tendler.

Also see Design Notes on System for the Analysis of Order- and Stepsize Changes for Cyclic Composite Multistep Methods.

1. The BDF.

The two first BDF are A-stable. Therefore the Widlund wedge α must be 90°, as can be seen by δ being zero.

You can hover over the diagrams and see the values.

2. Tendler's cyclic formulas.

3. Program stabregion. The C program to generate these values uses two LAPACK routines:

extern double dlamch_(char *cmach);
extern void zggev_(char *jobvl, char *jobvr, int *n, complex *a, int *lda, complex *b,
    int *ldb, complex *alpha, complex *beta, complex *vl, int *ldvl, complex *vr, int *ldvr,
    complex *work, int *lwork, double *rwork, int *info);

It stores the formulas like so:

typedef struct {
    const char *name;
    const int p;	// order (for information only)
    const int k;	// number of start steps, 2*(k+l) = number of rows of a[]
    const int l;	// cycle length, i.e., column count of a[]
    double *a;	// transpose of A_i and B_i matrices in a single matrix
} formula_t;

formula_t F[] = {
    {
        "BDF1", 1, 1, 1,
        (double[]){
        -1, 1,
        0, 1 }
    },
    {
        "BDF2", 2, 2, 1,
        (double[]){
        1, -4, 3,
        0, 0, 2 }
    },
    . . .
    {
        "Tendler3", 3, 3, 3,	// name, p, k, l
        (double[]){
        -2, 0, 0,
        9, -2, 0,
        -18, 9, 0,
        11, -18, 9,
        0, 11, -12,
        0, 0, 3,
        // --------
        0, 0, 0,
        0, 0, 0,
        0, 0, 0,
        6, 0, -4,
        0, 6, -4,
        0, 0, 2 }
    },
    . . .
};

It "transposes" above formulas in to matrices A1, B1, A0, and B0:

if (fm->k >= fm->l) {
    rest = fm->k - fm->l,
    n = fm->k,
    nrest = fm->l,
    nsq = fm->k * fm->k;
} else {
    rest = fm->l - fm->k,
    n = fm->l,
    nrest = fm->k,
    nsq = fm->l * fm->l;
}

memset(a0,0,sizeof(*a0) * nsq);	// set a0[] to zero matrix
memset(b0,0,sizeof(*b0) * nsq);	// set b0[] to zero matrix
memset(a1,0,sizeof(*a1) * nsq);	// set a1[] to zero matrix
memset(b1,0,sizeof(*b1) * nsq);	// set b1[] to zero matrix

for (i=0; i<rest; ++i) {
    a1[i+i*n] = 1;	// fill with parts of identity matrix
    if (i + nrest < n) a0[i+(i+nrest)*n] = -1;	// cancel out
}
for (i=rest; i<n; ++i) {
    for (j=0; j<n; ++j)
        a0[i+j*n] = fm->a[(i-rest)+j*fm->l];
    for (j=rest; j<n; ++j)
        a1[i+j*n] = fm->a[(i-rest)+((j-rest)+n)*fm->l],
        b1[i+j*n] = fm->a[(i-rest)+((j-rest+fm->l)+2*n)*fm->l];
}

Finally, it iterates through the unit circle and computes all eigenvalues.

for (i=0,phiinc=2*M_PI/nr,phi=0; i<nr; phi+=phiinc,++i) {
    eiphi = cexp(I * phi);
    for (j=0; j<nsq; ++j)	// A = A1 e^{i\phi} + A0, same for B
        a[j] = a1[j] * eiphi + a0[j],
        b[j] = b1[j] * eiphi + b0[j];
    info = 0;
    zggev_(jobvl, jobvr, &n, a, &lda, b, &ldb, alpha, beta, vl, &ldvl, vr, &ldvr, work, &lwork, rwork, &info);
    for (j=0; j<n; ++j) {
        if (beta[j] == 0) continue;
        lambda = alpha[j] / beta[j];
        rlambda = creal(lambda),
        ilambda = cimag(lambda),
        xmin = fmin(xmin,rlambda);	// Widlund distance delta
        if (rlambda < 0) {
            if (fabs(rlambda) > eps && fabs(ilambda) > eps)
                widlw = fmin(widlw,fabs(atan(ilambda/rlambda)));	// Widlund wedge in radians
        }
        prtf(i,lambda,xmin,180/M_PI*widlw);
}

The function pointer prtf() points to a function that prints JavaScript suitable for using in Apache ECharts.

3. Running program. The program stabregion is controlled by command line arguments.

1. With the output flag -oj we generate JavaScript.
2. Command line flags -b are BDF, -t are for Tendler's formulas of the respective order.
3. Command line flag -r specifies how many subdivisions we do for the unit circle.

rm /tmp/formula
for i in `seq 1 6`; do stabregion -oj -b$i -r600 >> /tmp/formula; done
for i in `seq 3 7`; do stabregion -oj -t$i -r600 >> /tmp/formula; done