Construction of New Generalizations of Wynn's Epsilon and Rho Algorithm by Solving Finite Difference Equations in the Transformation Order

TL;DR

This paper introduces new sequence transformation algorithms that generalize Wynn's epsilon and rho methods, unifying existing algorithms and demonstrating improved convergence acceleration through numerical experiments.

Contribution

The paper develops a new class of sequence transformations that encompass Wynn's epsilon, rho, and Osada's generalized rho algorithms, offering a unified framework.

Findings

Algorithms effectively accelerate convergence of various sequences.

Numerical results show improved performance over existing methods.

The new transformations handle both convergent and divergent series.

Abstract

We construct new sequence transformations based on Wynn's epsilon and rho algorithms. The recursions of the new algorithms include the recursions of Wynn's epsilon and rho algorithm and of Osada's generalized rho algorithm as special cases. We demonstrate the performance of our algorithms numerically by applying them to some linearly and logarithmically convergent sequences as well as some divergent series.

Tables

Table 1. Table 1: Numerical results of example 4.1

$n$	$\left|S_n-S\right|$	$\left|{\tilde{T}}_{2\lfloor \left(n-1\right)/2\rfloor }^{\left(n-1-2\lfloor \left(n-1\right)/2\rfloor \right)}-S\right|$	$\left|{ϵ}_{2\lfloor n/2\rfloor }^{\left(n-2\lfloor n/2\rfloor \right)}-S\right|$	$\left|{\vartheta }_{2\lfloor n/3\rfloor }^{\left(n-3\lfloor n/3\rfloor \right)}-S\right|$	$\left|{𝒥}_{\lfloor n/3\rfloor }^{\left(n-3\lfloor n/3\rfloor \right)}-S\right|$
	Eq. 4.4	Eq. 3.20	Eq. 2.6	Eq. 2.18	Eq. 2.23
7	0.059	6.768$\times {10}^{-7}$	1.437$\times {10}^{-6}$	5.702$\times {10}^{-7}$	1.527$\times {10}^{-8}$
8	0.052	3.485$\times {10}^{-8}$	1.518$\times {10}^{-7}$	1.961$\times {10}^{-7}$	7.793$\times {10}^{-9}$
9	0.048	5.907$\times {10}^{-9}$	3.807$\times {10}^{-8}$	8.881$\times {10}^{-10}$	6.034$\times {10}^{-10}$
10	0.043	1.920$\times {10}^{-9}$	4.402$\times {10}^{-9}$	2.764$\times {10}^{-10}$	2.115$\times {10}^{-10}$
11	0.040	8.728$\times {10}^{-11}$	1.042$\times {10}^{-9}$	9.280$\times {10}^{-11}$	7.530$\times {10}^{-11}$
12	0.037	4.244$\times {10}^{-11}$	1.282$\times {10}^{-10}$	1.171$\times {10}^{-13}$	4.954$\times {10}^{-13}$
13	0.034	1.396$\times {10}^{-11}$	2.909$\times {10}^{-11}$	1.665$\times {10}^{-14}$	3.331$\times {10}^{-14}$
14	0.032	3.936$\times {10}^{-13}$	3.745$\times {10}^{-12}$	1.266$\times {10}^{-14}$	2.331$\times {10}^{-15}$

Table 2. Table 2: Numerical results of example 4.2

$n$	$\left|S_n-S\right|$	$\left|{\tilde{T}}_{2\lfloor \left(n-1\right)/2\rfloor }^{\left(n-1-2\lfloor \left(n-1\right)/2\rfloor \right)}-S\right|$	$\left|{ϵ}_{2\lfloor n/2\rfloor }^{\left(n-2\lfloor n/2\rfloor \right)}-S\right|$	$\left|{\vartheta }_{2\lfloor n/3\rfloor }^{\left(n-3\lfloor n/3\rfloor \right)}-S\right|$	$\left|{𝒥}_{\lfloor n/3\rfloor }^{\left(n-3\lfloor n/3\rfloor \right)}-S\right|$
	Eq. 4.6	Eq. 3.20	Eq. 2.6	Eq. 2.18	Eq. 2.23
13	9.657$\times {10}^{-3}$	7.687$\times {10}^{-4}$	3.637$\times {10}^{-6}$	3.533$\times {10}^{-6}$	9.972$\times {10}^{-10}$
14	7.312$\times {10}^{-3}$	1.226$\times {10}^{-4}$	1.272$\times {10}^{-6}$	3.560$\times {10}^{-6}$	2.092$\times {10}^{-10}$
15	5.552$\times {10}^{-3}$	3.558$\times {10}^{-4}$	5.260$\times {10}^{-7}$	3.561$\times {10}^{-6}$	2.026$\times {10}^{-9}$
16	4.228$\times {10}^{-3}$	4.880$\times {10}^{-5}$	1.860$\times {10}^{-7}$	3.585$\times {10}^{-6}$	2.479$\times {10}^{-10}$
17	3.227$\times {10}^{-3}$	1.608$\times {10}^{-5}$	7.597$\times {10}^{-8}$	3.573$\times {10}^{-6}$	6.144$\times {10}^{-11}$
18	2.469$\times {10}^{-3}$	1.990$\times {10}^{-5}$	2.769$\times {10}^{-8}$	3.582$\times {10}^{-6}$	3.903$\times {10}^{-11}$
19	1.892$\times {10}^{-3}$	8.434$\times {10}^{-6}$	1.025$\times {10}^{-8}$	1.837$\times {10}^{-6}$	5.445$\times {10}^{-11}$
20	1.453$\times {10}^{-3}$	4.816$\times {10}^{-6}$	5.256$\times {10}^{-9}$	5.225$\times {10}^{-9}$	4.459$\times {10}^{-11}$

Table 3. Table 3: Numerical results of example 4.3

$n$	$\left|S_n-S\right|$	$\left|{\tilde{T}}_{2\lfloor \left(n-1\right)/2\rfloor }^{\left(n-1-2\lfloor \left(n-1\right)/2\rfloor \right)}-S\right|$	$\left|{\rho }_{2\lfloor n/2\rfloor }^{\left(n-2\lfloor n/2\rfloor \right)}-S\right|$	$\left|{\vartheta }_{2\lfloor n/3\rfloor }^{\left(n-3\lfloor n/3\rfloor \right)}-S\right|$	$\left|{𝒥}_{\lfloor n/3\rfloor }^{\left(n-3\lfloor n/3\rfloor \right)}-S\right|$
	LABEL:{logex1}	Eq. 3.20	Eq. 2.12	Eq. 2.18	Eq. 2.23
3	0.120	4.051$\times {10}^{-4}$	2.942$\times {10}^{-4}$	4.051$\times {10}^{-4}$	4.051$\times {10}^{-4}$
4	0.097	3.190$\times {10}^{-4}$	6.774$\times {10}^{-5}$	3.190$\times {10}^{-4}$	3.190$\times {10}^{-4}$
5	0.081	6.518$\times {10}^{-4}$	7.156$\times {10}^{-6}$	1.875$\times {10}^{-4}$	1.875$\times {10}^{-4}$
6	0.070	2.724$\times {10}^{-5}$	8.286$\times {10}^{-8}$	2.391$\times {10}^{-4}$	2.393$\times {10}^{-4}$
7	0.061	1.521$\times {10}^{-8}$	5.091$\times {10}^{-9}$	2.221$\times {10}^{-5}$	2.238$\times {10}^{-5}$
8	0.055	1.188$\times {10}^{-8}$	1.870$\times {10}^{-9}$	4.735$\times {10}^{-6}$	4.822$\times {10}^{-6}$
9	0.049	5.295$\times {10}^{-9}$	6.288$\times {10}^{-10}$	2.316$\times {10}^{-7}$	1.974$\times {10}^{-7}$
10	0.045	6.931$\times {10}^{-10}$	1.289$\times {10}^{-11}$	4.116$\times {10}^{-8}$	7.269$\times {10}^{-9}$

Table 4. Table 4: Numerical results of example 4.4

$n$	$\left|S_n-S\right|$	$\left|{\tilde{T}}_{2\lfloor \left(n-1\right)/2\rfloor }^{\left(n-1-2\lfloor \left(n-1\right)/2\rfloor \right)}-S\right|$	$\left|{\overline{\rho }}_{2\lfloor n/2\rfloor }^{\left(n-2\lfloor n/2\rfloor \right)}-S\right|$	$\left|{\vartheta }_{2\lfloor n/3\rfloor }^{\left(n-3\lfloor n/3\rfloor \right)}-S\right|$	$|{𝒥}_{\lfloor n/3\rfloor }^{\left(n-3\lfloor n/3\rfloor \right)}-S|$
	Eq. 4.11	Eq. 3.20	Eq. 2.14, $\theta =1/2$	Eq. 2.18	Eq. 2.23
14	0.073	1.029$\times {10}^{-8}$	8.637$\times {10}^{-12}$	8.733$\times {10}^{-10}$	2.680$\times {10}^{-9}$
15	0.071	1.201$\times {10}^{-10}$	1.038$\times {10}^{-12}$	5.325$\times {10}^{-10}$	2.565$\times {10}^{-9}$
16	0.069	5.704$\times {10}^{-11}$	2.517$\times {10}^{-12}$	5.542$\times {10}^{-10}$	1.921$\times {10}^{-9}$
17	0.067	1.081$\times {10}^{-10}$	9.550$\times {10}^{-13}$	1.299$\times {10}^{-9}$	9.512$\times {10}^{-9}$
18	0.065	1.606$\times {10}^{-11}$	1.384$\times {10}^{-11}$	5.174$\times {10}^{-10}$	2.457$\times {10}^{-9}$
19	0.063	1.333$\times {10}^{-10}$	6.473$\times {10}^{-13}$	1.086$\times {10}^{-9}$	7.667$\times {10}^{-9}$
20	0.062	1.583$\times {10}^{-11}$	3.970$\times {10}^{-13}$	1.614$\times {10}^{-6}$	7.311$\times {10}^{-8}$
21	0.060	1.591$\times {10}^{-12}$	6.120$\times {10}^{-13}$	2.152$\times {10}^{-7}$	1.540$\times {10}^{-8}$

Table 5. Table 5: Numerical results of example 4.5

$n$	$\left|S_n-S\right|$	$\left|{\tilde{T}}_{2\lfloor \left(n-1\right)/2\rfloor }^{\left(n-1-2\lfloor \left(n-1\right)/2\rfloor \right)}-S\right|$	$\left|{\overline{\rho }}_{2\lfloor n/2\rfloor }^{\left(n-2\lfloor n/2\rfloor \right)}-S\right|$	$\left|{\vartheta }_{2\lfloor n/3\rfloor }^{\left(n-3\lfloor n/3\rfloor \right)}-S\right|$	$\left|{𝒥}_{\lfloor n/3\rfloor }^{\left(n-3\lfloor n/3\rfloor \right)}-S\right|$
	Eq. 4.12	Eq. 3.20	Eq. 2.14, $\theta =1/2$	Eq. 2.18	Eq. 2.23
