<<

. 13
( 83 .)



>>

Theorem 2.1. The entries in columns of even index in the epsilon table are values of PadÃ
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)

<<

. 13
( 83 .)



>>