, 2 min read

Computing Error Constants in High Precision

In Searching for Tendler-like formulas we presented new Tendler-like cyclic linear multistep formulas for the solution of stiff differential equations. Assume the cyclic linear multistep methods of the form

$$ \sum_{j=-k+1}^\ell \left[ \alpha_{ij} y_{m\ell+j} - h \beta_{ij} \dot y_{m\ell+j} \right] = 0, \quad i=1,\ldots,\ell. $$

When testing these formulas numerically, we are particularly interested in the unscaled error constant

$$ c_{p+1} = \frac{1}{(p+1)!} \sum_{i=0}^k\bigl(\alpha_ii^{p+1}-(p+1)\beta_ii^p\bigr). $$

The scaled error constant is

$$ \eta_{i,p+1} = \frac{1}{(p+1)!} \frac{1}{\alpha_{ii}} c_{i,p+1} $$

Unfortunately, for the new formulas eTendler8 and eTendler9 of order 8 and 9 the computation in double precision (double in the C programming language) gives downright wrong answers.

However, the standard Unix tool bc allows computation with arbitrary precision and is usually installed anyways.

Order 8.

$$ \begin{array}{r|rrrr} p=8 & 1 & 2 & 3 & 4\\ \hline -7 & 105 & 0 & 0 & 0\\ -6 & -960 & 10560 & 0 & 0\\ -5 & 3920 & -96740 & 4350 & 0\\ -4 & -9408 & 396116 & -40060 & 11580\\ -3 & 14700 & -954618 & 165256 & -106094\\ -2 & -15680 & 1501850 & -402822 & 434406\\ -1 & 11760 & -1623860 & 646450 & -1046346\\ 0 & -6720 & 1267140 & -731500 & 1640450\\ 1 & 2283 & -701166 & 591360 & -1801730\\ 2 & 0 & 200718 & -290706 & 1438794\\ 3 & 0 & 0 & 57672 & -782406\\ 4 & 0 & 0 & 0 & 211346\\ \hline -7 & 0 & 0 & 0 & 0\\ -6 & 0 & 0 & 0 & 0\\ -5 & 0 & 0 & 0 & 0\\ -4 & 0 & 0 & 0 & 0\\ -3 & 0 & 0 & 0 & 0\\ -2 & 0 & 0 & 0 & 0\\ -1 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0\\ 1 & 840 & -56280 & 25200 & 21000\\ 2 & 0 & 76440 & -64680 & 2520\\ 3 & 0 & 0 & 24360 & -81480\\ 4 & 0 & 0 & 0 & 81480\\ \end{array} $$

To compute the unscaled and the scaled error constant we used below bc script:

a1[0]=105
a1[1]=-960
a1[2]=3920
a1[3]=-9408
a1[4]=14700
a1[5]=-15680
a1[6]=11760
a1[7]=-6720
a1[8]=2283
b1[0]=0
b1[1]=0
b1[2]=0
b1[3]=0
b1[4]=0
b1[5]=0
b1[6]=0
b1[7]=0
b1[8]=840

a4[0]=11580
a4[1]=-106094
a4[2]=434406
a4[3]=-1046346
a4[4]=1640450
a4[5]=-1801730
a4[6]=1438794
a4[7]=-782406
a4[8]=211346
b4[0]=0
b4[1]=0
b4[2]=0
b4[3]=0
b4[4]=0
b4[5]=21000
b4[6]=2520
b4[7]=-81480
b4[8]=81480

With these constants in place we then ran below for-loop for the various columns:

scale=10
f=1
for (q=1; q<=9; ++q) {
    f *= q + 1;
    r=0
    for (k=0; k<=8; ++k)
        r += k^q * a4[k] - q * k^(q-1) * b4[k]
    r
    r / f / a4[8]
}

Order 9.