14	0.149	4.688$\times {10}^{-9}$	3.580$\times {10}^{-10}$	4.073$\times {10}^{-5}$	1.149$\times {10}^{-7}$
15	0.144	3.321$\times {10}^{-9}$	1.266$\times {10}^{-10}$	4.362$\times {10}^{-5}$	1.068$\times {10}^{-8}$
16	0.139	1.443$\times {10}^{-8}$	2.067$\times {10}^{-10}$	2.654$\times {10}^{-5}$	8.333$\times {10}^{-10}$
17	0.135	1.086$\times {10}^{-9}$	1.338$\times {10}^{-10}$	5.316$\times {10}^{-8}$	4.021$\times {10}^{-9}$
18	0.131	3.993$\times {10}^{-9}$	3.329$\times {10}^{-11}$	8.259$\times {10}^{-8}$	2.607$\times {10}^{-9}$
19	0.128	4.042$\times {10}^{-10}$	8.087$\times {10}^{-11}$	5.279$\times {10}^{-8}$	7.819$\times {10}^{-10}$
20	0.125	2.013$\times {10}^{-9}$	6.629$\times {10}^{-11}$	5.270$\times {10}^{-8}$	1.851$\times {10}^{-8}$
21	0.122	1.075$\times {10}^{-11}$	8.181$\times {10}^{-11}$	5.275$\times {10}^{-8}$	1.858$\times {10}^{-8}$

Table 6. Table 6: Summation of the asymptotic series F 0 2 ​ ( 1 , 1 ; − 1 / z ) = z ​ e z ​ E 1 ​ ( z ) subscript subscript 𝐹 0 2 1 1 1 𝑧 𝑧 superscript 𝑒 𝑧 subscript 𝐸 1 𝑧 {}_{2}F_{0}(1,1;-1/z)=ze^{z}E_{1}(z) for z = 3 𝑧 3 z=3

$n$	$\left|S_n-S\right|$	$\left|{\tilde{T}}_{2\lfloor \left(n-1\right)/2\rfloor }^{\left(n-1-2\lfloor \left(n-1\right)/2\rfloor \right)}-S\right|$	$\left|{\epsilon }_{2\lfloor n/2\rfloor }^{\left(n-2\lfloor n/2\rfloor \right)}-S\right|$	$\left|{\vartheta }_{2\lfloor n/3\rfloor }^{\left(n-3\lfloor n/3\rfloor \right)}-S\right|$	$\left|{𝒥}_{\lfloor n/3\rfloor }^{\left(n-3\lfloor n/3\rfloor \right)}-S\right|$
	Eq. 4.17, z=3	Eq. 3.20	Eq. 2.6	Eq. 2.18	Eq. 2.23
14	1.504$\times {10}^4$	2.587$\times {10}^{-8}$	9.326$\times {10}^{-7}$	5.357$\times {10}^{-11}$	4.735$\times {10}^{-11}$
15	7.609$\times {10}^4$	3.287$\times {10}^{-9}$	7.006$\times {10}^{-7}$	1.379$\times {10}^{-11}$	6.865$\times {10}^{-12}$
16	4.100$\times {10}^5$	1.568$\times {10}^{-9}$	2.877$\times {10}^{-7}$	4.047$\times {10}^{-11}$	6.424$\times {10}^{-11}$
17	2.344$\times {10}^6$	3.844$\times {10}^{-11}$	2.192$\times {10}^{-7}$	3.737$\times {10}^{-10}$	3.744$\times {10}^{-11}$
18	1.418$\times {10}^7$	8.860$\times {10}^{-9}$	9.445$\times {10}^{-8}$	4.069$\times {10}^{-11}$	1.258$\times {10}^{-12}$
19	9.048$\times {10}^7$	2.204$\times {10}^{-9}$	7.289$\times {10}^{-8}$	6.707$\times {10}^{-11}$	2.000$\times {10}^{-10}$
20	6.073$\times {10}^8$	4.265$\times {10}^{-9}$	3.272$\times {10}^{-8}$	3.673$\times {10}^{-13}$	1.386$\times {10}^{-12}$
21	4.277$\times {10}^9$	2.101$\times {10}^{-10}$	2.552$\times {10}^{-8}$	9.726$\times {10}^{-14}$	1.785$\times {10}^{-12}$

Table 7. Table 7: Summation of the asymptotic series F 0 2 ​ ( 1 , 1 ; − 1 / z ) = z ​ e z ​ E 1 ​ ( z ) subscript subscript 𝐹 0 2 1 1 1 𝑧 𝑧 superscript 𝑒 𝑧 subscript 𝐸 1 𝑧 {}_{2}F_{0}(1,1;-1/z)=ze^{z}E_{1}(z) for z = 1 / 2 𝑧 1 2 z=1/2

$n$	$\left|S_n-S\right|$	$\left|{\tilde{T}}_{2\lfloor \left(n-1\right)/2\rfloor }^{\left(n-1-2\lfloor \left(n-1\right)/2\rfloor \right)}-S\right|$	$\left|{\epsilon }_{2\lfloor n/2\rfloor }^{\left(n-2\lfloor n/2\rfloor \right)}-S\right|$	$\left|{\vartheta }_{2\lfloor n/3\rfloor }^{\left(n-3\lfloor n/3\rfloor \right)}-S\right|$	$\left|{𝒥}_{\lfloor n/3\rfloor }^{\left(n-3\lfloor n/3\rfloor \right)}-S\right|$
	Eq. 4.17, z=1/2	Eq. 3.20	Eq. 2.6	Eq. 2.18	Eq. 2.23
14	1.379$\times {10}^{15}$	3.012$\times {10}^{-4}$	2.065$\times {10}^{-3}$	3.963$\times {10}^{-5}$	1.862$\times {10}^{-5}$
15	4.147$\times {10}^{16}$	4.489$\times {10}^{-5}$	3.346$\times {10}^{-3}$	4.578$\times {10}^{-7}$	1.052$\times {10}^{-6}$
16	1.330$\times {10}^{18}$	1.498$\times {10}^{-5}$	1.257$\times {10}^{-3}$	1.314$\times {10}^{-7}$	4.844$\times {10}^{-7}$
17	4.529$\times {10}^{19}$	1.364$\times {10}^{-6}$	1.987$\times {10}^{-3}$	4.963$\times {10}^{-7}$	2.710$\times {10}^{-7}$
18	1.633$\times {10}^{21}$	7.763$\times {10}^{-6}$	7.876$\times {10}^{-4}$	4.044$\times {10}^{-7}$	2.207$\times {10}^{-6}$
19	6.214$\times {10}^{22}$	6.418$\times {10}^{-6}$	1.218$\times {10}^{-3}$	1.425$\times {10}^{-7}$	4.204$\times {10}^{-7}$
20	2.489$\times {10}^{24}$	2.846$\times {10}^{-6}$	5.052$\times {10}^{-4}$	8.537$\times {10}^{-8}$	2.591$\times {10}^{-7}$
21	1.047$\times {10}^{26}$	3.780$\times {10}^{-6}$	7.662$\times {10}^{-4}$	3.127$\times {10}^{-8}$	3.538$\times {10}^{-7}$

Equations

$$\mathcal{T}:\{S_n\}\to \{T_n\},n\in {\mathbb{N}}_0,$$

$$\underset{n\to \infty }{\lim }\frac{T_n-S}{S_n-S}=0.$$

$$\underset{n\to \infty }{\lim }\frac{S_{n+1}-S}{S_n-S}=\rho .$$

$$T_{k}^{(n)}\;=\;T_{k}^{(n)}\bigl{(}S_{n},S_{n+1},\ldots,S_{n+\ell(k)}\bigr{)}\,.$$

$$S_n=S+R_n,n\in {\mathbb{N}}_0.$$

$$S_n^{\prime }=S_n-{\tilde{R}}_n=S+R_n-{\tilde{R}}_n=S+R_n^{\prime }.$$

$$S_n=S+C{\lambda }^n,C\ne 0,|\lambda |<1,n\in {\mathbb{N}}_0.$$

$${\mathcal{A}}_1^{(n)}=S_n-\frac{[\Delta S_n{]}^2}{{\Delta }^2{S}_n},n\in {\mathbb{N}}_0,$$

$${\mathcal{A}}_0^{(n)}$$

$${\mathcal{A}}_{k+1}^{(n)}$$

$$S_n=S+\sum _{j=0}^{k-1}{C}_j({\lambda }_j{)}^n,k,n\in {\mathbb{N}}_0.$$

$$e_k(S_n)=\frac{\left|\begin{array}{cccc}{S}_n& S_{n+1}& \dots & S_{n+k}\\ \Delta S_n& \Delta S_{n+1}& \dots & \Delta S_{n+k}\\ ⋮& ⋮& \ddots & ⋮\\ \Delta S_{n+k-1}& \Delta S_{n+k}& \dots & \Delta S_{n+2k-1}\end{array}\right|}{\left|\begin{array}{cccc}1& 1& \dots & 1\\ \Delta S_n& \Delta S_{n+1}& \dots & \Delta S_{n+k}\\ ⋮& ⋮& \ddots & ⋮\\ \Delta S_{n+k-1}& \Delta S_{n+k}& \dots & \Delta S_{n+2k-1}\end{array}\right|},k,n\in {\mathbb{N}}_0.$$

$$f_n(z)=\sum _{\nu =0}^n{\gamma }_{\nu }{z}^{\nu },n\in {\mathbb{N}}_0,$$

$$e_{k}\bigl{(}f_{n}(z)\bigr{)}\;=\;[k+n/k]_{f}\,(z)\,,\qquad k,n\in\mathbb{N}_{0}\,.$$

$$H_0(u_n)=1,H_1(u_n)=u_n,n\in {\mathbb{N}}_0,$$

$$H_{k+2}(u_n)H_k(u_{n+2})=H_{k+1}(u_n)H_{k+1}(u_{n+2})-{\left[H_{k+1}(u_{n+1})\right]}^2,k,n\in {\mathbb{N}}_0.$$

$${\epsilon }_{-1}^{(n)}$$

$${\epsilon }_{k+1}^{(n)}$$

$${\epsilon }_{2k}^{(n)}=e_k(S_n),k,n\in {\mathbb{N}}_0,$$

$${\epsilon }_{2k+1}^{(n)}=1/e_k(\Delta S_n),k,n\in {\mathbb{N}}_0,$$

$${\epsilon }_{2k}^{(n)}=[k+n/k{]}_f(z),k,n\in {\mathbb{N}}_0.$$

$${\rho }_{-1}^{(n)}$$

$${\rho }_{k+1}^{(n)}$$

$$0<x_0<x_1<\cdots <x_m<x_{m+1}<\cdots ,$$

$$\underset{n\to \infty }{\lim }{x}_n=\infty .$$

$${\rho }_{-1}^{(n)}$$

$${\rho }_{k+1}^{(n)}$$

$$S_n=S+(n+\beta {)}^{-\theta }\sum _{j=0}^{\infty }{c}_j/(n+\beta {)}^j,n\in {\mathbb{N}}_0.$$

$${\bar{\rho }}_{-1}^{(n)}$$

$${\bar{\rho }}_{k+1}^{(n)}$$

$${\bar{\rho}}_{2k}^{(n)}\,-\,S\;=\;\mathrm{O}\bigl{(}n^{-\theta-2k}\bigr{)}\,,\qquad n\to\infty\,.$$

$$T_n=\frac{[{\Delta }^2{S}_n][{\Delta }^2{S}_{n+1}]}{[\Delta S_{n+1}][{\Delta }^2{S}_{n+1}]-[\Delta S_{n+2}][{\Delta }^2{S}_n]}-1,n\in {\mathbb{N}}_0.$$

$$\theta =T_n+O(1/n^2),n\to \infty ,$$

$${\vartheta }_{-1}^{(n)}$$

Full text

††journal: Numerical Algorithms

We dedicate this article to the memory of Peter Wynn (1931 - 2017)

1. 1 Introduction

2. 2 Review of Wynn’s epsilon and rho algorithm and of related transformations

3. 2.1 The Shanks transformation and Wynn’s epsilon algorithm

4. 2.2 Wynn’s rho algorithm

5. 2.3 Brezinski’s theta algorithm and its iteration

6. 3 Difference equations with respect to the transformation order

