<< Предыдущая стр. 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 ), 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 , 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 . 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 . 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 
eikr в€ћ
t в€’1=2 eв€’krt (1 + + it)
i (1)
в€љ
G(x; y) = H0 (kr) + dt; (3.13)
 << Предыдущая стр. 13(из 83 стр.)ОГЛАВЛЕНИЕ Следующая >>