. 4
( 83 .)


singularities. In particular, Havie began to study this question under Romberg (Romberg emigrated
to Norway and came to Trondheim in 1949). In 1978, Havie wrote a report, published one year
later [71], where he treated the most general case of an error expansion of the form
S(h) ’ S = a1 g1 (h) + a2 g2 (h) + · · · ;
where S(h) denotes the approximation obtained by the trapezoidal rule with step size h to the deÿnite
integral S and the gi are the known functions (forming an asymptotic sequence when h tends to
zero) appearing in the expansion of the error. Let h0 ¿ h1 ¿ · · · ¿ 0; Sn = S(hn ) and gi (n) = gi (hn ).
Havie set
g1 (n + 1)Sn ’ g1 (n)Sn+1
E1 = :
g1 (n + 1) ’ g1 (n)
Replacing Sn and Sn+1 by their expansions, he obtained
(n) (n) (n)
E1 = S + a2 g1; 2 + a3 g1; 3 + · · ·
g1 (n + 1)gi (n) ’ g1 (n)gi (n + 1)
g1; i = :
g1 (n + 1) ’ g1 (n)
(n) (n)
The same process can be repeated for eliminating g1; 2 in the the expansion of E1 , and so on. Thus,
once again we obtain the E-algorithm
(n+1) (n) (n) (n+1)
gk’1; k Ek’1 ’ gk’1; k Ek’1
Ek = (n+1) (n)
gk’1; k ’ gk’1; k
(n) (n) (n)
with E0 = Sn and g0; i = gi (n). The auxiliary quantities gk; i are recursively computed by the quite
similar rule
(n+1) (n) (n) (n+1)
gk’1; k gk’1; i ’ gk’1; k gk’1; i
gk; i = (n+1) (n)
gk’1; k ’ gk’1; k
with g0; i = gi (n).
C. Brezinski / Journal of Computational and Applied Mathematics 122 (2000) 1“21 11

Havie gave an interpretation of this algorithm in terms of the Gaussian elimination process for
solving the system
Ek + b1 g1 (n + i) + · · · + bk gk (n + i) = Sn+i ; i = 0; : : : ; k
for the unknown Ek .
In 1980, Brezinski took up the same problem, but from the point of view of extrapolation [19].
Let (Sn ) be the sequence to be accelerated. Interpolating it by a sequence of the form Sn = S +
a1 g1 (n) + · · · + ak gk (n), where the gi ™s are known sequences which can depend on the sequence (Sn )
itself, leads to
Sn+i = Sn+i ; i = 0; : : : ; k:
Solving this system directly for the unknown S (which, since it depends on n and k, will be denoted
by Ek ) gives
Sn · · · Sn+k
g1 (n) · · · g1 (n + k)
. .
. .
. .
gk (n) · · · gk (n + k)
Ek = :
1 ··· 1
g1 (n) · · · g1 (n + k)
. .
. .
. .
gk (n) · · · gk (n + k)
Thus Ek is given as a ratio of determinants which is very similar to the ratios previously mentioned.
Indeed, for the choice gi (n)= Sn+i , the ratio appearing in Shanks transformation results while, when
gi (n) = xin , we obtain the ratio expressing the quantities involved in the Richardson extrapolation
process. Other algorithms may be similarly derived.
Now the problem is to ÿnd a recursive algorithm for computing the Ek ™s. Applying Sylvester™s
determinantal identity, Brezinski obtained the two rules of the above E-algorithm. His derivation
of the E-algorithm is closely related to Havie™s since Sylvester™s identity can be proved by using
Gaussian elimination. Brezinski also gave convergence and acceleration results for this algorithm
when the (gi (n)) satisfy certain conditions [19]. These results show that, for accelerating the con-
vergence of a sequence, it is necessary to know the expansion of the error Sn ’ S with respect to
some asymptotic sequence (g1 (n)); (g2 (n)); : : : : The gi (n) are those to be used in the E-algorithm. It
can be proved that, ∀k
Ek+1 ’ S
lim = 0:
n’∞ E (n) ’ S

