, 27 min read
Stability Regions for Varah's Formulas
James Varah created linear multistep method with zero roots at infinity, i.e., $A_\infty^0[\alpha]$-stable methods. These methods share this property with the BDF. For the terminology see Tendler-like Formulas for Stiff ODEs. This post gives the details of the formulas and plots the stability regions.
References.
- Varah, James M.: “Stiffly Stable Linear Multistep Methods of Extended Order”, SIAM Journal on Numerical Analysis, 1978, Vol. 15, Issue 6, pp.123—1246
Varah (1978) gives his multistep methods in terms of the $r(z)$, with
$$
r(z) = \left({1-z\over2}\right)^k \, \rho\left({1+z\over1-z}\right) = \sum_{j=0}^k a_j z^j.
$$
$\rho()$ and $\sigma()$ being the generating polynomials of the linear multistep method.
| Coefficient | $k=7$ $p=5$ |
$k=8$ $p=6$ |
$k=9$ $p=7$ |
$k=10$ $p=7$ |
$k=11$ $p=8$ |
$k=12$ $p=8$ |
|---|---|---|---|---|---|---|
| $a_1$ | 2.0 | 2.0 | 2.0 | 2.0 | 2.0 | 2.0 |
| $a_2$ | 14.0 | 16.0 | 18.0 | 20.0 | 22.0 | 24.0 |
| $a_3$ | 42.667 | 56.667 | 72.667 | 90.667 | 110.67 | 132.67 |
| $a_4$ | 74.667 | 117.33 | 174.0 | 246.67 | 337.33 | 448.0 |
| $a_5$ | 84.4 | 159.07 | 276.4 | 450.4 | 697.06 | 1034.4 |
| $a_6$ | 46.35 | 152.53 | 311.6 | 588.0 | 1038.4 | 1735.5 |
| $a_7$ | 20.51 | 69.65 | 266.69 | 578.28 | 1166.3 | 2204.7 |
| $a_8$ | 31.29 | 113.32 | 325.0 | 1037.1 | 2203.4 | |
| $a_9$ | 50.25 | 162.3 | 531.2 | 1388.6 | ||
| $a_{10}$ | 20.0 | 270.0 | 873.6 | |||
| $a_{11}$ | 32.8 | 187.91 | ||||
| $a_{12}$ | 79.0 |
Mapping the left half-plane back to the unit circle gives below classical coefficients for the generating polynomial $\rho(\xi)$,
$$
\rho(\xi) = \sum_{i=0}^k \alpha_i \xi^i.
$$
I asked ChatGPT to do this for me. The generating polynomial $\sigma(\zeta) = \zeta^k$, i.e., $\beta_k=1$, all other $\beta_i=0$. So, only the $\alpha_i$ are needed.
The inverse transformation from $r(z)$ to $\rho(\xi)$ is
$$
\rho(\xi) = {(\xi+1)^k\over 2^k} \, r\left({\xi-1\over\xi+1}\right).
$$
| Coeff. | $k=7$ $p=5$ |
$k=8$ $p=6$ |
$k=9$ $p=7$ |
$k=10$ $p=7$ |
$k=11$ $p=8$ |
$k=12$ $p=8$ |
|---|---|---|---|---|---|---|
| $\alpha_0$ | -0.113750000 | 0.116262 | -0.099779 | -0.082009 | 0.080469 | 0.100883789 |
| $\alpha_1$ | 0.622546875 | -0.725523 | 0.714959 | 0.508000 | -0.575303 | -0.740019531 |
| $\alpha_2$ | -1.546562500 | 1.989430 | -2.270547 | -1.279249 | 1.713125 | 2.389287109 |
| $\alpha_3$ | 2.625796875 | -3.411867 | 4.423000 | 1.644664 | -2.682900 | -4.648496094 |
| $\alpha_4$ | -3.840625000 | 4.722578 | -6.522184 | -0.976564 | 2.127422 | 6.609213867 |
| $\alpha_5$ | 4.783265625 | -6.010945 | 8.509996 | -0.581992 | 0.127051 | -8.080507813 |
| $\alpha_6$ | -4.754062500 | 6.455258 | -9.798625 | 3.006525 | -3.453203 | 9.064160156 |
| $\alpha_7$ | 2.223390625 | -5.496664 | 8.967422 | -5.764508 | 7.356777 | -9.414804688 |
| $\alpha_8$ | 2.361473 | -6.433865 | 6.906183 | -10.130703 | 9.761733398 | |
| $\alpha_9$ | 2.509623 | -5.806164 | 9.613408 | -10.146035156 | ||
| $\alpha_{10}$ | 2.425114 | -6.737109 | 9.056708984 | |||
| $\alpha_{11}$ | 2.560967 | -6.470136719 | ||||
| $\alpha_{12}$ | 2.518012695 |
1. Consistency
The output of stabregion4 -df Varah7 gives:
Varah7, p=5, k=7, l=1
-0.1138
0.6225
-1.5466
2.6258
-3.8406
4.7833
-4.7541
2.2234
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
1.0000
rho_0 0.000000000
rho_1 0.000000000 <-----
rho_2 -0.000000000 <-----
rho_3 0.000018740 <-----
rho_4 0.000074961 <-----
rho_5 0.000159291 <-----
rho_6 -0.152846199 <-----
The non-zero error constants are due to rounding.
Varah8, p=6, k=8, l=1
0.1163
-0.7255
1.9894
-3.4119
4.7226
-6.0109
6.4553
-5.4967
2.3615
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
1.0000
rho_0 0.000002000 <-----
rho_1 0.000002964 <-----
rho_2 0.000008681 <-----
rho_3 0.000037194 <-----
rho_4 0.000017380 <-----
rho_5 -0.000102595 <-----
rho_6 -0.000303929 <-----
rho_7 -0.147744074 <-----
Varah9, p=7, k=9, l=1
-0.0998
0.7150
-2.2705
4.4230
-6.5222
8.5100
-9.7986
8.9674
-6.4339
2.5096
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
1.0000
rho_0 -0.000000000 <-----
rho_1 -0.000000000 <-----
rho_2 0.000001594 <-----
rho_3 0.000022580 <-----
rho_4 0.000087795 <-----
rho_5 0.000203550 <-----
rho_6 0.000344146 <-----
rho_7 0.000478091 <-----
rho_8 -0.122145352 <-----
Varah10, p=7, k=10, l=1
-0.0820
0.5080
-1.2792
1.6447
-0.9766
-0.5820
3.0065
-5.7645
6.9062
-5.8062
2.4251
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
1.0000
rho_0 0.000000000 <-----
rho_1 0.000000000 <-----
rho_2 -0.000001649 <-----
rho_3 0.000009347 <-----
rho_4 0.000148996 <-----
rho_5 0.000610968 <-----
rho_6 0.001504919 <-----
rho_7 0.002694619 <-----
rho_8 -0.211789593 <-----
Varah11, p=8, k=11, l=1
0.0805
-0.5753
1.7131
-2.6829
2.1274
0.1271
-3.4532
7.3568
-10.1307
9.6134
-6.7371
2.5610
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
1.0000
rho_0 0.000001000 <-----
rho_1 0.000002343 <-----
rho_2 0.000010152 <-----
rho_3 0.000196280 <-----
rho_4 0.000902003 <-----
rho_5 0.002305171 <-----
rho_6 0.004096907 <-----
rho_7 0.005605721 <-----
rho_8 0.006196791 <-----
rho_9 -0.188687540 <-----
Varah12, p=8, k=12, l=1
0.1009
-0.7400
2.3893
-4.6485
6.6092
-8.0805
9.0642
-9.4148
9.7617
-10.1460
9.0567
-6.4701
2.5180
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
1.0000
rho_0 -0.000000003 <-----
rho_1 -0.000000008 <-----
rho_2 -0.000000034 <-----
rho_3 0.000165369 <-----
rho_4 0.000992582 <-----
rho_5 0.003184829 <-----
rho_6 0.007403980 <-----
rho_7 0.014063941 <-----
rho_8 0.023046591 <-----
rho_9 -0.303630714 <-----
2. Parasitic roots
We will clearly see that the Varah formulas have dangerously high parasitic roots.
Varah7
parasitic roots of Varah7
nr real imag abs 1-th root
0 1.00000000 0.00000000 1.00000000 1.00000000
1 -0.24634479 -0.65178270 0.69678292 0.69678292
2 -0.24634479 0.65178270 0.69678292 0.69678292
3 0.34897127 0.56510842 0.66417503 0.66417503
4 0.34897127 -0.56510842 0.66417503 0.66417503
5 0.46647547 0.14587085 0.48875113 0.48875113
6 0.46647547 -0.14587085 0.48875113 0.48875113
By design all Varah methods have only zero as sole roots at infinity. For example, for Varah7:
radius at infinity of Varah7
nr real imag abs 1-th root
0 0.00000000 0.00000000 0.00000000 0.00000000
1 0.00000000 0.00000000 0.00000000 0.00000000
2 0.00000000 0.00000000 0.00000000 0.00000000
3 0.00000000 0.00000000 0.00000000 0.00000000
4 0.00000000 0.00000000 0.00000000 0.00000000
5 0.00000000 0.00000000 0.00000000 0.00000000
6 0.00000000 -0.00000000 0.00000000 0.00000000
Varah8
parasitic roots of Varah8
nr real imag abs 1-th root
0 0.99999800 0.00000000 0.99999800 0.99999800
1 0.30536203 -0.71462628 0.77713364 0.77713364
2 0.30536203 0.71462628 0.77713364 0.77713364
3 -0.34347859 0.69675076 0.77681347 0.77681347
4 -0.34347859 -0.69675076 0.77681347 0.77681347
5 0.46269176 -0.26123688 0.53134581 0.53134581
6 0.46269176 0.26123688 0.53134581 0.53134581
7 0.47849372 0.00000000 0.47849372 0.47849372
Varah9
parasitic roots of Varah9
nr real imag abs 1-th root
0 1.00000000 -0.00000000 1.00000000 1.00000000
1 0.22279293 -0.85616613 0.88467911 0.88467911
2 0.22279293 0.85616613 0.88467911 0.88467911
3 -0.36246291 -0.71938800 0.80554234 0.80554234
4 -0.36246291 0.71938800 0.80554234 0.80554234
5 0.44666741 -0.36226909 0.57510927 0.57510927
6 0.44666741 0.36226909 0.57510927 0.57510927
7 0.47484151 -0.10590790 0.48650894 0.48650894
8 0.47484151 0.10590790 0.48650894 0.48650894
Varah10
parasitic roots of Varah10
nr real imag abs 1-th root
0 1.00000000 0.00000000 1.00000000 1.00000000
1 0.30875638 -0.81372522 0.87033284 0.87033284
2 0.30875638 0.81372522 0.87033284 0.87033284
3 -0.24385210 0.76988452 0.80758035 0.80758035
4 -0.24385210 -0.76988452 0.80758035 0.80758035
5 -0.71900771 -0.00000000 0.71900771 0.71900771
6 0.48630477 0.34866237 0.59837930 0.59837930
7 0.48630477 -0.34866237 0.59837930 0.59837930
8 0.50538574 -0.10234409 0.51564432 0.51564432
9 0.50538574 0.10234409 0.51564432 0.51564432
Varah11
parasitic roots of Varah11
nr real imag abs 1-th root
0 0.99999900 0.00000000 0.99999900 0.99999900
1 0.23641736 -0.94102290 0.97026659 0.97026659
2 0.23641736 0.94102290 0.97026659 0.97026659
3 -0.28365549 -0.82745981 0.87472863 0.87472863
4 -0.28365549 0.82745981 0.87472863 0.87472863
5 -0.72761615 0.00000000 0.72761615 0.72761615
6 0.47293726 0.43778137 0.64445495 0.64445495
7 0.47293726 -0.43778137 0.64445495 0.64445495
8 0.50103277 -0.18679251 0.53471981 0.53471981
9 0.50103277 0.18679251 0.53471981 0.53471981
10 0.50484286 -0.00000000 0.50484286 0.50484286
Varah12
parasitic roots of Varah12
nr real imag abs 1-th root
0 1.00000000 0.00000000 1.00000000 1.00000000
1 0.31945233 -0.87188659 0.92856665 0.92856665
2 0.31945233 0.87188659 0.92856665 0.92856665
3 -0.68490651 -0.62630897 0.92809475 0.92809475
4 -0.68490651 0.62630897 0.92809475 0.92809475
5 -0.15767709 0.84680528 0.86136012 0.86136012
6 -0.15767709 -0.84680528 0.86136012 0.86136012
7 0.51151246 0.41736182 0.66017868 0.66017868
8 0.51151246 -0.41736182 0.66017868 0.66017868
9 0.53030903 0.17960482 0.55989781 0.55989781
10 0.53030903 -0.17960482 0.55989781 0.55989781
11 0.53216050 -0.00000000 0.53216050 0.53216050
3. Stability region
Below are all the stability regions for Varah's formulas.
for i in `seq 7 12`; do stabregion4 -f Varah$i -oj -r750 | grep -Fv '[0.00000000, 0.00000000],' >> /tmp/Varah7-12v; done
4. Comparison
For comparison we show the Widlund-wedge angle α, Widlung-distance δ for
- Varah 1978) formulas
- Original Tendler 1973) formulas
- Extended Tendler formulas
- Tischer formulas
We list the values given by Varah (1978) and not the computed ones from above. We assume that the difference is due to rounding the coefficients.
| Formula | Order p | α | δ | α(Tdl) | δ(Tdl) | α(eTdl) | δ(eTdl) | α(Tischer) | δ(Tischer) |
|---|---|---|---|---|---|---|---|---|---|
| Varah7 | 5 | 72° | 1.6 | 77.5° | 1.4 | 77.8° | 0.4 | 86.6° | 0.04 |
| Varah8 | 6 | 55° | 4.5 | 63.2° | 2.9 | 71.6° | 1.0 | 76.3° | 0.3 |
| Varah9 | 7 | 29° | 11.2 | 33.5° | 10.2 | 55.1° | 3.9 | 57.7° | 1.0 |
| Varah10 | 7 | 39° | 6.6 | 33.5° | 10.2 | 55.1° | 3.9 | 57.7° | 1.0 |
| Varah11 | 8 | 2° | 15.5 | n/a | n/a | n/a | 15.1 | 22.1° | 2.5 |
| Varah12 | 8 | 30° | 18.4 | n/a | n/a | n/a | 15.1 | 22.1° | 2.5 |
Outlook: It would be interesting to create extended Tendler-like formulas based on Varah's formulas. Tendler's and eTendler-formulas are based on BDF.