ñòð. 12 |

includes the complex case wherever possible. c 2000 Elsevier Science B.V. All rights reserved.

Keywords: Epsilon algorithm; qd algorithm; PadÃƒ ; Vector-valued approximant; Wynn; Cross rule; Star identity; Compass

e

identity; Designant

1. Introduction

A sequence with a limit is as basic a topic in mathematics as it is a useful concept in science

and engineering. In the applications, it is usually the limit of a sequence, or a Ã¿xed point of its

generator, that is required; the existence of the limit is rarely an issue, and rapidly convergent

sequences are welcomed. However, if one has to work with a sequence that converges too slowly,

the epsilon algorithm is arguably the best all-purpose method for accelerating its convergence. The

algorithm was discovered by Wynn [54] and his review article [59] is highly recommended. The

epsilon algorithm can also be used for weakly diverging sequences, and for these the desired limit

is usually deÃ¿ned as being a Ã¿xed point of the operator that generates the sequence. There are

interesting exceptional cases, such as quantum well oscillators [51], where the epsilon algorithm is

not powerful enough and we refer to the companion paper by Homeier [33] in which the more

âˆ—

Corresponding author.

E-mail address: p.r.graves-morris@bradford.ac.uk (P.R. Graves-Morris).

0377-0427/00/$ - see front matter c 2000 Elsevier Science B.V. All rights reserved.

PII: S 0 3 7 7 - 0 4 2 7 ( 0 0 ) 0 0 3 5 5 - 1

52 P.R. Graves-Morris et al. / Journal of Computational and Applied Mathematics 122 (2000) 51â€“80

powerful Levin-type algorithms, etc., are reviewed. The connections between the epsilon algorithm

and similar algorithms are reviewed by Weniger [50,52].

This paper is basically a review of the application of the epsilon algorithm, with an emphasis on

the case of complex-valued, vector-valued sequences. There are already many reviews and books

which include sections on the scalar epsilon algorithm, for example [1,2,9,17,53]. In the recent past,

there has been progress with the problem of numerical breakdown of the epsilon algorithm. Most

notably, Cordellierâ€™s algorithm deals with both scalar and vector cases [13â€“16]. This work and its

theoretical basis has been extensively reviewed [26,27]. In this paper, we focus attention on how

the epsilon algorithm is used for sequences (si ) in which si âˆˆ Cd . The case d = 1 is the scalar case,

and the formulation for si âˆˆ C is essentially the same as that for si âˆˆ R. Not so for the vector case,

and we give full details of how the vector epsilon and vector qd algorithms are implemented when

si âˆˆ Cd , and of the connections with vector PadÃƒ approximation. Understanding these connections

e

is essential for specifying the range of validity of the methods. Frequently, the word â€œnormallyâ€

appears in this paper to indicate that the results may not apply in degenerate cases. The adaptations

for the treatment of degeneracy are almost the same for both real and complex cases, and so we

refer to [25 â€“27] for details.

In Section 2, we formulate the epsilon algorithm, and we explain its connection with PadÃƒ e

approximation and the continued fractions called C-fractions. We give an example of how the

epsilon algorithm works in ideal circumstances, without any signiÃ¿cant loss of numerical precision

(which is an unusual outcome).

In Section 3, we formulate the vector epsilon algorithm, and we review its connection with

vector-valued PadÃƒ approximants and with vector-valued C-fractions. There are two major general-

e

isations of the scalar epsilon algorithm to the vector case. One of them is Brezinskiâ€™s topological

epsilon algorithm [5,6,35,48,49]. This algorithm has two principal forms, which might be called the

forward and backward versions; and the backward version has the orthogonality properties associated

with Lanczos methods [8]. The denominator polynomials associated with all forms of the topological

epsilon algorithm have degrees which are the same as those for the scalar case [2,5,8]. By contrast,

the other generalisation of the scalar epsilon algorithm to the vector case can be based on using

generalised inverses of vectors, and it is this generalisation which is the main topic of this paper. We

illustrate how the vector epsilon algorithm works in a two-dimensional real space, and we give a re-

alistic example of how it works in a high-dimensional complex space. The denominator polynomials

used in the scalar case are generalised both to operator polynomials of the same degree and to scalar

polynomials of double the degree in the vector case, and we explain the connections between these

twin generalisations. Most of the topics reviewed in Section 3 have a direct generalisation to the

rational interpolation problem [25]. We also note that the method of GIPAs described in Section 3

generalises directly to deal with sequences of functions in L2 (a; b) rather than vectors Cd ; in this

sense, the vectors are regarded as discretised functions [2].

In Section 4 we review the use of the vector qd algorithm for the construction of vector-valued

C-fractions, and we note the connections between vector orthogonal polynomials and the vector

epsilon algorithm. We prove the cross-rule (4.18), (4.22) using a Cli ord algebra. For real-valued

vectors, we observe that it is really an overlooked identity amongst Hankel designants. Here, the