These results were reÿned by Avram Sidi [134 “136]. Thus the study of the asymptotic expansion of
the error of the sequences to be accelerated is of primary importance, see Walz [144]. For example,
Mohammed Kzaz [79,80] and Pierre Verlinden [142] applied this idea to the problem of accelerating
the convergence of Gaussian quadrature formulae [79] and Pedro Lima and Mario Graca to boundary
value problems with singularities [88,87] (see also the works of Lima and Diogo [87], and Lima and
Carpentier [86]). Other acceleration results were obtained by Matos and Marc Prà vost [94], Prà vost
e e
12 C. Brezinski / Journal of Computational and Applied Mathematics 122 (2000) 1“21

[116] and Pascal Mortreux and Prà vost [102]. An algorithm, more economical than the E-algorithm,
was given by William F. Ford and Avram Sidi [60]. The connection between the E-algorithm and
the -algorithm was studied by Bernhard Beckermann [4]. A general -algorithm connected to the
E-algorithm was given by Carsten Carstensen [48]. See [27] for a more detailed review on the
Convergence acceleration algorithms can also be used for predicting the unknowns terms of a
series or a sequence. This idea, introduced by Jacek Gilewicz [64], was studied by Sidi and Levin
[137], Brezinski [22] and Denis Vekemans [141].

2.5. A new approach

Over the years, a quite general framework was constructed for the theory of extrapolation algo-
rithms. The situation was quite di erent for the practical construction of extrapolation algorithms
and there was little systematic research in their derivation. However, thanks to a formalism due to
Weniger [145], such a construction is now possible, see Brezinski and Matos [38]. It is as follows.
Let us assume that the sequence (Sn ) to be accelerated satisÿes, ∀n; Sn ’ S = an Dn where (Dn ) is a
known sequence, called a remainder (or error) estimate for the sequence (Sn ), and (an ) an unknown
sequence. It is possible to construct a sequence transformation such that its kernel is precisely this
set of sequences. For that purpose, we have to assume that a di erence operator L (that is a linear
mapping of the set of sequences into itself) exists such that ∀n; L(an ) = 0. This means that the
sequence obtained by applying L to the sequence (an ) is identically zero. Such a di erence operator
is called an annihilation operator for the sequence (an ). We have
Sn S
’ = an :
Dn Dn
Applying L and using linearity leads to
Sn 1
L ’ SL = L(an ) = 0:
Dn Dn
We solve for S and designate it by the sequence transformation
L(Sn =Dn )
Tn = :
L(1=Dn )
The sequence (Tn ) is be such that ∀n; Tn = S if and only if ∀n; Sn ’ S = an Dn . This approach is highly
All the algorithms described above and the related devices such as error control, composite se-
quence transformations, least squares extrapolation, etc., can be put into this framework. Moreover,
many new algorithms can be obtained using this approach. The E-algorithm can also be put into this
framework which provides a deeper insight and leads to new properties [41]. Matos [93], using re-
sults from the theory of di erence equations, obtained new and general convergence and acceleration
results when (an ) has an asymptotic expansion of a certain form.

2.6. Integrable systems

The connection between convergence acceleration algorithms and discrete integrable systems is a
subject whose interest is rapidly growing among physicists. When a numerical scheme is used for
C. Brezinski / Journal of Computational and Applied Mathematics 122 (2000) 1“21 13

