e

approximants given by

( j)

= [ j + k=k](1) (2.19)

2k

provided (i) zero divisors do not occur in the construction of the epsilon table; and (ii) the corre-

sponding PadÃ approximants identiÿed by (2:19) exist.

e

Proof. The entries W; C; E in the PadÃ table of Figs. 3 and 4 may be taken to correspond to entries

e

( j’1) ( j) ( j+1)

; 2k ; 2k , respectively, in the epsilon table. They neighbour other elements in columns of odd

2k

( j) ( j+1) ( j) ( j’1)

index in the epsilon table, nw := 2k’1 ; ne := 2k’1 ; se := 2k+1 and sw := 2k+1 . By re-pairing, we

have

(nw ’ sw) ’ (ne ’ se) = (nw ’ ne) ’ (sw ’ se): (2.20)

By applying the epsilon algorithm to each term in (2.20), we obtain the compass identity (2.18).

With our conventions, the approximants of type [˜=0] lie in the ÿrst row (m = 0) of the PadÃ e

table. This is quite natural when we regard these approximants as MacLaurin sections of f(z).

However, it must be noted that the row sequence ([˜=m](1); ˜ = m + j; m + j + 1; : : : ; m ÿxed)

( j)

corresponds to the column sequence of entries ( 2m ; j = 0; 1; 2; : : : ; m ÿxed); this identiÿcation follows

from (2.19).

A key property of PadÃ approximants that is an axiom of their deÿnition is that of accuracy-through-

e

order, also called correspondence. Before PadÃ approximants were known as such, attention had

e

rightly been focused on the particular sequence of rational fractions which are truncations of the

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

continued fraction

c0 za1 za2 za3

f(z) = ··· : (2.21)

1’ 1 ’ 1 ’ 1 ’

The right-hand side of (2.21) is called a C-fraction (for instance, see [36]), which is short for

corresponding fraction, and its truncations are called its convergents. Normally, it can be constructed

by successive reciprocation and re-expansion. The ÿrst stage of this process is

1 ’ c0 =f(z) a1 za2 za3

= ··· : (2.22)

1’ 1 ’ 1 ’

z

By undoing this process, we see that the convergents of the C-fraction are rational fractions in the

variable z.

By construction, we see that these convergents agree order by order with f(z), provided all ai = 0,

and this property is called correspondence.

Example 2.2. We truncate (2.21) after a2 and obtain

A2 (z) c0 za1

= : (2.23)

1 ’ 1 ’ za2

B2 (z)

This is a rational fraction of type [1=1], and we take

A2 (z) = c0 (1 ’ za2 ); B2 (z) = 1 ’ z(a1 + a2 ):

Provided all the ai = 0, the convergents of (2.21) are well deÿned. The equality in (2.21) is not

to be understood in the sense of pointwise convergence for each value of z, but in the sense of

correspondence order by order in powers of z.

The numerators and denominators of the convergents of (2.21) are usually constructed using

Euler™s recursion. It is initialised, partly artiÿcially, by

A’1 (z) = 0; A0 (z) = c0 ; B’1 (z) = 1; B0 (z) = 1 (2.24)

and the recursion is

Ai+1 (z) = Ai (z) ’ ai+1 zAi’1 (z); i = 0; 1; 2; : : : ; (2.25)

Bi+1 (z) = Bi (z) ’ ai+1 zBi’1 (z); i = 0; 1; 2; : : : : (2.26)

Euler™s formula is proved in many texts, for example, [1,2,36]. From (2.24) to (2.26), it follows by

induction that

i i+1

˜ = deg{Ai (z)}6 ; m = deg{Bi (z)}6 ; (2.27)

2 2

where [ · ] represents the integer part function and the Baker normalisation is built in:

Bi (0) = 1; i = 0; 1; 2; : : : : (2.28)

The sequence of approximants generated by (2.24) “ (2.26) is shown in Fig. 5.

From (2.19) and (2.27), we see that the convergents of even index i = 2k correspond to PadÃ e

(0)

approximants of type [k=k]; when they are evaluated at z = 1, they are values of 2k on the leading

diagonal of the epsilon table.

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

Fig. 5. A staircase sequence of approximants indexed by i, as in (2.27).

The epsilon algorithm was introduced in (2.1) “ (2.4) as a numerical algorithm. Eq. (2.19) states

its connection with values of certain PadÃ approximants. However, the epsilon algorithm can be

e

given a symbolic interpretation if it is initialised with

j

( j) ( j)

ci z i

= 0; = (2.29)

’1 0

i=0

instead of (2.2) and (2.3). In this case, (2.19) would become

(j)

2k (z) = [ j + k=k](z): (2.30)

The symbolic implementation of the iterative process (2.4) involves considerable cancellation of

polynomial factors, and so we regard this procedure as being primarily of conceptual value.

We have avoided detailed discussions of normality and degeneracy [1,2,25] in this paper so as to

focus on the algorithmic aspects. The case of numerical breakdown associated with zero divisors is

treated by Cordellier [14,15] for example. Refs. [1,2] contain formulae for the di erence between

PadÃ approximants occupying neighbouring positions in the PadÃ table. Using these formulae, one

e e

can show that condition (i) of Theorem 2.1 implies that condition (ii) holds, and so conditions (ii)

can be omitted.

It is always worthwhile to consider the case in which an approximation method gives exact results

at an intermediate stage so that the algorithm is terminated at that stage. For example, let

k

„

f(z) = + (2.31)

0

1 ’ z‚„

„=1

with „ ; ‚„ ∈ C, each |‚„ | ¡ 1, each „ = 0 and all ‚„ distinct. Then f(z) is a rational function