7. 3.1 General Considerations

8. 3.2 First order difference equations

9. 3.3 Second order difference equations

10. 3.4 Initial conditions for the new transformation

11. 4 Numerical examples

12. 4.1 Asymptotic model sequences

13. 4.2 Linearly convergent sequences

14. 4.3 Logarithmically convergent sequences

15. 4.4 Divergent series

16. 5 Conclusions and discussions

17. A Reduced Bessel functions

18. B Levin’s transformation

19. B.1 The general Levin transformation

20. B.2 Levin’s $u$ transformation

21. B.3 Levin’s $v$ transformation

1 Introduction

Sequences and series are extremely important mathematical tools. They appear naturally in many methods for solving differential and integral equations, in discretization methods, in quadrature schemes, in perturbation techniques, or in the evaluation of special functions. Unfortunately, the resulting sequences or series often do not converge fast enough to be practically useful, or they even diverge, which means that summation techniques have to be employed to obtain numerically useful information. In such situations, it is often helpful to employ so-called sequence transformations.

Formally, a sequence transformation $\mathcal{T}$ is simply a map

[TABLE]

which transforms a given sequence $\{S_{n}\}$, whose convergence may be unsatisfactory, to a new sequence $\{T_{n}\}$ with hopefully better numerical properties. Of course, there is the minimal requirement that the transformed sequence $\{T_{n}\}$ must have the same (generalized) limit $S$ as the original sequence $\{S_{n}\}$, but otherwise we have a lot of freedom. It is common to say that $\mathcal{T}$ accelerates convergence if $\{T_{n}\}$ converges more rapidly to $S$ than $\{S_{n}\}$ according to

[TABLE]

Many practically relevant sequences can be classified according to some simple convergence types. Let us assume that the elements of a sequence $\{S_{n}\}$ with (generalized) limit $S$ satisfy the following criterion resembling the well known ratio test:

[TABLE]

If $0<|\rho|<1$ holds, $\{S_{n}\}$ converges linearly, if $\rho=1$ holds, $\{S_{n}\}$ converges logarithmically, and if $\rho=0$ holds, it converges hyperlinearly. Obviously, $|\rho|>1$ implies divergence. The partial sums $\sum_{k=0}^{n}z^{k}=[1-z^{n+1}]/[1-z]$ of the geometric series converge linearly for $|z|<1$ as $n\to\infty$, and they diverge for $|z|\geq 1$. The partial sums $\sum_{k=0}^{n}(k+1)^{-s}$ of the Dirichlet series $\zeta(s)=\sum_{\nu=0}^{\infty}(\nu+1)^{-s}$ for the Riemann zeta function converge logarithmically for $\Re(s)>1$ (see for example (Weniger and Kirtman, 2003, Eq. (2.6))), and the partial sums $\sum_{k=0}^{n}z^{k}/k!$ of the power series for $\exp(z)$ converge hyperlinearly for $z\in\mathbb{C}$.

A sequence transformation $\mathcal{T}$ corresponds to an infinite set of doubly indexed quantities $T_{k}^{(n)}$ with $k,n\in\mathbb{N}_{0}$. Each $T_{k}^{(n)}$ is computed from a finite subset \bigl{\{}S_{n},S_{n+1},\dots,S_{n+\ell(k)}\bigr{\}} of consecutive elements, where $\ell=\ell(k)$. In our notation, $T_{0}^{(n)}$ always corresponds to an untransformed sequence element $T_{0}^{(n)}=S_{n}$.

The transformation order $k$ is a measure for the complexity of $T_{k}^{(n)}$, and the superscript $n$ corresponds to the minimal index $n$ of the string \bigl{\{}S_{n},S_{n+1},\dots,S_{n+\ell(k)}\bigr{\}} of input data. An increasing value of $k$ implies that the complexity of the transformation process as well as $\ell=\ell(k)$ increases. Thus, the application of a sequence transformation $\mathcal{T}$ to $\{S_{n}\}$ produces for every $k,n\geq\mathbb{N}_{0}$ a new transformed element

[TABLE]

Sequence transformations are at least as old as calculus, and in rudimentary form they are even much older. For the older history of sequence transformations, we recommend a monograph by Brezinski (1991b), which discusses earlier work on continued fractions, Padé approximants, and related topics starting from the 17th century until 1945, as well as an article by Brezinski (2009), which emphasizes pioneers of extrapolation methods. More recent developments are treated in articles by Brezinski (1996, 2000a, 2019) and Brezinski and Redivo Zaglia (2019). A very complete list of older references up to 1991 can be found in a bibliography by Brezinski (1991a).

Ever since the pioneering work of Shanks (1949) and Wynn (1956a), respectively, nonlinear and in general also nonregular sequence transformations have dominated both research and practical applications. For more details, see the monographs by Brezinski (1977, 1978, 1980), Brezinski and Redivo Zaglia (1991), Cuyt and Wuytack (1987), Delahaye (1988), Liem, Lü, and Shih (1995), Marchuk and Shaidurov (1983), Sidi (2003), Walz (1996), and Wimp (1981), or articles by Caliceti, Meyer-Hermann, Ribeca, Surzhykov, and Jentschura (2007), Homeier (2000a), and Weniger (1989). Sequence transformations are also treated in the book by Baker and Graves-Morris (1996) on Padé approximants, which can be viewed to be just a special class of sequence transformations since they convert the partial sums of a power series to a doubly indexed sequence of rational functions.

Practical applications of sequence transformations are emphasized in a book by Bornemann, Laurie, Wagon, and Waldvogel (2004, Appendix A) on extreme digit hunting, in the most recent (third) edition of the book Numerical Recipes (Press et al, 2007) (a review of the treatment of sequence transformations in Numerical Recipes by Weniger (2007) can be downloaded from the Numerical Recipes web page), in a book by Gil, Segura, and Temme (2007) on the evaluation of special functions (compare also the related review articles by Gil, Segura, and Temme (2011) and by Temme (2007)), in the recently published NIST Handbook of Mathematical Functions (Olver et al, 2010, Chapter 3.9 Acceleration of Convergence), and in a recent book by Trefethen (2013) on approximation theory. Readers interested in current research both on and with sequence transformations might look into the proceedings of a recent conference (Brezinski et al, 2010; Weniger, 2010).

Certain sequence transformations are closely related to dynamical systems. For example, Wynn’s epsilon (Wynn, 1956a) and rho (Wynn, 1956b) algorithm can be viewed to be just fully discrete integrable systems (Nagai and Satsuma, 1995; Papageorgiou et al, 1993). New sequence transformations were also derived via this connection with dynamical systems (Brezinski et al, 2012; He et al, 2011; Chang et al, 2018; Sun et al, 2013).

It is the basic assumption of all the commonly used sequence transformations that the elements $S_{n}$ of a given sequence $\{S_{n}\}$ can be partitioned into a (generalized) limit $S$ and a remainder $R_{n}$ according to

[TABLE]

How do sequence transformations accomplish an acceleration of convergence or a summation of a divergent sequence? Let us assume that some sufficiently good approximations $\{\tilde{R}_{n}\}$ to the actual remainders $\{R_{n}\}$ of the elements of the sequence $\{S_{n}\}$ are known. Elimination of these approximations $\tilde{R}_{n}$ yields a new sequence

[TABLE]

If $\tilde{R}_{n}$ is a good approximation to $R_{n}$, the transformed remainder $R^{\prime}_{n}=R_{n}-\tilde{R}_{n}$ vanishes more rapidly than the original remainders $R_{n}$ as $n\to\infty$. At least conceptually, this is what a sequence transformation tries to accomplish, although it may actually use a completely different computational algorithm. Thus, when trying to construct a new sequence transformation, we have to look for arithmetic operations that lead in each step to an improved approximate elimination of the truncation errors.

The idea of viewing a sequence transformation as an approximate elimination procedure for the truncation errors becomes particularly transparent in the case of model sequences whose remainders $R_{n}$ consists of a single exponential term:

[TABLE]

This is an almost trivially simple model problem. Nevertheless, many practically relevant sequences are known where exponential terms of the type of $C\lambda^{n}$ form the dominant parts of their truncation errors.

A short computation shows that the so-called $\Delta^{2}$ formula (see for example (Weniger, 1989, Eq. (5.1-4)))

[TABLE]

is exact for the model sequence (1.7), i.e., we have $\mathcal{A}_{1}^{(n)}=S$. Here, $\Delta$ is the finite difference operator defined by $\Delta f(n)=f(n+1)-f(n)$.

Together with its numerous mathematically, but not necessarily numerically equivalent explicit expressions (see for example (Weniger, 1989, Eqs. (5.1-6) - (5.1-12))), the $\Delta^{2}$ formula (1.8) is commonly attributed to Aitken (1926), although it is in fact much older. Brezinski (1991b, pp. 90 - 91) mentioned that in 1674 Seki Kowa, the probably most famous Japanese mathematician of that period, tried to obtain better approximations to $\pi$ with the help of the $\Delta^{2}$ formula (see also Osada’s article (Osada, 2012, Section 5)), and according to Todd (1962, p. 5), the $\Delta^{2}$ formula was in principle known to Kummer (1837) already in 1837.

The power and practical usefulness of the $\Delta^{2}$ formula (1.8) is obviously limited since it only eliminates a single exponential term. However, the $\Delta^{2}$ formula (1.8) can be iterated by using the output data $\mathcal{A}_{1}^{(n)}$ as input data in the $\Delta^{2}$ formula. Repeating this process yields the following nonlinear recursive scheme (see for example (Weniger, 1989, Eq. (5.1-15))):

[TABLE]

Here, $\Delta$ only acts on the superscript $n$ but not on the subscript $k$ according to $\Delta\mathcal{A}_{k}^{(n)}=\mathcal{A}_{k}^{(n+1)}-\mathcal{A}_{k}^{(n)}$ (if necessary, we write $\Delta_{n}\mathcal{A}_{k}^{(n)}$ for the sake of clarity).

The iterated $\Delta^{2}$ process (1.9) is a fairly powerful sequence transformation (see for example (Weniger, 1989, Table 13-1 on p. 328 or Section 15.2)). In (Weniger, 1989, p. 228), one finds the listing of a FORTRAN 77 program that computes the iterated $\Delta^{2}$ process using a single one-dimensional array.

There is an extensive literature on both the $\Delta^{2}$ formula (1.8) and its iteration (1.9) (see for example (Brezinski and Redivo Zaglia, 2019; Weniger, 2000a) and references therein). But it is even more important that Aitken paved the way for the Shanks transformation (Shanks, 1955) and for Wynn’s epsilon algorithm (Wynn, 1956a), which permits a convenient recursive computation of the Shanks transformation and thus also of Padé approximants. In 1956, Wynn also derived his so-called rho algorithm (Wynn, 1956b), whose recursive scheme closely resembles that of the epsilon algorithm. These well known facts will be reviewed in more detail in Section 2.

The aim of our article is the clarification of the relationships among known generalizations of Wynn’s epsilon (Wynn, 1956a) and rho (Wynn, 1956b) algorithms, and the construction of new generalizations. We accomplish this by considering finite difference equations in the transformation order $k$. To the best of our knowledge, this is a novel approach. It leads to a better understanding of the epsilon and rho algorithm and of many of their generalizations. In addition, our approach also produces several new sequence transformations.

It would be overly optimistic to assume that this article, which is the first one to describe and apply our novel approach, could provide definite answers to all questions occurring in this context. We also do not think that the new transformations, which are presented in this article, exhaust the potential of our novel approach.

The details of our novel approach will be discussed in Section 3. In Section 4, we will show the performance of our algorithms by some numerical examples. Finally, Section 5 is devoted to discussions and conclusions.

