Construction of new generalizations of Wynn’s epsilon and rho algorithm by solving finite difference equations in the transformation order

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.

Appendices

Appendix 1. Reduced Bessel functions

In this Appendix, we briefly discuss the most important properties of the so-called reduced Bessel functions [76, Eqs. (3.1) and (3.2)] and their anisotropic generalizations, the so-called B functions [42, Eqs. (3.3) and (3.4)]. These functions have gained some prominence in molecular electronic structure theory (see, for example, [96] 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 [86, 87]. But in this article, we are predominantly interested in reduced Bessel functions because of the series expansion (A1.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 (2.13) with a non-integral decay parameter.

Based on previous work by Shavitt [71, Eq. (55) on p. 15], the reduced Bessel function with in general complex order ν and complex argument z was introduced by Steinborn and Filter [76, Eqs. (3.1) and (3.2)]:

$$ \widehat{k}_{\nu} (z) = (2/\pi)^{1/2} z^{\nu} K_{\nu} (z) . $$

(A1.1)

Here, K_ν(z) is a modified Bessel functions of the second kind [59, Eqs. (10.27.4) and (10.27.5)].

Many properties of the reduced Bessel functions follow directly from the corresponding properties of K_ν(z). One example is the recursion

$$ \widehat{k}_{\nu+ 1} (z) = 2 \nu \widehat{k}_{\nu} (z) + z^{2} \widehat{k}_{\nu-1} (z) , $$

(A1.2)

which is stable upwards and which follows directly from the three-term recursion K_ν+ 1(z) = (2ν/z)K_ν(z) + K_ν− 1(z) [59, 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_−ν(z) = K_ν(z) [59, Eq. (10.27.3)].

If the order ν of $\widehat {k}_{\nu } (z)$ is half-integral, ν = 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 ₁F₁ (see, for example, [99, Eq. (3.7)]):

$$ \widehat{k}_{n + 1/2} (z) = 2^{n} (1/2)_{n} \mathrm{e}^{-z} {~}_{1} F_{1} (-n; -2n; 2z) , \qquad n \in \mathbb{N}_{0} . $$

(A1.3)

The functions $\widehat {k}_{n + 1/2} (z)$ can be computed conveniently and reliably with the help of the recursion (A1.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

$$ \theta_{n} (z) = \mathrm{e}^{z} \widehat{k}_{n + 1/2} (z) = 2^{n} (1/2)_{n} \mathrm{e}^{-z} {~}_{1} F_{1} (-n; -2n; 2z) $$

(A1.4)

is used [48, Eq. (1) on p. 34]. Together with some other, closely related polynomials, the θ_n(z) are called Bessel polynomials [48]. According to Grosswald [48, 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é [63, 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 ₁F₁ Baker and Graves-Morris [2, Eq. (2.12)]:

$$ [n/m]_{\exp} (z) = \frac {{~}_{1} F_{1} (- n; - n - m; z)} {{~}_{1} F_{1} (- m; - n - m; - z)} , \qquad m, n \in \mathbb{N}_{0} . $$

(A1.5)

Thus, the diagonal Padé approximants [n/n]_exp(z) are ratios of two Bessel polynomials [88, Eq. (14.3-15)]:

$$ [n/n]_{\exp} (z) = \frac{\theta_{n} (z/2)}{\theta_{n} (-z/2)} , \qquad n \in \mathbb{N}_{0} . $$

(A1.6)

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

$$ 0 < \widehat{k}_{n + 1/2} (z) \le \widehat{k}_{n + 1/2} (0) = 2^{n} (1/2)_{n} , \qquad 0 \le z < \infty . \quad n \in \mathbb{N}_{0} . $$

(A1.7)

In Grosswald’s book [48], it is shown that for fixed and finite argument z the Bessel polynomials θ_n(z) satisfy the leading-order asymptotic estimate [48, p. 125]

$$ \theta_{n} (z) \sim \frac {(2 n)!} {2^{n} n!} \mathrm{e}^{z} = 2^{n} (1/2)_{n} \mathrm{e}^{z} , \qquad n \to \infty . $$

(A1.8)

Combination of (A1.4) and (A1.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 [99, Eq. (3.9)]:

$$ \widehat{k}_{n + 1/2} (z) = \widehat{k}_{n + 1/2} (0) \Bigl[1 + \mathrm{O} \bigl(n^{- 1} \bigr) \Bigr] = 2^{n} (1/2)_{n} \Bigl[1 + \mathrm{O} \bigl(n^{- 1} \bigr) \Bigr] , \qquad n \to \infty . $$

(A1.9)

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 [43, Eq. (6.5)]

$$ \frac{1}{z} = \sum\limits_{m = 0}^{\infty} \widehat{k}_{m-1/2} (z) / [2^{m} m!] , \qquad z > 0 . $$

(A1.10)

Equation (A1.9) implies that the terms $\widehat {k}_{m-1/2} (z) / [2^{m} m!]$ of this series decay like m^− 3/2 as m →∞ [49, 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 →∞. Thus, the series (A1.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 (A1.10) was demonstrated in [49, Table I]. After 10⁶ terms, only 3 decimal digits were correct, which is in agreement with the truncation error estimate given above.

Appendix 2. Levin’s transformation

2.1 A2.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, [88, 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 information into Wynn’s recursive epsilon algorithm (2.6) or into other sequence transformations with similar features.

In 1973, Levin [55] 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 remainder 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 [88, Eqs. (7.1-6) and (7.1-7)]:

$$ \begin{array}{@{}rcl@{}} \mathcal{L}_{k}^{(n)} (\beta, S_{n}, \omega_{n}) & =& \frac {{\Delta}^{k} [(\beta+n)^{k-1} S_{n} / \omega_{n}]} {{\Delta}^{k} [(\beta+n)^{k-1} / \omega_{n}]} \end{array} $$

(A2.1)

$$ \begin{array}{@{}rcl@{}} & = & \frac {\displaystyle \sum\limits_{j = 0}^{k} (-1)^{j} {\binom{k}{j}} \frac {(\beta+n+j)^{k-1}} {(\beta+n+k)^{k-1}} \frac {S_{n+j}} {\omega_{n+j}}} {\displaystyle \sum\limits_{j = 0}^{k} (-1)^{j} {\binom{k}{j}} \frac {(\beta+n+j)^{k-1}} {(\beta+n+k)^{k-1}} \frac {1} {\omega_{n+j}} } , \quad \!k, n \in \mathbb{N}_{0} , \quad \!\beta \!>\! 0 .\\ \end{array} $$

(A2.2)

Here, β > 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 (A2.1) and (A2.2) it is essential that the input sequence {S_n} starts with the sequence element S₀. The reason is that both the definitions according to (A2.1) and (A2.2) as well as the recurrence formulas for its numerator and denominator sums (see, for example, [93, Eq. (3.11)]) depend explicitly on n as well as on β.

Levin’s transformation, which is also discussed in the NIST Handbook [59, §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, [13, 28, 73,74,75, 88, 93] 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, [36, 52, 88, 90, 93] 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 [88, Eq. (7.3-1)]:

$$ R_{n} = S_{n} - S = \omega_{n} \left[ C + \mathrm{O} \bigl(1/n \bigr) \right] , \qquad C \ne 0 , \quad n \to \infty . $$

(A2.3)

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 →∞ 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 [55] and Smith and Ford [74], 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 ω_n = (β + n)ΔS_n− 1, and ω_n = −[ΔS_n− 1][ΔS_n]/[Δ²S_n− 1], which lead to Levin’s u and v variants [55, 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 [93, Section IV]). Therefore, Levin’s u and v transformations are principally suited to represent the remaining initial condition $T_{2}^{(n)}$ in (3.19).

2.2 A2.2. Levin’s u transformation

Levin’s u transformation, which first appeared already in 1936 in the work of Bickley and Miller [10], is characterized by the remainder estimate ω_n = (β + n)ΔS_n− 1 [55, Eq. (58)]. Inserting this into the finite difference representation (A2.1) yields:

$$ u_{k}^{(n)} (\beta, S_{n}) = \frac {{\Delta}^{k} [(\beta+n)^{k-2} S_{n} / {\Delta} S_{n-1}]} {{\Delta}^{k} [(\beta+n)^{k-2} / {\Delta} S_{n-1}]} . $$

(A2.4)

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 β. Thus, we obtain:

$$ u_{2}^{(n)} (\beta, S_{n}) = \frac {{\Delta}^{2} \bigl[ S_{n} / {\Delta} S_{n-1} \bigr]} {{\Delta}^{2} \bigl[ 1 / {\Delta} S_{n-1} \bigr]} . $$

(A2.5)

If we now use S_n = S_n− 1 + ΔS_n− 1, we obtain S_n/ΔS_n− 1 = S_n− 1/ΔS_n− 1 + 1, which implies Δ[S_n/ΔS_n− 1] = Δ[S_n− 1/ΔS_n− 1]. Thus, we obtain the alternative expression,

$$ u_{2}^{(n)} (\beta, S_{n}) = \frac {{\Delta}^{2} [ S_{n-1} / {\Delta} S_{n-1} ]} {{\Delta}^{2} [ 1 / {\Delta} S_{n-1} ]} , $$

(A2.6)

which can be reformulated as follows:

$$ u_{2}^{(n)} (\beta, S_{n}) = \frac {S_{n} \bigl[ {\Delta} S_{n + 1} \bigr] \bigl[ {\Delta}^{2} S_{n-1} \bigr] - S_{n + 1} \bigl[ {\Delta} S_{n-1} \bigr] \bigl[ {\Delta}^{2} S_{n} \bigr]} {\bigl[ {\Delta} S_{n + 1} \bigr] \bigl[ {\Delta}^{2} S_{n-1} \bigr] - \bigl[ {\Delta} S_{n-1} \bigr] \bigl[ {\Delta}^{2} S_{n} \bigr]} . $$

(A2.7)

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

$$ u_{2}^{(n)} (\beta, S_{n}) = S_{n} - \frac {\bigl[ {\Delta} S_{n-1} \bigr] \bigl[ {\Delta} S_{n} \bigr] \bigl[ {\Delta}^{2} S_{n} \bigr]} {\bigl[ {\Delta} S_{n + 1} \bigr] \bigl[ {\Delta}^{2} S_{n-1} \bigr] - \bigl[ {\Delta} S_{n-1} \bigr] \bigl[ {\Delta}^{2} S_{n} \bigr]} . $$

(A2.8)

2.3 A2.3. Levin’s v transformation

Levin’s v transformation is characterized by the remainder estimate [55, Eq. (67)] which is inspired by Aitken’s Δ² formula (1.8):

$$ \omega_{n} = \frac{[{\Delta} S_{n-1}] [{\Delta} S_{n}]} {{\Delta} S_{n-1} - {\Delta} S_{n}} = - \frac{[{\Delta} S_{n-1}] [{\Delta} S_{n}]}{{\Delta}^{2} S_{n-1}} = \frac{[{\Delta} S_{n-1}] [{\Delta} S_{n}]}{{\Delta}^{2} S_{n-1}} , $$

(A2.9)

In (A2.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 ω_n, ω_n+ 1, …. Thus, we can use any of the expressions in (A2.9) without affecting the value of the transformation.

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

$$ v_{k}^{(n)} (\beta, S_{n}) = \frac {\displaystyle {\Delta}^{k} \frac{\displaystyle (\beta+n)^{k-1} S_{n} [{\Delta}^{2} S_{n-1}]}{\displaystyle [{\Delta} S_{n-1}] [{\Delta} S_{n}]}} {\displaystyle {\Delta}^{k} \frac{\displaystyle (\beta+n)^{k-1} [{\Delta}^{2} S_{n-1}]}{\displaystyle [{\Delta} S_{n-1}] [{\Delta} S_{n}]}} . $$

(A2.10)

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 β. We obtain:

$$ v_{1}^{(n)} (\beta, S_{n}) = \frac {\displaystyle S_{n + 1} [{\Delta} S_{n-1}] [{\Delta}^{2} S_{n}] - S_{n} [{\Delta} S_{n + 1}] [{\Delta}^{2} S_{n-1}]} {\displaystyle [{\Delta} S_{n-1}] [{\Delta}^{2} S_{n}] - [{\Delta} S_{n + 1}] [{\Delta}^{2} S_{n-1}]} . $$

(A2.11)

Comparison of (A2.7) and (A2.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 (A2.11) as follows:

$$ v_{1}^{(n)} (\beta, S_{n}) = S_{n} + \frac {\displaystyle [{\Delta} S_{n}] [{\Delta} S_{n-1}] [{\Delta}^{2} S_{n}]} {\displaystyle [{\Delta} S_{n-1}] [{\Delta}^{2} S_{n}] - [{\Delta} S_{n + 1}] [{\Delta}^{2} S_{n-1}]}. $$

(A2.12)