of precise type [k/k]. It is the generating function of the generalised geometric sequence S with

elements

k j+1

1 ’ ‚„

sj = + ; j = 0; 1; 2; : : : : (2.32)

0 „

1 ’ ‚„

„=1

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

This sequence is sometimes called a Dirichlet series and it converges to s∞ = f(1) as j ’ ∞. Its

elements can also be expressed as

k

j

sj = s∞ ’ w „ ‚„ (2.33)

„=1

if

k k

w„ = ‚„ „ (1 ’ ‚„ )’1 :

s∞ = + w„ and

„

„=0 „=1

Then (2.33) expresses the fact that S is composed of exactly k non-trivial, distinct geometric com-

ponents. Theorem 2.1 shows that the epsilon algorithm yields

( j)

= s∞ ; j = 0; 1; 2; : : :

2k

which is the ˜exact result™ in each row of the column of index 2k, provided that zero divisors have

not occurred before this column is constructed. The algorithm should be terminated at this stage via

a consistency test, because zero divisors necessarily occur at the next step. Remarkably, the epsilon

algorithm has some smoothing properties [59], which may (or may not) disguise this problem when

rounding errors occur.

In the next sections, these results will be generalised to the vector case. To do that, we will

also need to consider the paradiagonal sequences of PadÃ approximants given by ([m + J /m](z); m =

e

0; 1; 2; : : : ; J ¿0; J ÿxed). After evaluation at z = 1, we ÿnd that this is a diagonal sequence ( 2m) ; m =

(J

0; 1; 2; : : : ; J ¿0; J ÿxed) in the epsilon table.

3. The vector epsilon algorithm

The epsilon algorithm acquired greater interest when Wynn [57,58] showed that it has a useful

and immediate generalisation to the vector case. Given a sequence

S = (s0 ; s1 ; s2 ; : : : : si ∈ Cd ); (3.1)

the standard implementation of the vector epsilon algorithm (VEA) consists of the following initial-

isation from S followed by its iteration phase:

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

( j)

=0 (artiÿcially); (3.2)

’1

( j)

= sj : (3.3)

0

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

( j) ( j+1) ( j+1) ( j) ’1

= +[ ’ k] : (3.4)

k+1 k’1 k

The iteration formula (3.4) is identical to (2.4) for the scalar case, except that it requires the

speciÿcation of an inverse (reciprocal) of a vector. Usually, the Moore“Penrose (or Samelson)

inverse

d

’1 — H —

|vi |2

C = C =(C C) = C (3.5)

i=1

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

Fig. 6. Columns k = 0; 2 and 4 of the vector epsilon table for Example 3.1 are shown numerically and graphically.

(where the asterisk denotes the complex conjugate and H the Hermitian conjugate) is the most

useful, but there are exceptions [39]. In this paper, the vector inverse is deÿned by (3.5). The vector

epsilon table can then be constructed column by column from (3.2) to (3.4), as in the scalar case,

and as shown in Fig. 6.

Example 3.1. The sequence S is initialised by

s0 := b := (’0:1; 1:5)T (3.6)

(where T denotes the transpose) and it is generated recursively by

sj+1 := b + Gsj ; j = 0; 1; 2; : : : (3.7)

with

0:6 0:5

G= : (3.8)

’1 0:5

The ÿxed point of (3.7) is x = [1; 1], which is the solution of Ax = b with A = I ’ G.

Notice that

( j)

=x for j = 0; 1; 2

4

and this ˜exact™ result is clearly demonstrated in the right-hand columns of Fig. 6.

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

Fig. 7. Schematic view of the two components of u1 (x) and the boundary on the x1 -axis.

This elementary example demonstrates how the VEA can be a powerful convergence accelerator

in an ideal situation. With the same rationale as was explained in the scalar case, the vector epsilon

algorithm is used for sequences of vectors when their convergence is too slow. Likewise, the VEA

can ÿnd an accurate solution (as a ÿxed point of an associated matrix operator) even when the

sequence of vectors is weakly divergent. In applications, these vector sequences usually arise as

sequences of discretised functions, and the operator is a (possibly nonlinear) integral operator. An

example of this kind of vector sequence is one that arises in a problem of current interest. We

consider a problem in acoustics, which is based on a boundary integral equation derived from the

Helmholtz equation [12]. Our particular example includes impedance boundary conditions (3.12)

relevant to the design of noise barriers.

Example 3.2. This is an application of the VEA for the solution of

u(x) = u1 (x) + ik G(x; y)[ÿ(y) ’ 1]u(y) dy (3.9)

for the acoustic ÿeld u(x) at the space point x = (x1 ; x2 ). This ÿeld is conÿned to the half-space

x2 ¿0 by a barrier shown in Fig. 7. The inhomogeneous term in (3.9) is

u1 (x) = eik(x1 sin ‚’x2 cos ‚) + R:eik(x1 sin ‚+x2 cos ‚) (3.10)

which represents an incoming plane wave and a “partially re ected” outgoing plane wave with wave

number k. The re ection coe cient in (3.10) is given by

‚

R = ’tan2 ; (3.11)

2

so that u1 (x) and u(x) satisfy the impedance boundary conditions

@u1 @u

= ’iku1 and = ’ikÿu on : (3.12)

@x2 @x2

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

Notice that u(x1 ; 0)=u1 (x1 ; 0) if ÿ(x1 ; 0) ≡ 1. Then a numerically useful form of the Green™s function

in (3.9) is [10]

eikr ∞

t ’1=2 e’krt (1 + + it)

i (1)

√

G(x; y) = H0 (kr) + dt; (3.13)