$$ \begin{array}{r|rrrrr} p=9 & 1 & 2 & 3 & 4 & 5\\ \hline -8 & -280 & 0 & 0 & 0 & 0\\ -7 & 2835 & -5285 & 0 & 0 & 0\\ -6 & -12960 & 53730 & -13715 & 0 & 0\\ -5 & 35280 & -246960 & 138885 & -24780 & 0\\ -4 & -63504 & 677376 & -634992 & 250764 & -22331\\ -3 & 79380 & -1233036 & 1728720 & -1145544 & 225768\\ -2 & -70560 & 1569960 & -3111108 & 3115434 & -1029642\\ -1 & 45360 & -1446480 & 3883740 & -5600364 & 2789808\\ 0 & -22680 & 1028160 & -3422160 & 6991530 & -4946214\\ 1 & 7129 & -486351 & 2295792 & -6110664 & 6531756\\ 2 & 0 & 88886 & -1194345 & 3889494 & -5933718\\ 3 & 0 & 0 & 329183 & -2019384 & 3364992\\ 4 & 0 & 0 & 0 & 653514 & -1609983\\ 5 & 0 & 0 & 0 & 0 & 629564\\ \hline -8 & 0 & 0 & 0 & 0 & 0\\ -7 & 0 & 0 & 0 & 0 & 0\\ -6 & 0 & 0 & 0 & 0 & 0\\ -5 & 0 & 0 & 0 & 0 & 0\\ -4 & 0 & 0 & 0 & 0 & 0\\ -3 & 0 & 0 & 0 & 0 & 0\\ -2 & 0 & 0 & 0 & 0 & 0\\ -1 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0\\ 1 & 2520 & -98280 & -80640 & -40320 & -241920\\ 2 & 0 & 35280 & -63000 & -73080 & -168840\\ 3 & 0 & 0 & 118440 & 35280 & 171360\\ 4 & 0 & 0 & 0 & 229320 & 171360\\ 5 & 0 & 0 & 0 & 0 & 216720\\ \end{array} $$

Setting the constants in bc:

a1[0]=-280
a1[1]=2835
a1[2]=-12960
a1[3]=35280
a1[4]=-63504
a1[5]=79380
a1[6]=-70560
a1[7]=45360
a1[8]=-22680
a1[9]=7129
b1[0]=0
b1[1]=0
b1[2]=0
b1[3]=0
b1[4]=0
b1[5]=0
b1[6]=0
b1[7]=0
b1[8]=0
b1[9]=2520

a2[0]=-5285
a2[1]=53730
a2[2]=-246960
a2[3]=677376
a2[4]=-1233036
a2[5]=1569960
a2[6]=-1446480
a2[7]=1028160
a2[8]=-486351
a2[9]=88886
b2[0]=0
b2[1]=0
b2[2]=0
b2[3]=0
b2[4]=0
b2[5]=0
b2[6]=0
b2[7]=0
b2[8]=-98280
b2[9]=35280

a3[0]=-13715
a3[1]=138885
a3[2]=-634992
a3[3]=1728720
a3[4]=-3111108
a3[5]=3883740
a3[6]=-3422160
a3[7]=2295792
a3[8]=-1194345
a3[9]=329183
b3[0]=0
b3[1]=0
b3[2]=0
b3[3]=0
b3[4]=0
b3[5]=0
b3[6]=0
b3[7]=-80640
b3[8]=-63000
b3[9]=118440

a4[0]=-24780
a4[1]=250764
a4[2]=-1145544
a4[3]=3115434
a4[4]=-5600364
a4[5]=6991530
a4[6]=-6110664
a4[7]=3889494
a4[8]=-2019384
a4[9]=653514
b4[0]=0
b4[1]=0
b4[2]=0
b4[3]=0
b4[4]=0
b4[5]=0
b4[6]=-40320
b4[7]=-73080
b4[8]=35280
b4[9]=229320

a5[0]=-22331
a5[1]=225768
a5[2]=-1029642
a5[3]=2789808
a5[4]=-4946214
a5[5]=6531756
a5[6]=-5933718
a5[7]=3364992
a5[8]=-1609983
a5[9]=629564
b5[0]=0
b5[1]=0
b5[2]=0
b5[3]=0
b5[4]=0
b5[5]=-241920
b5[6]=-168840
b5[7]=171360
b5[8]=171360
b5[9]=216720

Now run below for-loop for the five columns.

scale=30
f=1
for (q=1; q<=10; ++q) {
    f *= q;
    r=0
    for (k=0; k<=9; ++k)
        r += k^q * a5[k] - q * k^(q-1) * b5[k]
    r
    r / f / a5[9]
}

These are the scaled error constants given in Tendler-like Formulas for Stiff ODEs.