, 71 min read

Comparing BDF vs. Tendler vs. Tischer formulas

Original post is here eklausmeier.goip.de/blog/2025/12-17-comparing-bdf-vs-tendler-vs-tischer.


1. Dahlquist's test equation

We tested the BDF, Tendler's formulas, new Tendler-like formulas, and Tischer's formulas on the classical test equation

$$ y´(t) = \lambda y \in \mathbb{C}, \quad y(t_0) = e^{\lambda t_0}, \quad t\in[t_0, t_1]. $$

We used

$$ \lambda = r e^{i \varphi}, \quad t_0=0, \quad t1=-40, \quad h\lt 0. $$

The radius is fixed to $r=100$, and $\varphi$ varies from 5° to 90° in steps of 5°. $\varphi$ directly tests the Widlund wedge angle of the formula.

The stepsize $h$ of each formula is varied repeatedly from -0.1, -0.01, and -0.001.

Above differential equation is now solved with the following formulas:

  1. BDF order 1–6
  2. Tendler order 3–7
  3. new Tendler-like formulas order 3–9, see Tendler-like formulas
  4. Tischer's formulas of order 2–8 using $S=0$

This creates 1350 data points. Multiply by 2 machines, and two precisions, gives 5400 records.

We report the summed global error computed as

$$ g_\hbox{err} = \sum_{i=0}^n |y(t_i) - y_i|, \quad n={t_1-t_0\over h}. $$

Computations were done in double precision (double complex). Results are here. As the global error $|y(t_i) - y_i|$ per step varies wildly between $10^{-324}$ to $10^{+304}$, we computed

$$ \hat g_\hbox{err} = \begin{cases} 5, & \hbox{if } g_\hbox{err} \gt 30 \hbox{ or NaN}\\ \log_{10} g_\hbox{err}, & \hbox{else} \end{cases} $$

That's what is shown in below graphic.

Computations were done on below machine.

Hard- and software Machine 1 Machine 2
CPU AMD Ryzen 7 5700G with Radeon Graphics ARM Cortex-A77
CPU max MHz 4673.8232 3187.2000
C compiler gcc version 15.2.1 20251112 (GCC) clang version 21.1.7
Floating point arithmetic IEEE-754 standard compliant floating-point operator (with only minor documented deviations) IEEE 754 standard

2. Logarithmic error in double precision

Below 3D chart can be rotated and zoomed in or out. Clicking on any entry shows the log global error $\hat g_\hbox{err}$. Double precision on AMD Ryzen.

Once can clearly see that the higher order methods quickly lose precision when a higher Widlund wedge angle is required. This is the reason why GEAR, EPISODE, LSODE and CVODE all do not use BDF6. However, in reality, this is mainly a problem for the stepsize and order control segment to properly switch order.

For the new Tendler-like formulas we intend to use even the higher order methods of order 7, 8, and 9 as we aim to provide a type-insensitive code, which just switches between fixed point iteration and modified Newton method. It is well known, see Montenbruck/Gill (2000) §4.1.6, that higher order methods are indeed required for certain precisions.

3. Logarithmic error in single precision and other machines

As Gaffney (1984) noted a sensitivity regarding the used precision, we repeated above test in single precision (float complex) and on another machine. The qualitative results didn't differ in any way. Results are here.

1. Single precision on AMD Ryzen. Results are here.

2. Double precision on ARM Cortex. Results are here.

3. Single precision on ARM Cortex. Results are here.

4. Scripts to generate results

We used the C program stabregion2.c to generate above results. Below shell script stabregion2-batch was once run with F=t (lower case) and once with F=T (capital), while varying $\varphi$ and $h$.

#!/bin/bash
# Run stabregion2 with various parameters for multiple formulas

F=T

echo BDF1-6
for i in 1 2 3 4 5 6; do stabregion2 -f BDF$i -${F}100:$1:0:-40:$2 | tail -1; done
echo Tendler3-7
for i in 3 4 5 6 7; do stabregion2 -f Tendler$i -${F}100:$1:0:-40:$2 | tail -1; done
echo eTendler3-9
for i in 3 4 5 6 7 8 9; do stabregion2 -f eTendler$i -${F}100:$1:0:-40:$2 | tail -1; done
echo Tischer2-8
for i in 2 3 4 5 6 7 8 ; do stabregion2 -f Tischer$i -${F}100:$1:0:-40:$2 | tail -1; done

Above script was called like so:

for phi in `seq 5 5 90`; do echo phi=$phi; for h in -0.1 -0.01 -0.001; do echo h=$h; stabregion2-batch $phi $h; done; done

The output of above script was then processed with below Perl script to generate JSON for Apache ECharts:

#!/bin/perl -W
# Generate JSON + ECharts commands from stabregion2-test-phi1-phi2-h1-h2.txt
# That file was generated via
#    for phi in `seq 5 5 90`; do echo phi=$phi; for h in -0.1 -0.01 -0.001; do echo h=$h; stabregion2-batch $phi $h; done; done

use strict;
my ($log10,$phi,$h,$formula,$p,$gerr,@F) = (log(10),0,0,"",0,0,());

print "[\n\t[\"formula\", \"h_phi\", \"gerr\"],\n";
while (<>) {
    if (/^phi=(\d+)/) { $phi = $1; next; }
    if (/^h=([+-]\d+\.\d+)/) { $h = $1; next; }
    if (/^(BDF|Tendler|eTendler|Tischer)(\d)-(\d)/) { $formula = $1; $p = $2 - 1; next; }
    if (/^\s*\d+\s+[+-]\d+\.\d+/) {
        @F = split(/\s+/);
        $gerr = $F[6];
        $p += 1;
        if ($gerr =~ /(-nan|nan)/ || $gerr > 30) { $gerr = 5; }
        else { $gerr = log($gerr) / $log10; }
    }
    #printf("phi=%d, h=%g, formula=%s, p=%d, gerr=%g\n",$phi,$h,$formula,$p,$gerr);
    printf("\t[ \"%s%d\", \"%g/%g\", %g],\n",$formula,$p,$h,$phi,$gerr);
}
print "];\n";