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.

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.

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
        if (fabs(rlambda) > eps && fabs(ilambda) > eps) widlw = fmin(widlw,carg(lambda));	// Widlund wedge in radians
        prtf(i,lambda,xmin,180+180/M_PI*widlw);
    }
}

The function pointer prtf() 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