Cross Rule is proved as an identity amongst complex-valued vectors using Mooreâ€“Penrose inverses.

The importance of studying the vector epsilon algorithm lies partly in its potential [20] for

application to the acceleration of convergence of iterative solution of discretised PDEs. For

P.R. Graves-Morris et al. / Journal of Computational and Applied Mathematics 122 (2000) 51â€“80 53

example, Gaussâ€“Seidel iteration generates sequences of vectors which often converge too slowly

to be useful. SOR, multigrid and Lanczos methods are alternative approaches to the problem which

are currently popular, but the success of the techniques like CGS and LTPMs (see [31] for an expla-

nation of the techniques and the acronyms) indicates the need for continuing research into numerical

methods for the acceleration of convergence of vector-valued sequences.

To conclude this introductory section, we recall that all algorithms have their domains of validity.

The epsilon algorithm fails for logarithmically convergent sequences (which converge too slowly) and

it fails to Ã¿nd the Ã¿xed point of the generator of sequences which diverge too fast. For example, if

C

+ O(nâˆ’2 ); C = 0;

sn âˆ’ s =

n

the sequence (sn ) is logarithmically convergent to s. More precisely, a sequence is deÃ¿ned to con-

verge logarithmically to s if it converges to s at a rate governed by

sn+1 âˆ’ s

lim = 1:

nâ†’âˆž sn âˆ’ s

Not only does the epsilon algorithm usually fail for such sequences, but Delahaye and Germain-Bonne

[18,19] have proved that there is no universal accelerator for logarithmically convergent

sequences.

Reviews of series transformations, such as those of the energy levels of the quantum-mechanical

harmonic oscillator [21,50,51], and of the Riemann zeta function [34], instructively show the in-

adequacy of the epsilon algorithm when the series coe cients diverge too fast. Information about

the asymptotic form of the coe cients and scaling properties of the solution is exploited to create

purpose-built acceleration methods. Exotic applications of the -algorithm appear in [55].

2. The epsilon algorithm

The epsilon algorithm was discovered by Wynn [54] as an e cient implementation of Shanksâ€™

method [47]. It is an algorithm for acceleration of convergence of a sequence

S = (s0 ; s1 ; s2 ; : : : ; si âˆˆ C) (2.1)

and it comprises the following initialisation and iterative phases:

Initialisation: For j = 0; 1; 2; : : :

( j)

= 0 (artiÃ¿cially); (2.2)

âˆ’1

( j)

= sj : (2.3)

0

Iteration: For j; k = 0; 1; 2; : : :

( j) ( j+1) ( j+1) ( j) âˆ’1

= +[ âˆ’ k] : (2.4)

k+1 kâˆ’1 k

The entries k j) are displayed in the epsilon table on the left-hand side of Fig. 1, and the initialisation

(

has been built in.

54 P.R. Graves-Morris et al. / Journal of Computational and Applied Mathematics 122 (2000) 51â€“80

Fig. 1. The epsilon table, and a numerical example of it.

Example 2.1. Gregoryâ€™s series for tanâˆ’1 z is

z3 z5 z7

tanâˆ’1 z = z âˆ’ + âˆ’ + Â·Â·Â· : (2.5)

3 5 7

This series can be used to determine the value of by evaluating its MacLaurin sections at z = 1:

2j+1

sj := [4 tanâˆ’1 (z)]0 ; j = 0; 1; 2; : : : : (2.6)

z=1

Nuttallâ€™s notation is used here and later on. For a function whose MacLaurin series is

2

(z) = + 1z + 2z + Â·Â·Â·;

0

its sections are deÃ¿ned by

k

(z)]k i

[ = iz for 06j6k: (2.7)

j

i=j

In fact, sj â†’ as j â†’ âˆž [2] but sequence (2.6) converges slowly, as is evidenced in the column

k = 0 of entries sj = 0j) in Fig. 1. The columns of odd index have little signiÃ¿cance, whereas the