integrating a partial di erential evolution equation, it is important that it preserves the quantities that
are conserved by the partial di erential equation itself. An important character is the integrability
of the equation. Although this term has not yet received a completely satisfactory deÿnition (see
[66]), it can be understood as the ability to write the solution explicitly in terms of a ÿnite number
of functions or as the conÿnement of singularities in ÿnite domains. The construction of integrable
discrete forms of integrable partial di erential equations is highly nontrivial. A major discovery in
the ÿeld of integrability was the occurrence of a solitary wave (called a soliton) in the Korteweg“
de Vries (KdV) equation. Integrability is a rare phenomenon and the typical dynamical system is
nonintegrable. A test of integrability, called singularity conÿnement, was given by B. Grammaticos,
A. Ramani and V. Papageorgiou [67]. It turns out that this test is related to the existence of singular
rules for avoiding a division by zero in convergence acceleration algorithms (see Section 1.2).
The literature on this topic is vast and we cannot enter into the details of it. We only want to
give an indication of the connection between these two subjects since both domains could beneÿt
from it.
In the rule for the -algorithm, V. Papageorgiou, B. Grammaticos and A. Ramani set m = k + n
and replaced k by u(n; m) + mp + nq, where p and q satisfy p2 ’ q2 = 1. They obtained [111]
[p ’ q + u(n; m + 1) ’ u(n + 1; m)][p + q + u(n + 1; m + 1) ’ u(n; m)] = p2 ’ q2 :
This is the discrete lattice KdV equation. Since this equation is integrable, one can expect integrability
to hold also for the -algorithm, and, thanks to the singular rules of Wynn and Cordellier mentioned
at the end of Subsection 1.2, this is indeed the case.
In the rule of the -algorithm, making the change of variable k = t= 3 and n ’ 1=2 = x= ’ ct= 3
and replacing k by p + 2 u(x ’ =2; t) where c and p are related by 1 ’ 2c = 1=p2 , A. Nagai and
J. Satsuma obtained [105]
1 1
u(x ’ =2 + c ; t + 3 ) ’ 2 u(x + =2 ’ c ; t ’ 3 ) = ’ :
p + 2 u(x + =2; t) p + 2 u(x ’ =2; t)
We have, to terms of order , the KdV equation
1 1
(1 ’ p’4 )uxxx = 0:
ut ’ uux +
p3 48p2
Other discrete numerical algorithms, such as the qd, LR, and -algorithms are connected to other
discrete or continuous integrable equations (see, for example [112]). Formal orthogonal polynomials,
continued fractions, Padà approximation also play a rˆ le in this topic [113].
e o
By replacing the integer n in the -algorithm by a continuous variable, Wynn derived the con uent
form of the -algorithm [149]
k+1 (t) = k’1 (t) +
k (t)

with ’1 (t) ≡ 0 and 0 (t) = f(t). This algorithm is the continuous counterpart of the -algorithm
and its aim is to compute limt’∞ f(t). Setting Nk (t) = k (t) k+1 (t), A. Nagai, T. Tokihiro and
J. Satsuma [106] obtained
Nk (t) = Nk (t)[Nk’1 (t) ’ Nk+1 (t)]:
The above equation is the Backlund transformation of the discrete Toda molecule equation [139].
14 C. Brezinski / Journal of Computational and Applied Mathematics 122 (2000) 1“21

So, we see that some properties of integrable systems are related to properties of convergence
acceleration algorithms. On the other hand, discretizing integrable partial di erential equations leads
to new sequence transformations which have to be studied from the point of view of their algebraic
and acceleration properties. Replacing the second integer k in the con uent form of the -algorithm by
a continuous variable, Wynn obtained a partial di erential equation [152]. Its relation with integrable
systems is an open question.
The connection between integrable systems and convergence acceleration algorithms needs to be
investigated in more details to fully understand its meaning which is not clear yet.

3. The vector case