2 Review of Wynn’s epsilon and rho algorithm and of related

transformations

2.1 The Shanks transformation and Wynn’s epsilon algorithm

An obvious generalization of the model sequence (1.7), which leads to the $\Delta^{2}$ formula (1.8), is the following one, which contains $k$ exponential terms instead of one:

[TABLE]

It is usually assumed that the $\lambda_{j}$ are distinct ($\lambda_{m}\neq\lambda_{n}$ for $m\neq n$), and that they ordered according to magnitude ($|\lambda_{0}|>|\lambda_{1}|>\dots>|\lambda_{k-1}|$).

Although the $\Delta^{2}$ formula (1.8) is by construction exact for the model sequence (1.7), its iteration (1.9) with $k>1$ is not exact for the model sequence (2.1). Instead, the Shanks transformation, which is defined by the following ratio of Hankel determinants (Shanks, 1955, Eq. (2)),

[TABLE]

is exact for the model sequence (2.1). Actually, the Shanks transformation had originally been derived in 1941 by Schmidt (1941), but Schmidt’s article had largely been ignored.

Shanks’ discovery of his transformation (2.2) had an enormous impact, in particular since he could show that it produces Padé approximants. Let

[TABLE]

stand for the partial sums of the (formal) power series for some function $f$. Then (Shanks, 1955, Theorem VI on p. 22),

[TABLE]

We use the notation of Baker and Graves-Morris (1996, Eq. (1.2)) for Padé approximants.

