ñòð. 4 |

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

(n)

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 + Â· Â· Â·

with

g1 (n + 1)gi (n) âˆ’ g1 (n)gi (n + 1)

(n)

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

(n)

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

(n)

gk; i = (n+1) (n)

gkâˆ’1; k âˆ’ gkâˆ’1; k

(n)

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

(n)

Ek + b1 g1 (n + i) + Â· Â· Â· + bk gk (n + i) = Sn+i ; i = 0; : : : ; k

(n)

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

(n)

by Ek ) gives

Sn Â· Â· Â· Sn+k

g1 (n) Â· Â· Â· g1 (n + k)

. .

. .

. .

gk (n) Â· Â· Â· gk (n + k)

(n)

Ek = :

1 Â·Â·Â· 1

g1 (n) Â· Â· Â· g1 (n + k)

. .

. .

. .

gk (n) Â· Â· Â· gk (n + k)

(n)

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.

(n)

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

(n)

Ek+1 âˆ’ S

lim = 0:

nâ†’âˆž E (n) âˆ’ S

k

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,

e

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

E-algorithm.

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

versatile.

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

(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

(n)

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

2

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)

5

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]

1

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

(n)

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

(n)

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

(n)

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

(n)

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

e

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 |