In numerical analysis, many iterative methods lead to vector sequences. To accelerate the conver-
gence of such sequences, it is always possible to apply a scalar algorithm componentwise. However,
vector sequence transformations, specially built for that purpose, are usually more powerful. The
ÿrst vector algorithm to be studied was the vector -algorithm. It was obtained by Wynn [150] by
(n) (n)
replacing, in the rule of the scalar -algorithm, 1= k by ( k )’1 where the inverse y’1 of a
vector y is deÿned by y’1 = y=(y; y). Thus, with this deÿnition, the rule of the -algorithm can be
applied to vector sequences. Using Cli ord algebra, J.B. McLeod proved in 1971 [96] that ∀n; 2k =S
if the sequence (Sn ) satisÿes a0 (Sn ’ S) + · · · + ak (Sn+k ’ S) = 0; ∀n with a0 ak = 0; a0 + · · · + ak = 0.
This result is only valid for real sequences (Sn ) and real ai ™s. Moreover, contrary to the scalar case,
this condition is only su cient. In 1983, Peter R. Graves“Morris [68] extended this result to the
complex case using a quite di erent approach.
A drawback to the development of the theory of the vector -algorithm was that it was not known
whether a corresponding generalization of Shanks transformation was underlying the algorithm, that
is whether the vectors k obtained by the algorithm could be expressed as ratios of determinants
(or some kind of generalization of determinants). This is why Brezinski [17], following the same
path as Shanks, tried to construct a vector sequence transformation with the kernel a0 (Sn ’ S) + · · · +
ak (Sn+k ’ S) = 0. He obtained a transformation expressed as a ratio of determinants. He then had
to develop a recursive algorithm for avoiding their computation. This was the so-called topological
-algorithm. This algorithm has many applications, in particular, to the solution of systems of linear
equations (it is related to the biconjugate gradient algorithm [18, pp. 185 ]). In the case of a
system of nonlinear equations, it gave rise to a generalization of Ste ensen™s method [13]. That
algorithm has a quadratic convergence under some assumptions as established by Hervà Le Ferrande
[83] following the ideas presented by Khalide Jbilou and Sadok [75]. The denominator of the vector
2k obtained by the vector -algorithm was ÿrst written as a determinant of dimension 2k + 1 by
Graves“Morris and Chris Jenkins in [69]. The numerator follows immediately by modifying the ÿrst
row of the denominator, a formula given by Ahmed Salam and Graves“Morris [126]. However, the
dimension of the corresponding determinants in the scalar case is only k + 1. It was proved by
Salam [124] that the vectors 2k computed by the vector -algorithm can be expressed as a ratio of
two designants of dimension k + 1. A designant is a generalization of a determinant when solving
a system of linear equations in a noncommutative algebra. An algebraic approach to this algorithm
was given in [125]. This approach, which involves the use of a Cli ord algebra, was used in [45]
for extending the mechanism given in [41] to the vector and matrix cases. The vector generalization
C. Brezinski / Journal of Computational and Applied Mathematics 122 (2000) 1“21 15

of the E-algorithm [19] can be explained similarly. This algorithm makes use of a ÿxed vector
y. Jet Wimp [146, pp. 176 “177] generalized it using a sequence (yn ) instead of y. Jeannette van
Iseghem [140] gave an algorithm for accelerating vector sequences based on the vector orthogonal
polynomials she introduced for generalizing Padà approximants to the vector case. Other vector
sequence transformations are due to Osada [109] and Jbilou and Sadok [76]. Benchiboun [6] and
Abderrahim Messaoudi [100] studied matrix extrapolation algorithms.
We have seen that, in the scalar case, the kernels of sequence transformations may be expressed
as relationships with constant coe cients. This is also the case for the vector and the topological
-algorithms and the vector E-algorithm. The ÿrst (and, to my knowledge, only) transformation
treating a relationship with varying coe cients was introduced in [42]. The theory developed there
also explains why the case of a relationship with non-constant coe cients is a di cult problem in
the scalar case and why it could be solved, on the contrary, in the vector case. The reason is that
the number of unknown coe cients appearing in the expression of the kernel must be strictly less
than the dimension of the vectors. Brezinski in [34] proposed a general methodology for constructing
vector sequence transformations. It leads to a uniÿed presentation of several approaches to the subject


. 4
( 83 .)