The Hankel determinants in Shanks’ transformation (2.2) can be computed recursively via the following non-linear recursion (see for example (Brezinski and Redivo Zaglia, 1991, p. 80):

[TABLE]

This recursive scheme is comparatively complicated. Therefore, Wynn’s discovery of his celebrated epsilon algorithm was a substantial step forwards (Wynn, 1956a, Theorem on p. 91):

[TABLE]

Wynn (1956a, Theorem on p. 91) showed that the epsilon elements with even subscripts give Shanks’ transformation

[TABLE]

whereas the elements with odd subscripts are only auxiliary quantities

[TABLE]

which diverge if the transforms (2.7) converge to the (generalized) limit $S$ of $\{S_{n}\}$.

The epsilon algorithm is a close relative of Aitken’s iterated $\Delta^{2}$ process. We have $\mathcal{A}_{1}^{(n)}=\varepsilon_{2}^{(n)}$, but for $k>1$ we in general have $\mathcal{A}_{k}^{(n)}\neq\varepsilon_{2k}^{(n)}$. Nevertheless, the iterated $\Delta^{2}$ process and the epsilon algorithm often have similar properties in convergence acceleration and summation processes.

Because of the rhombus structure of the four-term recursion (2.6b), it appears that a program for Wynn’s epsilon algorithm would need either a single two-dimensional or at least two one-dimensional arrays. However, Wynn (1965) could show with the help of his moving lozenge technique, which is described in a very detailed way in (Brezinski, 1978, Chapter 4.2.1.2), that a single one-dimensional array plus two auxiliary variables suffice. In (Brezinski, 1978, Chapter 4.3.2), one finds listings of several FORTRAN 66 programs for the epsilon algorithm.

As shown in (Weniger, 1989, Section 4.3), a slight improvement of Wynn’s moving lozenge technique is possible. The listing of a compact FORTRAN 77 program for the epsilon algorithm using this modification can be found in (Weniger, 1989, p. 222). In the book by Press, Teukolsky, Vetterling, and Flannery (2007, p. 213), one finds a translation of this FORTRAN 77 program to C.

It follows from Eqs. 2.4 and 2.7 that the epsilon algorithm produces Padé approximants if its input data $S_{n}$ are the partial sums (2.3) of a power series:

[TABLE]

Consequently, the epsilon algorithm is treated in books on Padé approximants such as the one by Baker and Graves-Morris (1996), but there is also an extensive literature dealing directly with it. According to Wimp (1981, p. 120), over 50 articles on the epsilon algorithm were published by Wynn alone, and at least 30 articles by Brezinski. As a fairly complete source of references, Wimp recommends Brezinski’s first book Brezinski (1977). However, this book was published in 1977, and since then many more books and articles dealing with the epsilon algorithm have been published. Any attempt of providing a reasonably complete bibliography would be beyond the scope of this article.

The epsilon algorithm is not limited to input data that are the partial sums of a (formal) power series. Therefore, it is actually more general than Padé approximants. For example, sequences of vectors can be used as input data. For a relatively recent review, see (Graves-Morris et al, 2000a). Both Shanks’ transformation and the epsilon algorithm are discussed in the NIST Handbook of Mathematical Functions (Olver, Lozier, Boisvert, and Clark, 2010, Chapter 3.9(iv) Shanks’ Transformation).

Several generalizations of the epsilon algorithm are known. In 1972, Brezinski (1972, pp. 72 and 78) introduced what he called his first and second generalization of the epsilon algorithm (see also (Brezinski and Redivo Zaglia, 1991, Chapter 2.6)). Sablonniere (1991, Algorithm A1 on p. 180) and Sedogbo (1990, Algorithm 1 and Theorem 1 on p. 255) introduced other generalizations of the epsilon algorithms that were specially designed to be effective for certain logarithmically convergent interation sequences. Other generalizations were considered by Vanden Broeck and Schwartz (1979, Eq. (2.1)) and by Barber and Hamer (1982, Eqs. (43) - (44)).

2.2 Wynn’s rho algorithm

The iterated $\Delta^{2}$ process (1.9) and the epsilon algorithm (2.6) are powerful accelerators for linearly convergent sequences, and they are also able to sum many alternating divergent series (see for example (Weniger, 1989, Section 15.2) or the rigorous convergence analysis of the summation of the Euler series in (Borghi and Weniger, 2015)). However, both fail in the case of logarithmic convergence (Wynn, 1966, Theorem 12). Fortunately, in 1956b Wynn also derived his rho algorithm (Wynn, 1956b, Eq. (8))

[TABLE]

which is – as for example emphasized by Osada (1996) – an effective accelerator for many logarithmically convergent sequences.

As in the case of the formally almost identical epsilon algorithm, only the elements $\rho_{2k}^{(n)}$ with even subscripts provide approximations to the limit of the input sequence. The elements $\rho_{2k+1}^{(n)}$ with odd subscripts are only auxiliary quantities which diverge if the transformation process converges. Actually, the elements $\rho_{2k}^{(n)}$ correspond to the terminants of an interpolating continued fraction with interpolation points $\{x_{n}\}$ that are extrapolated to infinity (see for example (Cuyt and Wuytack, 1987, Chapter IV.1.4)). Consequently, the interpolation points $\{x_{n}\}$ must be positive, strictly increasing, and unbounded:

[TABLE]

In the vast majority of all applications, Wynn’s rho algorithm (2.10) is used in combination with the interpolation points $x_{n}=n+1$, yielding its standard form (see for example (Weniger, 1989, Eq. (6.2-4))):

[TABLE]

The standard form (2.12) is not able to accelerate the convergence of all logarithmically convergent sequences of interest. This can be demonstrated by considering the following class of logarithmically convergent model sequences:

[TABLE]

Here, $\theta$ is a positive decay parameter and $\beta$ is a positive shift parameter. The elements of numerous practically relevant logarithmically convergent sequences $\{S_{n}\}$ can at least in the asymptotic domain of large indices $n$ be represented by series expansions of that kind.

Osada (1990, Theorem 3.2) showed that the standard form (2.12) accelerates the convergence of sequences of the type of Eq. 2.13 if the decay parameter $\theta$ is a positive integer, but it fails if $\theta$ is non-integral.

If the decay parameter $\theta$ of a sequence of the type of Eq. 2.13 is explicitly known, Osada’s variant of the rho algorithm can be used (Osada, 1990, Eq. (3.1)):

[TABLE]

Osada showed that his variant (2.14) accelerates the convergence of sequences of the type of Eq. 2.13 for known, but otherwise arbitrary, i.e., also for non-integral $\theta>0$. He proved the following asymptotic estimate for fixed $k$ (Osada, 1990, Theorem 4):

[TABLE]

As remarked above, the decay parameter $\theta$ must be explicitly known if Osada’s transformation (2.14) is to be applied to a sequence of the type of Eq. 2.13. An approximation to $\theta$ can be obtained with the help of the following nonlinear transformation, which was first derived in a somewhat disguised form by Drummond (1976, p. 419) and later rederived by Bjørstad, Dahlquist, and Grosse (1981, Eq. (4.1)):

[TABLE]

Bjørstad, Dahlquist, and Grosse (1981, Eq. (4.1)) showed that

[TABLE]

if the input data are the elements of a sequence of the type of Eq. 2.13.

The variants of Wynn’s rho algorithm can also be iterated in the spirit of Aitken’s iterated $\Delta^{2}$ process (1.9), which can also be viewed to be an iteration of $\varepsilon_{2}^{(n)}$. By iterating the expression for $\rho_{2}^{(n)}$ involving unspecified interpolation points $x_{n}$, a rho analog of Aitken’s iterated $\Delta^{2}$ was constructed (Weniger, 1991, Eq. (2.10)). An alternative iteration was derived by Bhowmick, Bhattacharya, and Roy (1989, Eq. (2.25)), which is, however, significantly less efficient than the iteration (Weniger, 1991, Eq. (2,10)) involving unspecified interpolation points (compare (Weniger, 1991, Table I)). This is caused by the fact that Bhowmick, Bhattacharya, and Roy (1989, Eq. (2.25)) had started from the standard form (2.12) and not from the general form (2.10). Osada’s variant (2.14) of the rho algorithm can also be iterated. It yields a recursive scheme (Weniger, 1991, Eq. (2.29)), which was originally derived by Bjørstad, Dahlquist, and Grosse (1981, Eq. (2.4)) who called the resulting algorithm a modified $\Delta^{2}$ formula.

Wynn’s epsilon algorithm (2.6) is undoubtedly the currently most important and most often employed sequence transformation. Wynn’s rho algorithm (2.10) and its standard form (2.12), Osada’s generalized rho algorithm (2.14), as well as the iterations mentioned above are also very important and practically useful sequence transformations. Accordingly, there are numerous references describing successful applications (far too many to be cited here). In addition, these transformations were very consequential from a purely theoretical point of view. They inspired the derivation of numerous other sequence transformations and ultimately also this article.

2.3 Brezinski’s theta algorithm and its iteration

Soon after their derivation it had become clear that epsilon and rho have largely complementary properties. The epsilon algorithm is a powerful accelerator of linear convergence and also frequently works well in the case of divergent series, but it does not accomplish anything substantial in the case of logarithmic convergence. In contrast, the rho algorithm is ineffective in the case of linear convergence and even more so in the case of divergent series, but is able to speed up the convergence of many logarithmically convergent sequences.

This observation raised the question whether it would be possible to modify the recursions of epsilon and/or rho in such a way that the resulting transformation would possess the advantageous features of both epsilon and rho. This goal was achieved in 1971 by Brezinski who suitably modified the epsilon recursion (2.6) to obtain his theta algorithm (Brezinski, 1971a, p. 729) (compare also (Brezinski and Redivo Zaglia, 1991, Chapter 2.3) or (Weniger, 1989, Section 10.1)):

[TABLE]

We again assume that the finite difference operator $\Delta$ only acts on $n$ but not on $k$.

As in the case of Wynn’s epsilon and rho algorithm, only the elements $\vartheta_{2k}^{(n)}$ with even subscripts provide approximations to the (generalized) limit of the sequence $\{S_{n}\}$ to be transformed. The elements $\vartheta_{2k+1}^{(n)}$ with odd subscripts are only auxiliary quantities which diverge if the whole process converges.

Brezinski’s derivation of his theta algorithm was purely experimental, but it was certainly a very successful experiment. Extensive numerical studies performed by Smith and Ford (1979, 1982) showed that among virtually all sequence transformations known at that time only the theta algorithm together with Levin’s $u$ and $v$ transformations (Levin, 1973, Eqs. (59) and (68)) provided consistently good results for a very wide range of test problems. This positive verdict was also confirmed by later numerical studies.

In Brezinski’s second book (Brezinski, 1978, pp. 368 - 370), one finds the listing of a FORTRAN IV program which computes Brezinski’s theta algorithm 2.18 using three one-dimensional arrays, and in (Weniger, 1989, pp. 279 - 281), one finds the listing of a FORTRAN 77 program which computes theta using two one-dimensional arrays.

It is possible to proceed in the spirit of Aitken’s iterated $\Delta^{2}$ process (1.9) and to construct iterations of Brezinski’s theta algorithm (2.18). A natural starting point for such an iteration would be $\vartheta_{2}^{(n)}$ which possesses many alternative expressions (see for example (Weniger, 1989, Eqs. (10.3-1) - (10.3-3)) or (Weniger, 2000a, Eq. (4.10))):

[TABLE]

As in the case of $\mathcal{A}_{2}^{(n)}$, these expressions are obviously mathematically equivalent, but they may differ – possibly substantially – with respect to their numerical stability.

By iterating the explicit expression (2.19) for $\vartheta_{2}^{(n)}$ in the spirit of Aitken’s iterated $\Delta^{2}$ process (1.9), the following nonlinear recursive system was derived (Weniger, 1989, Eq. (10.3-6)):

[TABLE]

In (Weniger, 1989, pp. 284 - 285), one finds the listing a FORTRAN 77 program that computes $\mathcal{J}_{k}^{(n)}$ using a single one-dimensional array.

Based on the explicit expression (2.22) for $\vartheta_{2}^{(n)}$, an alternative recursive scheme for $\mathcal{J}_{k}^{(n)}$ was derived in (Weniger, 2000a, Eq. (4.11)) and used for the prediction of unknown power series coefficients.

The iterated transformation $\mathcal{J}_{k}^{(n)}$ has similar properties as the theta algorithm from which it was derived. Both are very powerful as well as very versatile. $\mathcal{J}_{k}^{(n)}$ is not only an effective accelerator for linear convergence as well as able to sum divergent series, but it is also able to accelerate the convergence of many logarithmically convergent sequences and series (see for example (Sidi, 2002; Weniger, 1989, 1991, 1994, 2000a) and references therein).

As discussed in Appendix B, the special case $\vartheta_{2}^{(n)}=\mathcal{J}_{1}^{(n)}$ of the iterated transformation $\mathcal{J}_{k}^{(n)}$ can be related to the $u$ and $v$ variants of Levin’s general transformation (Levin, 1973, Eq. (22)) defined by Eq. B.2. For Levin’s $u$ and $v$ variants (Levin, 1973, Eqs. (59) and (68)), we use the notation (Weniger, 1989, Eqs. (7.3-5) and (7.3-11)):

[TABLE]

Equations 2.19, B.8 and B.12 imply (see also (Weniger, 1989, Eq. (10.3-4))):

[TABLE]

3 Difference equations with respect to the transformation order

3.1 General Considerations

In this Section, we want to clarify the relationship between known generalizations of Wynn’s epsilon and rho algorithm, and – more important – we also want to derive new generalizations. For that purpose, we pursue a novel approach and look at their recursions from the perspective of finite difference equations with respect to the transformation order $k$.

Our approach is based on the observation that the recursions (2.6b), (2.10b), (2.12b), and (2.14b) for Wynn’s epsilon and rho algorithm and its special cases all possess the following general structure:

[TABLE]

The quantities $F_{k}^{(n)}$, which in general depend on both $k$ and $n$, fully characterize the corresponding transformations. The algebraic properties of sequence transformations satisfying a recursion of the type of Eq. 3.1 had already been studied by Brezinski in his thesis (Brezinski, 1971b, pp. 120 - 127) (see also (Brezinski and Redivo Zaglia, 1991, pp. 106 - 107)).

As shown later, many known variants of epsilon and rho can be classified according to the properties of $F_{k}^{(n)}$. This alone makes our subsequent analysis based on Eq. 3.1 useful. But it is our ultimate goal to construct new generalizations of epsilon and rho.

In principle, we could try to achieve our aim of constructing new transformations by simply choosing new and hopefully effective expressions for $F_{k}^{(n)}$ in Eq. 3.1. In combination with suitable initial conditions, this would produce recursive schemes for new sequence transformations. To the best of our knowledge, questions of that kind have not been treated in the literature yet. So – frankly speaking – we simply do not know how a practically useful $F_{k}^{(n)}$ should look like.

In order to get some help in this matter, we will make some specific assumptions about the dependence of the $F_{k}^{(n)}$ as functions of $k$ and $n$. These assumptions will allow us to investigate whether and under which conditions the $F_{k}^{(n)}$ satisfy certain finite difference equations in the transformation order $k$. It will become clear later that our analysis of the resulting finite difference equations in $k$ leads to some insight which ultimately makes it possible to construct new sequence transformations.

For an analysis of the known variants of epsilon and rho, the full generality of Eq. 3.1 is only needed in the case of the general form (2.10) of Wynn’s rho algorithm with unspecified interpolation points $\{x_{n}\}$. In all other cases, the apparent $n$-dependence of the right-hand side of Eq. 3.1 cancels out and we only have to deal with a simpler, $n$-independent expression of the following general structure:

[TABLE]

In the case of the epsilon algorithm (2.6), we have $F_{k}=1$, in the case of the standard form (2.12) of the rho algorithm, we have $F_{k}=k+1$, and in the case of Osada’s generalization (2.14), we have $F_{k}=k+\theta$.

3.2 First order difference equations

The simplest case $F_{k}=1$ in Eq. 3.2, which corresponds to Wynn’s epsilon algorithm (2.6), satisfies the following first order finite difference equation in $k$:

[TABLE]

Equation 3.3 implies that the recursion (3.2) can also be expressed as follows:

[TABLE]

This reformulation of Eq. 3.2 does not look like a great achievement. It also appears to be a bad idea to use the more complicated recursion (3.2) instead of the much simpler recursion (3.2). However, the more complicated recursion (3.2) will help us to understand better certain features of the simpler recursion (3.2) and of the epsilon algorithm (2.6). In addition, the recursion (3.2) as well as related expressions occurring later open the way for the construction of new sequence transformations.

The first order finite difference equation $\Delta_{m}f(m)=f(m+1)-f(m)=0$ with $m\in\mathbb{N}_{0}$ possesses the general solution $f(m)=f(0)$. Thus, Eq. 3.3 implies $F_{k}=F_{0}$, and Eq. 3.2 can be rewritten as follows:

[TABLE]

If we choose $F_{0}=1$, we obtain Wynn’s epsilon algorithm (2.6).

Alternatively, we could also assume that $F_{0}$ is a non-zero constant. By combining Eq. 3.5 with the epsilon-type initial conditions $T_{-1}^{(n)}=0$ and $T_{0}^{(n)}=S_{n}$, we obtain $T_{2k}^{(n)}=\varepsilon_{2k}^{(n)}$ and $T_{2k+1}^{(n)}=F_{0}\,\varepsilon_{2k+1}^{(n)}$. But only the transforms $T_{2k}^{(n)}=\varepsilon_{2k}^{(n)}$ with even subscripts provide approximations to the limit of the input sequence, whereas the transforms $T_{2k+1}^{(n)}=F_{0}\varepsilon_{2k+1}^{(n)}$ with odd subscripts are only auxiliary quantities. Thus, the modified recursion (3.5) with $F_{0}\neq 0$ does not produce anything new, and the epsilon algorithm remains invariant if we replace in Eq. 2.6b $1/[\varepsilon_{k}^{(n+1)}-\varepsilon_{k}^{(n)}]$ by $F_{0}/[\varepsilon_{k}^{(n+1)}-\varepsilon_{k}^{(n)}]$. To the best of our knowledge, this invariance of the epsilon algorithm had not been observed before.

A similar invariance exists in the case of the general form (2.10) of Wynn’s rho algorithm and its variants. If we replace the interpolation points $\{x_{n}\}$ by the scaled interpolation points $\{\alpha x_{n}\}$ with $\alpha>0$, we obtain $\rho_{2k}^{(n)}(\alpha x_{n})=\rho_{2k}^{(n)}(x_{n})$ and $\rho_{2k+1}^{(n)}(\alpha x_{n})=\alpha\rho_{2k+1}^{(n)}(x_{n})$. This invariance can also be deduced from the determinantal representation for $\rho_{2k}^{(n)}$ (Brezinski, 1978, Theoreme 2.1.52 on p. 97).

Our assumption, that $F_{k}$ satisfies the first order finite difference equation (3.3) in $k$ and is therefore constant as a function of $k$, is apparently too restrictive to allow the construction of useful alternative transformations beyond Wynn’s epsilon algorithm.

Let us now assume that the $n$-dependence of $F_{k}^{(n)}$ is such that the first order finite difference equation

[TABLE]

is satisfied for all $n\in\mathbb{N}_{0}$. Then, the finite difference equation (3.6) implies $F_{k}^{(n)}=F_{0}^{(n)}$ and

[TABLE]

If we choose $F_{0}^{(n)}=x_{n+1}-x_{n}$ in Eq. 3.7, we obtain Brezinski’s first generalization of the epsilon algorithm (Brezinski, 1972, p. 72), and if we choose $F_{0}^{(n)}=x_{n+k+1}-x_{n+k}$, we formally obtain Brezinski’s second generalization (Brezinski, 1972, p. 78). However, Brezinski’s second generalization is only a special case of the recursion (3.7) if the original $k$-dependence of the interpolation points $x_{n+k}$ is eliminated by forming the differences $\Delta x_{n+k}=x_{n+k+1}-x_{n+k}$.

By combining the recursion (3.7) with suitable initial conditions and suitable choices for $F_{0}^{(n)}$, we obtain new sequence transformations. Possible choices of the initial conditions will be be discussed in Section 3.4.

3.3 Second order difference equations

We have $F_{k}=k+1$ in the case of the standard form (2.12) of Wynn’s rho algorithm, and in the case of Osada’s generalization (2.14) we have $F_{k}=k+\theta$. Thus, it looks like a natural idea to assume that $F_{k}$ is linear in $k$, which means that $F_{k}$ is annihilated by $\Delta_{k}^{2}$:

[TABLE]

The second order finite difference equation (3.8) is equivalent to the recursion

[TABLE]

which is even more complicated than the analogous recursion (3.2) based on the first order difference equation (3.3).

The second order finite difference equation $\Delta_{m}^{2}f(m)=f(m+2)-2f(m+1)+f(m)=0$ with $m\in\mathbb{N}_{0}$ possesses the general solution $f(m)=f(0)+m[f(1)-f(0)]$. Thus, Eq. 3.8 implies $F_{k}=F_{0}+k[F_{1}-F_{0}]$, which means that Eq. 3.2 can be reformulated as follows:

[TABLE]

This recursion contains several variants of Wynn’s rho algorithm as special cases. If we choose $F_{0}=1$ and $F_{1}=2$, we obtain the standard form (2.12) of Wynn’s rho algorithm, and if we choose $F_{0}=\theta$ and $F_{1}=\theta+1$, we obtain Osada’s variant (2.14) of the rho algorithm. Finally, $F_{0}=F_{1}\neq 0$ yields epsilon. Other choices for $F_{0}$ and $F_{1}$ are obviously possible, but to the best of our knowledge they have not been explored yet in the literature.

Our approach based on finite difference equations in the transformation order $k$ was successful in the case of the variants of rho mentioned above, but the general form (2.10) of Wynn’s rho algorithm with essentially arbitrary interpolation points $\{x_{n}\}$ remains elusive. In the case of essentially arbitrary interpolation points $\{x_{n}\}$, we cannot express the numerator $x_{n+k+1}-x_{n}$ of the ratio in Eq. 2.10b as an $n$-independent $F_{k}$ according to Eq. 3.2.

At least, we have to employ the more general expression (3.1) involving a general $n$-dependent $F_{k}^{(n)}$. But even this does not guarantee success. Of course, $F_{k}^{(n)}=x_{n+k+1}-x_{n}$ with unspecified and essentially arbitrary interpolation points $\{x_{n}\}$ implies the finite difference equation $\Delta_{k}^{m}F_{k}^{(n)}=\Delta_{k}^{m}x_{n+k+1}$ in $k$. But to achieve any progress, we would have to make specific assumptions on the index dependence of the interpolation points $\{x_{n}\}$.

Therefore, let us now assume that the $n$-dependence of $F_{k}^{(n)}$ is such that

[TABLE]

is satisfied for arbitrary $n\in\mathbb{N}_{0}$. This yields the general solution $F_{k}^{(n)}=F_{0}^{(n)}+k[F_{1}^{(n)}-F_{0}^{(n)}]$, which implies that the recursion (3.10) can be generalized as follows:

[TABLE]

By combining this recursion with suitable initial conditions, we obtain a whole class of new sequence transformations. Possible choices of the initial conditions are discussed in Section 3.4.

3.4 Initial conditions for the new transformation

Let us first analyze the recursion (3.2) which can be reformulated as follows:

[TABLE]

We have to investigate for which values of $k\in\mathbb{N}_{0}$ this recursion can hold, and how its initial conditions could or should be chosen.

It is an obvious idea to choose the initial conditions according to $T_{-1}^{(n)}=0$ and $T_{0}^{(n)}=S_{n}$ as in the case of epsilon. However, these assumptions do not suffice here. Setting $k=-1$ in Eq. 3.13 yields:

[TABLE]

The numerator on the right-hand side contains the undefined quantity $T_{-2}^{(n+1)}$ and the undefined finite difference $\Delta T_{-1}^{(n)}=T_{-1}^{(n+1)}-T_{-1}^{(n)}$. To remove these ambiguities, we could for example assume $\Delta T_{-1}^{(n)}=T_{-1}^{(n+1)}-T_{-1}^{(n)}=0$ and $T_{-2}^{(n+1)}[T_{-1}^{(n+1)}-T_{-1}^{(n)}]=-1$. But in our opinion it would be simpler and more natural to use the recursion (3.13) in combination with suitable initial conditions only for $k\geq 0$.

If we set $k=0$ in Eq. 3.13 and use the epsilon initial conditions $T_{-1}^{(n)}=0$, $T_{0}^{(n)}=S_{n}$, and $T_{1}^{(n)}=1/[T_{0}^{(n+1)}-T_{0}^{(n)}]$, we obtain

[TABLE]

Comparison with Eq. 2.6b shows that $T_{2}^{(n)}$ is nothing but $\varepsilon_{2}^{(n)}$ in disguise. Next, we set $k=1$ in Eq. 3.13:

[TABLE]

At first sight, this looks like a new expression $T_{3}^{(n)}\neq\varepsilon_{3}^{(n)}$. Unfortunately, Eq. 3.15 implies $T_{2}^{(n)}-T_{0}^{(n+1)}=1/[T_{1}^{(n+1)}-T_{1}^{(n)}]$. Thus, we obtain $T_{3}^{(n)}=\varepsilon_{3}^{(n)}$, and with the help of complete induction in $k$ we deduce $T_{k}^{(n)}=\varepsilon_{k}^{(n)}$ for $k\geq 3$, or actually for all $k\in\mathbb{N}_{0}$. These considerations show that the recursion (3.13) produces nothing new if we use it in combination with the epsilon initial conditions $T_{-1}^{(n)}=0$, $T_{0}^{(n)}=S_{n}$, and $T_{1}^{(n)}=1/[T_{0}^{(n+1)}-T_{0}^{(n)}]$.

We will now discuss the choice of alternative initial conditions for the recursion (3.12), which is based on the assumption that $F_{k}^{(n)}$ satisfies the second order finite difference equation (3.11). This is a crucial step, and it is also the point where our approach becomes experimental.

We choose parts of the initial values of the recursion (3.3) as in Wynn’s epsilon algorithm:

[TABLE]

With these initial conditions, epsilon and rho are recovered by choosing $T_{2}^{(n)}=\varepsilon_{2}^{(n)}$ and $T_{2}^{(n)}=\rho_{2}^{(n)}$, respectively.

If we want to employ the recursion (3.3), we first have to make a choice about the initial value $F_{1}^{(n)}$. Equation 3.1 implies $F_{1}^{(n)}=[T_{2}^{(n)}-T_{0}^{(n+1)}][T_{1}^{(n+1)}-T_{1}^{(n)}]$. If $T_{-1}^{(n)}$, $T_{0}^{(n)}$, and $T_{1}^{(n)}$ are chosen according to Eq. 3.17, we only need one further choice for $T_{2}^{(n)}$ to fix $F_{1}^{(n)}$.

It is a relatively natural idea to choose $T_{2}^{(n)}$ in such a way that it corresponds to the special case of a sequence transformation with low transformation order. If we want the resulting new transformation to be more versatile than either epsilon or rho, it makes sense to choose a transformation which is known to accelerate both linear and logarithmic convergence effectively. Obvious candidates would be the spacial case $\vartheta_{2}^{(n)}$ of Brezinski’s theta algorithm discussed in Section 2.3, or suitable variants of Levin’s transformation discussed in Sections B.2 and B.3.

There is another decision, which we have to make. Should the expression for $T_{2}^{(n)}$ depend on $n$ only implicitly via the input data $\{S_{n}\}$, or can it depend on $n$ also explicitly (examples of low order transformations of that kind can be found in (Weniger, 1991, Section 3)). We think that it would be more in agreement with the spirit of the other initial conditions (3.17) if we require that $T_{2}^{(n)}$ depends on $n$ only implicitly via the input data.

Under these assumptions, it looks like an obvious idea to identify $T_{2}^{(n)}$ either with Brezinski’s theta transformation $\vartheta_{2}^{(n)}$ given by Eq. 2.19 or with the explicit expressions for the Levin $u$ and $v$ transformations $u_{2}^{(n)}(\beta,S_{n})$ and $v_{1}^{(n)}(\beta,S_{n})$ given by Eqs. 2.26 and 2.27.

We checked the performance of the resulting algorithms numerically. For $T_{2}^{(n)}=\vartheta_{2}^{(n)}=u_{2}^{(n+1)}(\beta,S_{n+1})=v_{1}^{(n+1)}(\beta,S_{n+1})$ (compare Eqs. 2.26 and 2.27), we obtained consistently good results, but our results for $T_{2}^{(n)}=u_{2}^{(n)}(\beta,S_{n})=v_{1}^{(n)}(\beta,S_{n})$ (compare Eqs. B.8 and B.12) were less good. Thus, we will exclusively use the following non-epsilon initial condition:

[TABLE]

Then, Eq. 3.1 implies

[TABLE]

and we obtain the following new recursive scheme:

[TABLE]

In Eq. 3.20, we wrote $\tilde{T}$ instead of $T$ in order to distinguish our new algorithm from others mentioned before. The performance of our new algorithm (3.20) will be studied in detail in the following Section 4.

4 Numerical examples

In this section, we apply our new algorithm (3.20) to some linearly and logarithmically convergent sequences and also to some divergent series. We compare the numerical results of our new algorithm with those obtained by applying Wynn’s epsilon and rho algorithms (2.6) and (2.12), Osada’s generalized rho algorithm (2.14), Brezinski’s theta algorithm (2.18) and the iteration (2.23) of $\vartheta_{2}^{(n)}$. All the numerical computation are obtained with the help of MATLAB 7.0 in double precision (15 decimal digits).

4.1 Asymptotic model sequences

In this Section, we want to analyze the numerical performance of our new sequence transformation (3.20). As input data, we use some model sequences which provide reasonably good approximations to the elements of numerous practically relevant sequences in the limit of large indices. Essentially the same model sequences had been used by Sidi in (Sidi, 2002, Eqs. (1.1) - (1.3)) and in (Sidi, 2003, Definition 15.3.2 on p. 285), and by Weniger (2004, Section IV):

Logarithmically convergent model sequence

[TABLE] 2. 2.

Linearly convergent model sequence

[TABLE] 3. 3.

Hyperlinearly convergent model sequence

[TABLE]

The series expansions Eqs. 4.1 – 4.3 are to be interpreted as purely formal expansions, i.e., we do not tacitly assume that they converge. It is actually more realistic to assume that they are only asymptotic to $S_{n}$ as $n\to\infty$. There are, however, restrictions on the parameters $\eta$ and $\xi$ if we require that the sequences $\left\{S_{n}\right\}$ defined by the series expansions Eqs. 4.1 – 4.3 are able to represent something finite.

The sequence $\{S_{n}\}$ defined by Eq. 4.1 can only converge to its limit $S$ if $\eta<0$ holds. For $\eta>0$, Sidi (2002, Text following Eq. (1.1)) calls $S$ the antilimit of this sequence. While this terminology is very useful in the case of divergent, but summable series, it is highly questionable here. The divergence of this sequence for $\eta>0$ is genuine, and the normally used summation techniques are not able to produce anything finite. Thus, the terminology antilimit does not make sense here.

The sequence $\left\{S_{n}\right\}$ defined by Eq. 4.2 converges for all $\eta\in\mathbb{R}$ linearly to its limit $S$ provided that $|\xi|<1$ holds. For $|\xi|>1$, this sequence diverges, but it should normally be summable to $S$ for all $\xi$ belonging to the cut complex plane $\mathbb{C}\setminus[1,\infty)$.

The most convenient situation occurs if the sequence elements $S_{n}$ satisfy Eq. 4.3. Then, this sequence converges for all $\xi,\eta\in\mathbb{R}$ hyperlinearly to its limit $S$. In Eq. 4.3, the restriction $r\in\mathbb{R}_{+}$ is essential, because for $r<0$, we would obtain a sequence which diverges (hyper)factorially for all $\xi\neq 0$. Divergent sequences of that kind can according to experience often be summed to their generalized limits $S$ for all $\xi\in\mathbb{C}\setminus[0,\infty)$, but a rigorous theoretical analysis of such a summation process is normally extremely difficult (compare the rigorous theoretical analysis of the summation of a comparatively simple special case in (Borghi and Weniger, 2015) and references therein).

4.2 Linearly convergent sequences

In this Subsection, we apply our new algorithm to some linearly convergent alternating and monotone series, which are special cases of the model sequence (4.1). The numerical results obtained by our new algorithm (3.20), Wynn’s epsilon algorithm (2.6), Brezinski’s theta algorithm (2.18) and the iteration (2.23) of $\vartheta_{2}^{(n)}$ are presented in Tables 1 and 2.

In the captions of the following Tables, $\lfloor x\rfloor$ denotes the integral part of $x$, which is the largest integer $n$ satisfying $n\leq x$.

Example 4.1: We consider the alternating partial sums

[TABLE]

which converge to $S=\ln(2)$. The Boole summation formula (Olver, Lozier, Boisvert, and Clark, 2010, Eq. (24.17.1)) yields

[TABLE]

This truncated asymptotic expansion is a special case of the model sequence (4.2) with $\xi=-1$ and $\eta=-1$.

Example 4.2: The linearly convergent monotone sequence

[TABLE]

converges to $S=\ln(5)$ as $n\to\infty$. As shown by Sidi (2003, p. 84), $S_{n}$ possesses the asymptotic expansion

[TABLE]

which is a special case of the model sequence (4.2) with $\xi=4/5$ and $\eta=-1$.

4.3 Logarithmically convergent sequences

First, we consider a logarithmically convergent sequence of the type of the model sequence (2.13) with integral decay parameter $\theta$. In Table 3, we present the numerical results obtained by applying our new algorithm (3.20), the standard form (2.12) of Wynn’s rho algorithm, Brezinski’s theta algorithm (2.18), and the iteration (2.23) of $\vartheta_{2}^{(n)}$. Other logarithmically convergent sequences with non-integral decay parameters $\theta$ will be discussed later.

Example 4.3: We consider the logarithmically convergent sequence

[TABLE]

which converges to $S=0.57721~{}56649~{}01532\cdots$ (Euler’s constant).

The Euler-Maclaurin formula (see for example (Olver, Lozier, Boisvert, and Clark, 2010, Eq. (2.10.1))) yields:

[TABLE]

where $B_{n}$ is a Bernoulli number (Olver, Lozier, Boisvert, and Clark, 2010, Eq. (24.2.1)).

Next, we consider logarithmically convergent sequences of the form of Eq. 4.1 with non-integral decay $\theta$. In Table 4, we present the numerical results for the test problem (4.10) obtained by our new algorithm (3.20), Osada’s generalized rho algorithm (2.14) with $\theta=1/2$, Brezinski’s theta algorithm (2.18) and the iteration (2.23) of $\vartheta_{2}^{(n)}$.

Example 4.4: A very useful test problem is the series expansion for the so-called lemniscate constant $A$ (Todd, 1975, Theorem 5):

[TABLE]

Standard asymptotic techniques show that the terms $(2m-1)!!/[(2m)!!(4m+1)]$ of this series decay like $m^{-3/2}$ as $m\to\infty$. This implies that the remainders

[TABLE]

decay like $n^{-1/2}$ as $n\to\infty$. Consequently, the series expansion (4.10) constitutes an extremely challenging convergence acceleration problem. Todd (1975, p. 14) stated in 1975 that this series expansion is of no practical use for the computation of $A$, which is certainly true if no convergence acceleration techniques are available. In 1979, the series expansion (4.10) was used by Smith and Ford (1979, Entry 13 of Table 6.1) to test the ability of sequence transformations of accelerating logarithmic convergence. They observed that the standard form (2.12) of Wynn’s rho algorithm fails to accelerate the convergence of this series (Smith and Ford, 1979, p. 235) which is in agreements with Osada’s later convergence analysis (Osada, 1990, Theorem 3.2). The series expansion (4.10) was also extensively used by Weniger (1989, Tables 14-2 - 14-4) as a test case.

Example 4.5: The series expansion (A.10) of $1/z$ in terms of reduced Bessel functions, whose basic properties are reviewed in Appendix A, is a challenging test problem for the ability of a sequence transformation to accelerate the convergence of a logarithmically convergent sequence with a non-integral decay parameter.

We consider the sequence of partial sum of the series (A.10) with $z=1$,

[TABLE]

which corresponds to a special case of the model sequence (4.1) with $\eta=-1/2$. In Table 5, we present the numerical results for the test problem (4.12) obtained by our new algorithm (3.20), Osada’s generalized rho algorithm (2.14) with $\theta=1/2$, Brezinski’s theta algorithm (2.18) and the iteration (2.23) of $\vartheta_{2}^{(n)}$.

Our numerical Examples 4.3 - 4.5 show that our new algorithm is apparently able to accelerate the convergence of logarithmically convergent sequences of the type of (2.13) without having to know the value of the decay parameter $\theta$ explicitly, as it is necessary in the case of Osada’s variant (2.14) of the rho algorithm.

4.4 Divergent series

In this Section, we apply our new algorithm for the summation of the factorially divergent Euler series. Since the early days of calculus, divergent series have been a highly controversial topic in mathematics (Burkhardt, 1911; Ferraro, 2008; Tucciarone, 1973), and to some extend they still are (see for example (Weniger, 2008, Appendix D) and references therein).

The exponential integral (Olver, Lozier, Boisvert, and Clark, 2010, Eq. (6.2.1))

[TABLE]

possesses the asymptotic expansion (Olver, Lozier, Boisvert, and Clark, 2010, Eq. (6.12.1))

[TABLE]

The exponential integral can be expressed as a Stieltjes integral (Olver, Lozier, Boisvert, and Clark, 2010, Eq. (6.7.1 )),

[TABLE]

which implies that the divergent inverse power series (4.14) is a Stieltjes series.

Asymptotic series as $z\to\infty$ involving a factorially divergent generalized hypergeometric series ${}_{2}F_{0}$ occur also among other special functions. Examples are the asymptotic expansion of the modified Bessel function $K_{\nu}(z)$ of the second kind (Olver, Lozier, Boisvert, and Clark, 2010, Eq. (10.40.2)), or the asymptotic expansion of the Whittaker function $W_{\kappa,\mu}(z)$ of the second kind (Olver, Lozier, Boisvert, and Clark, 2010, Eq. (13.19.3)). Similar divergent series occur also in statistics, as discussed in a book by Bowman and Shenton (1989).

If we replace in Eq. 4.14 $z$ by $1/z$, we obtain the so-called Euler series, which had already been studied by Euler (see for example the books by Bromwich (1991, pp. 323 - 324) and Hardy (1949, pp. 26 - 29) or the articles by Barbeau (1979) and Barbeau and Leah (1976)). For $|\arg(z)|<\pi$, the Euler series is the asymptotic series of the Euler integral:

[TABLE]

The Euler series is the most simple prototype of the factorially divergent perturbation series that occur abundantly in quantum physics. Already Dyson (1952) had argued that perturbation expansions in quantum electrodynamics must diverge factorially. Around 1970, Bender and Wu (1969, 1971, 1973) showed in their seminal work on anharmonic oscillators that factorially divergent perturbation expansions occur in nonrelativistic quantum mechanics. Later, it was found that factorially divergent perturbation expansions are actually the rule in quantum physics rather than the exception (see for example the articles by Fischer (1997, Table 1) and Suslov (2005), or the articles reprinted in the book by Le Guillou and Zinn-Justin (1990)). Applications of factorially divergent series in optics are discussed in a recent review by Borghi (2016).

The Euler integral $\mathcal{E}(z)$ is a Stieltjes function. Thus, the corresponding Euler series ${}_{2}F_{0}(1,1;-z)$ is its associated Stieltjes series (see for example (Baker and Graves-Morris, 1996, Chapter 5)). This Stieltjes property guarantees the convergence of certain subsequences of the Padé table of ${}_{2}F_{0}(1,1;-z)$ to the Euler integral (see for example (Baker and Graves-Morris, 1996, Chapter 5.2) or the discussion in (Borghi and Weniger, 2015)).

It has been shown in many calculations that sequence transformations are principal tools that can accomplish an efficient summation of the factorially divergent expansions of the type of the Euler series (see for example (Borghi and Weniger, 2015) and references therein). However, in the case of most sequence transformations, no rigorous theoretical convergence proofs are known (this applies also to Padé approximants if the series to be transformed is not Stieltjes). An exception is the delta transformation (Weniger, 1989, Eq. (8.4-4)) of the Euler series, whose convergence to the Euler integral for all $z\in\mathbb{C}\setminus[0,\infty)$ was proved rigorously in (Borghi and Weniger, 2015).

In the following Tables 6 and 7, we use as input data the partial sums

[TABLE]

of the factorially divergent asymptotic series (4.14) for the exponential integral $E_{1}(z)$. We present the summation results obtained for $z=3$ and $z=1/2$, respectively, by our new algorithm (3.20), Wynn’s epsilon algorithm (2.6), Brezinski’s theta algorithm (2.18), and the iteration (2.23) of $\vartheta_{2}^{(n)}$.

5 Conclusions and discussions

Starting from the known recursion relations of Wynn’s epsilon and rho algorithm or of Osada’s generalized rho algorithm, we constructed a new convergence acceleration algorithm (3.20) which is – as shown in our numerical examples in Section 4 – not only successful in the case of linearly convergent sequences, but works also in the case of many logarithmically convergent sequences. Our numerical results also showed that our new transformations sums factorially divergent series. Thus, our new transformation (3.20) is clearly more versatile than either epsilon or rho, having similar properties as Brezinski’s theta algorithm (2.18) or the iteration (2.23) of $\vartheta_{2}^{(n)}$.

For the derivation of our new algorithm (3.20) in Section 3, we pursued a novel approach. We analyzed the finite difference equations with respect to the transformation order, which the known recurrence formulas of Wynn’s epsilon and rho algorithms or of Osada’s generalized rho algorithm satisfy. This approach yielded a generalized recursive scheme (3.12) which contains the known recursions of Wynn’s epsilon and rho algorithms and Osada’s rho algorithm as special cases.

To complete the derivation of our new algorithm (3.20), we only had to choose appropriate initial conditions. Parts of the initial conditions were chosen as in Wynn’s epsilon algorithm according to Eq. 3.17. In addition to these epsilon-type initial conditions (3.17), we then needed only one further initial condition for $T_{2}^{(n)}$, which we identified with Brezinski’s theta transformation $\vartheta_{2}^{(n)}$ according to Section 3.4. Thus, our new algorithm (3.20) can be viewed to be essentially a blend of Wynn’s epsilon and Brezinski’s theta algorithm.

Of course, it should be possible to construct – starting from the generalized recursive scheme (3.12) – alternative new sequence transformations by modifying our choice for $T_{2}^{(n)}$. The feasibility of this approach will be investigated elsewhere.

Acknowledgements

Y. He was supported by the National Natural Science Foundation of China (Grants no. 11571358), the China Postdoctoral Science Foundation funded project (Grants no. 2012M510186 and 2013T60761), and the Youth Innovation Promotion Association CAS. X.K. Chang was supported by the National Natural Science Foundation of China (Grants no. 11701550, 11731014), and the Youth Innovation Promotion Association CAS. X.B. Hu was supported by the National Natural Science Foundation of China (Grants no. 11331008, 11571358). J.Q. Sun was supported by the National Natural Science Foundation of China (Grants no. 11401546, 11871444), and the Shandong Provincial Natural Science Foundation, China (Grant no. ZR2013AQ024).

\appendixpage\addappheadtotoc

Appendix A Reduced Bessel functions

In this Appendix, we briefly discuss the most important properties of the so-called reduced Bessel functions (Steinborn and Filter, 1975, Eqs. (3.1) and (3.2)) and their anisotropic generalizations, the so-called $B$ functions (Filter and Steinborn, 1978a, Eqs. (3.3) and (3.4)). These functions have gained some prominence in molecular electronic structure theory (see for example Weniger (2009) and references therein), but have been largely ignored in mathematics. Reduced Bessel and $B$ functions had also been the topic of Weniger’s Diploma and PhD theses Weniger (1977, 1982). But in this article, we are predominantly interested in reduced Bessel functions because of the series expansion (A.10) of $1/z$ in terms of reduced Bessel functions, which is well suited to test the ability of a sequence transformation to accelerate the logarithmic convergence of model sequences of the type of Eq. 2.13 with a non-integral decay parameter.

Based on previous work by Shavitt (1963, Eq. (55) on p. 15), the reduced Bessel function with in general complex order $\nu$ and complex argument $z$ was introduced by Steinborn and Filter (1975, Eqs. (3.1) and (3.2)):

[TABLE]

Here, $K_{\nu}(z)$ is a modified Bessel functions of the second kind (Olver, Lozier, Boisvert, and Clark, 2010, Eqs. (10.27.4) and (10.27.5)).

Many properties of the reduced Bessel functions follow directly from the corresponding properties of $K_{\nu}(z)$. One example is the recursion

[TABLE]

which is stable upwards and which follows directly from the three-term recursion $K_{\nu+1}(z)=(2\nu/z)K_{\nu}(z)+K_{\nu-1}(z)$ (Olver, Lozier, Boisvert, and Clark, 2010, Eq. (10.29.1)). Another example is the symmetry relationship $\widehat{k}_{-\nu}(z)=z^{-2\nu}\widehat{k}_{\nu}(z)$, which follows from the symmetry relationship $K_{-\nu}(z)=K_{\nu}(z)$ (Olver, Lozier, Boisvert, and Clark, 2010, Eq. (10.27.3)).

If the order $\nu$ of $\widehat{k}_{\nu}(z)$ is half-integral, $\nu=n+1/2$ with $n\in\mathbb{N}_{0}$, a reduced Bessel function can be expressed as an exponential multiplied by a terminating confluent hypergeometric series ${}_{1}F_{1}$ (see for example (Weniger and Steinborn, 1983, Eq. (3.7))):

[TABLE]

The functions $\widehat{k}_{n+1/2}(z)$ can be computed conveniently and reliably with the help of the recursion (A.2) in combination with the starting values $\widehat{k}_{-1/2}(z)=\mathrm{e}^{-z}/z$ and $\widehat{k}_{1/2}(z)=\mathrm{e}^{-z}$.

The polynomial part of $\widehat{k}_{n+1/2}(z)$ had been considered independently in the mathematical literature. There, the notation

[TABLE]

is used (Grosswald, 1978, Eq. (1) on p. 34). Together with some other, closely related polynomials, the $\theta_{n}(z)$ are called Bessel polynomials Grosswald (1978). According to Grosswald (1978, Section 14), they have been applied in such diverse and seemingly unrelated fields like number theory, statistics, and the analysis of complex electrical networks.

Bessel polynomials occur also in the theory of Padé approximants. Padé (1892, p. 82) had shown in his thesis that the Padé approximant $[n/m]_{\exp}(z)$ to the exponential function $\exp(z)$ can be expressed in closed form, and Baker and Graves-Morris showed that Padé’s expression corresponds in modern mathematical notation to the following ratio of two terminating confluent hypergeometric series ${}_{1}F_{1}$ Baker and Graves-Morris (1996, Eq. (2.12)):

[TABLE]

Thus, the diagonal Padé approximants $[n/n]_{\exp}(z)$ are ratios of two Bessel polynomials (Weniger, 1989, Eq. (14.3-15)):

[TABLE]

The known monotonicity properties of the modified Bessel function $K_{\nu}(z)$ imply that the reduced Bessel functions ${\widehat{k}}_{\nu}(z)$ are for $\nu>0$ and $z\geq 0$ positive and bounded by their values at the origin (Weniger and Steinborn, 1983, Eq. (3.1)). In the case of reduced Bessel functions with half-integral orders, this yields the following bound:

[TABLE]

In Grosswald’s book (Grosswald, 1978), it is shown that for fixed and finite argument $z$ the Bessel polynomials $\theta_{n}(z)$ satisfy the leading order asymptotic estimate (Grosswald, 1978, p. 125)

[TABLE]

Combination of Eqs. A.4 and A.8 shows that the dominant term of the Poincaré-type asymptotic expansion of $\widehat{k}_{n+1/2}(z)$ with fixed and finite argument $z$ corresponds to its value at the origin (Weniger and Steinborn, 1983, Eq. (3.9)):

[TABLE]

For several functions, finite or infinite expansions in terms of reduced Bessel functions are known. As a challenging problem for sequence transformations, we propose to use the series expansion (Filter and Steinborn, 1978b, Eq. (6.5))

[TABLE]

Equation A.9 implies that the terms $\widehat{k}_{m-1/2}(z)/[2^{m}m!]$ of this series decay like $m^{-3/2}$ as $m\to\infty$ (Grotendorst et al, 1986, p. 3709), or that the remainders $\sum_{m=n+1}^{\infty}\,\widehat{k}_{m-1/2}(z)/[2^{m}m!]$ decay like $n^{-1/2}$ as $n\to\infty$. Thus, the series (A.10) converges as slowly as the Dirichlet series $\zeta(s)=\sum_{m=0}^{\infty}(m+1)^{s}$ with $s=1/2$, which is notorious for extremely poor convergence. The slow convergence of the infinite series (A.10) was demonstrated in (Grotendorst et al, 1986, Table I). After $10^{6}$ terms, only 3 decimal digits were correct, which is in agreement with the truncation error estimate given above.

Appendix B Levin’s transformation

B.1 The general Levin transformation

A sequence transformations of the type of Wynn’s epsilon algorithm (2.6) only requires the input of a finite sub-string of a sequence $\{S_{n}\}_{n=0}^{\infty}$. No additional information is needed for the computation of approximations to the (generalized) limit $S$ of an input sequence $\{S_{n}\}_{n=0}^{\infty}$. Obviously, this is a highly advantageous feature. However, in more fortunate cases, additional information on the index dependence of the truncation errors $R_{n}=S_{n}-S$ is available. For example, the truncation errors of Stieltjes series are bounded in magnitude by the first terms neglected in the partial sums (see for example (Weniger, 1989, Theorem 13-2)), and they also have the same sign pattern as the first terms neglected. The utilization of such a structural information in a transformation process should enhance its efficiency. Unfortunately, there is no obvious way of incorporating such an information into Wynn’s recursive epsilon algorithm (2.6) or into other sequence transformations with similar features.

In 1973, Levin (1973) introduced a sequence transformation which overcame these limitations and which now bears his name. It uses as input data not only the elements of the sequence $\{S_{n}\}_{n=0}^{\infty}$, which is to be transformed, but also the elements of another sequence $\{\omega_{n}\}_{n=0}^{\infty}$ of explicit estimates of the remainders $R_{n}=S_{n}-S$. These remainders estimates, which must be explicitly known, make it possible to incorporate additional information into the transformation process and are thus ultimately responsible for the remarkable versatility of Levin’s transformation.

For Levin’s transformation, we use the notation introduced in (Weniger, 1989, Eqs. (7.1-6) and (7.1-7)):

[TABLE]

Here, $\beta>0$ is a shift parameter, $\{S_{n}\}_{n=0}^{\infty}$ is the input sequence, and $\{\omega_{n}\}_{n=0}^{\infty}$ is a sequence of remainder estimates.

In Eqs. B.1 and B.2 it is essential that the input sequence $\{S_{n}\}$ starts with the sequence element $S_{0}$. The reason is that both the definitions according to Eqs. B.1 and B.2 as well as the recurrence formulas for its numerator and denominator sums (see for example(Weniger, 2004, Eq. (3.11))) depend explicitly on $n$ as well as on $\beta$.

Levin’s transformation, which is also discussed in the NIST Handbook (Olver, Lozier, Boisvert, and Clark, 2010, §3.9(v) Levin’s and Weniger’s Transformations), is generally considered to be both a very powerful and a very versatile sequence transformation (see for example (Borghi and Weniger, 2015; Brezinski and Redivo Zaglia, 1991; Sidi, 2003; Smith and Ford, 1979, 1982; Weniger, 1989, 2004) and references therein). The undeniable success of Levin’s transformation inspired others to construct alternative sequence transformations that also use explicit remainder estimates (see for example (Čížek, Zamastil, and Skála, 2003; Homeier, 2000a; Weniger, 1989, 1992, 2004) and references therein).

We still have to discuss the choice of the so far unspecified remainder estimates $\{\omega_{n}\}_{n=0}^{\infty}$. A principal approach would be to look for remainder estimates that reproduce the leading order asymptotics of the actual remainders (Weniger, 1989, Eq. (7.3-1)):

[TABLE]

This is a completely valid approach, but there is the undeniable problem that for every input sequence $\{S_{n}\}_{n=0}^{\infty}$ the leading order asymptotics of the corresponding remainders $R_{n}=S_{n}-S$ as $n\to\infty$ has to be determined. Unfortunately, such an asymptotic analysis may well lead to some difficult technical problem and be a non-trivial research problem in its own right.

In practice, it is much more convenient to use simple explicit remainder estimates introduced by Levin (1973) and Smith and Ford (1979), respectively, which are known to work well even in the case of purely numerical input data. In this article, we only consider Levin’s remainder estimates $\omega_{n}=(\beta+n)\Delta S_{n-1}$, and $\omega_{n}=-[\Delta S_{n-1}][\Delta S_{n}]/[\Delta^{2}S_{n-1}]$, which lead to Levin’s $u$ and $v$ variants (Levin, 1973, Eqs. (58) and (67)). These estimates lead to transformations that are known to be powerful accelerators for both linear and logarithmic convergence (compare also the asymptotic estimates in (Weniger, 2004, Section IV)). Therefore, Levin’s $u$ and $v$ transformations are principally suited to represent the remaining initial condition $T_{2}^{(n)}$ in Eq. 3.19.

B.2 Levin’s $u$ transformation

Levin’s $u$ transformation, which first appeared already in 1936 in the work of Bickley and Miller (1936), is characterized by the remainder estimate $\omega_{n}=(\beta+n)\Delta S_{n-1}$ (Levin, 1973, Eq. (58)). Inserting this into the finite difference representation (B.1) yields:

[TABLE]

For $k=2$, the right-hand side of this expression depends on $n$ only implicitly via the input data $\{S_{n}\}$ and also does not depend on the scaling parameter $\beta$. Thus, we obtain:

[TABLE]

If we now use $S_{n}=S_{n-1}+\Delta S_{n-1}$, we obtain $S_{n}/\Delta S_{n-1}=S_{n-1}/\Delta S_{n-1}+1$, which implies $\Delta[S_{n}/\Delta S_{n-1}]=\Delta[S_{n-1}/\Delta S_{n-1}]$. Thus, we obtain the alternative expression,

[TABLE]

which can be reformulated as follows:

[TABLE]

In order to enhance numerical stability, it is recommendable to rewrite this expression as follows:

[TABLE]

B.3 Levin’s $v$ transformation

Levin’s $v$ transformation is characterized by the remainder estimate (Levin, 1973, Eq. (67)) which is inspired by Aitken’s $\Delta^{2}$ formula (1.8):

[TABLE]

In Eq. B.9, we made use of the fact that Levin’s transformation $\mathcal{L}_{k}^{(n)}(\beta,S_{n},\omega_{n})$ is a homogeneous function of degree zero of the $k+1$ remainder estimates $\omega_{n}$, $\omega_{n+1}$, …. Thus, we can use any of the expressions in Eq. B.9 without affecting the value of the transformation.

Inserting the remainder estimate (B.9) into the finite difference representation (B.1) yields:

[TABLE]

For $k=1$, the right-hand side of this expression depends on $n$ only implicitly via the input data $\{S_{n}\}$ and also does not depend on the scaling parameter $\beta$. We obtain:

[TABLE]

Comparison of Eqs. B.7 and B.11 yields $v_{1}^{(n)}(\beta,S_{n})=u_{2}^{(n)}(\beta,S_{n})$ for arbitrary $\beta\in\mathbb{R}$ and $n\in\mathbb{N}_{0}$.

Again, it is recommendable to rewrite Eq. B.11 as follows:

[TABLE]