(

columns of even index can be seen to converge to , which is the correct limit [2], increasingly

( j)

fast, as far as the table goes. Some values of 2k are also shown on the bar chart (Fig. 2). Notice

(2) (0)

that 2 = 3:145 and 4 = 3:142 cannot be distinguished visually on this scale.

In Example 2.1, convergence can be proved and the rate of convergence is also known [2]. From

the theoretical viewpoint, Example 2.1 is ideal for showing the epsilon algorithm at its best. It

is noticeable that the entries in the columns of odd index are large, and this e ect warns us to

beware of possible loss of numerical accuracy. Like all algorithms of its kind (which use reciprocal

di erences of convergent sequences) the epsilon algorithm uses (and usually uses up) numerical

precision of the data to do its extrapolation. In this case, there is little loss of numerical precision

(0)

using 16 decimal place (MATLAB) arithmetic, and 22 = almost to machine precision. In this

case, the epsilon algorithm converges with great numerical accuracy because series (2.5) is a totally

oscillating series [4,7,17,59].

To understand in general how and why the epsilon algorithm converges, whether we are referring

( j) ( j)

to its even columns ( 2k ; j = 0; 1; 2; : : : ; k Ã¿xed) or its diagonals ( 2k ; k = 0; 1; 2; : : : ; j Ã¿xed) or any

P.R. Graves-Morris et al. / Journal of Computational and Applied Mathematics 122 (2000) 51â€“80 55

( j)

Fig. 2. Values of 2k for Example 2.1, showing the convergence rate of the epsilon algorithm using n + 1 = 1; 2; 3; 4; 5

terms of the given sequence.

other sequence, the connection with PadÃƒ approximation is essential [1,2,56]. Given a (possibly

e

formal) power series

f(z) = c0 + c1 z + c2 z 2 + Â· Â· Â· ; (2.8)

the rational function

A(z)B(z)âˆ’1 â‰¡ [â€˜=m](z) (2.9)

is deÃ¿ned as a PadÃƒ approximant for f(z) of type [â€˜=m] if

e

(i) deg{A(z)}6â€˜; deg{B(z)}6m; (2.10)

(ii) f(z)B(z) âˆ’ A(z) = O(z â€˜+m+1 ); (2.11)

(iii) B(0) = 0: (2.12)

The Baker condition

B(0) = 1 (2.13)

is often imposed for reliability in the sense of (2.14) below and for a deÃ¿nite speciÃ¿cation of A(z)

and B(z). The deÃ¿nition above contrasts with the classical (Frobenius) deÃ¿nition in which axiom

(iii) is waived, and in this case the existence of A(z) and B(z) is guaranteed, even though (2.14)

below is not. Using speciÃ¿cation (2.10) â€“ (2.13), we Ã¿nd that

f(z) âˆ’ A(z)B(z)âˆ’1 = O(z â€˜+m+1 ); (2.14)

56 P.R. Graves-Morris et al. / Journal of Computational and Applied Mathematics 122 (2000) 51â€“80

Fig. 3. Relative location of PadÃƒ approximants.

e

provided that a solution of (2.15) below can be found. To Ã¿nd B(z), the linear equations corre-

sponding to accuracy-through-orders z â€˜+1 ; z â€˜+2 ; : : : ; z â€˜+m in (2.11) must be solved. They are

ï£® ï£¹ï£® ï£¹ ï£® ï£¹

câ€˜âˆ’m+1 ::: câ€˜ bm câ€˜+1

. . ï£ºï£¯ . ï£º ï£¯.ï£º

ï£¯ . . ï£»ï£° . ï£» = âˆ’ï£° . ï£»: (2.15)

ï£° . . . .

câ€˜ ::: câ€˜+mâˆ’1 b1 câ€˜+m

The coe cients of B(z) = m bi z i are found using an accurate numerical solver of (2.15). By

i=0

contrast, for purely theoretical purposes, Cramerâ€™s rule is applied to (2.15). We are led to deÃ¿ne

câ€˜âˆ’m+1 câ€˜âˆ’m+2 ::: câ€˜+1

câ€˜âˆ’m+2 câ€˜âˆ’m+3 ::: câ€˜+2

. . .

. . .

q[â€˜=m] (z) = (2.16)

. . .

câ€˜ câ€˜+1 ::: câ€˜+m

zm z mâˆ’1 ::: 1

and then we Ã¿nd that

B[â€˜=m] (z) = q[â€˜=m] (z)=q[â€˜=m] (0) (2.17)

is the denominator polynomial for the PadÃƒ approximation problem (2.9) â€“ (2.15) provided that

e

[â€˜=m]

q (0) = 0.

The collection of PadÃƒ approximants is called the PadÃƒ table, and in Fig. 3 we show Ã¿ve neigh-

e e

bouring approximants in the table.

These approximants satisfy a Ã¿ve-point star identity,

[N (z) âˆ’ C(z)]âˆ’1 + [S(z) âˆ’ C(z)]âˆ’1 = [E(z) âˆ’ C(z)]âˆ’1 + [W (z) âˆ’ C(z)]âˆ’1 ; (2.18)

called Wynnâ€™s identity or the compass identity. The proof of (2.18) is given in [1,2], and it is

also a corollary (in the case d = 1) of the more general result (3.59) that we prove in the next

section. Assuming (2.18) for the moment, the connection between PadÃƒ approximation and the epsilon

e

P.R. Graves-Morris et al. / Journal of Computational and Applied Mathematics 122 (2000) 51â€“80 57

Fig. 4. Some artiÃ¿cial entries in the PadÃƒ table are shown circled.

e

algorithm is given by connecting the coe cients of f(z) with those of S with

c0 = s0 ; ci = si âˆ’ siâˆ’1 ; i = 1; 2; 3; : : : ;

and by

ñòð. 12 |