Table Of ContentJournal of Computational and Applied Mathematics
Copyright © 2003 Elsevier B.V. All rights reserved
Volume 122, Issues 1-2, Pages 1-357 (1 October 2000)
Display Checked Docs | E-mail Articles | Export Citations View: Citations
Convergence acceleration during the 20th century, Pages 1-21
1. gfedc
C. Brezinski
SummaryPlus | Full Text + Links | PDF (139 K)
On the history of multivariate polynomial interpolation, Pages 23-35
2. gfedc
Mariano Gasca and Thomas Sauer
SummaryPlus | Full Text + Links | PDF (109 K)
Elimination techniques: from extrapolation to totally positive matrices and CAGD, Pages 37-50
3. gfedc
M. Gasca and G. Mühlbach
Abstract | PDF (116 K)
The epsilon algorithm and related topics, Pages 51-80
4. gfedc
P. R. Graves-Morris, D. E. Roberts and A. Salam
SummaryPlus | Full Text + Links | PDF (256 K)
Scalar Levin-type sequence transformations, Pages 81-147
5. gfedc
Herbert H. H. Homeier
Abstract | PDF (428 K)
Vector extrapolation methods. Applications and numerical comparison, Pages 149-165
6. gfedc
K. Jbilou and H. Sadok
SummaryPlus | Full Text + Links | PDF (125 K)
Multivariate Hermite interpolation by algebraic polynomials: A survey, Pages 167-201
7. gfedc
R. A. Lorentz
SummaryPlus | Full Text + Links | PDF (222 K)
Interpolation by Cauchy–Vandermonde systems and applications, Pages 203-222
8. gfedc
G. Mühlbach
SummaryPlus | Full Text + Links | PDF (165 K)
The E-algorithm and the Ford–Sidi algorithm, Pages 223-230
9. gfedc
Naoki Osada
SummaryPlus | Full Text + Links | PDF (76 K)
Diophantine approximations using Padé approximations, Pages 231-250
10. gfedc
M. Prévost
SummaryPlus | Full Text + Links | PDF (157 K)
The generalized Richardson extrapolation process GREP(1) and computation of derivatives of limits of sequences
11. gfedc
with applications to the d(1)-transformation, Pages 251-273
Avram Sidi
SummaryPlus | Full Text + Links | PDF (185 K)
Matrix Hermite–Padé problem and dynamical systems, Pages 275-295
12. gfedc
Vladimir Sorokin and Jeannette Van Iseghem
SummaryPlus | Full Text + Links | PDF (156 K)
Numerical analysis of the non-uniform sampling problem, Pages 297-316
13. gfedc
Thomas Strohmer
SummaryPlus | Full Text + Links | PDF (523 K)
Asymptotic expansions for multivariate polynomial approximation, Pages 317-328
14. gfedc
Guido Walz
SummaryPlus | Full Text + Links | PDF (95 K)
Prediction properties of Aitken's iterated 2 process, of Wynn's epsilon algorithm, and of Brezinski's iterated
15. gfedc
theta algorithm, Pages 329-356
Ernst Joachim Weniger
SummaryPlus | Full Text + Links | PDF (189 K)
Index, Page 357
16. gfedc
Unknown
PDF (31 K)
Numerical Analysis 2000 Vol. II: Interpolation and extrapolation, Pages ix-xi
17. gfedc
C. Brezinski
SummaryPlus | Full Text + Links | PDF (34 K)
Journal of Computational and Applied Mathematics 122 (2000) ix{xi
www.elsevier.nl/locate/cam
Preface
Numerical Analysis 2000
Vol. II: Interpolation and extrapolation
C. Brezinski
Laboratoire d’Analyse Nume(cid:19)rique et d’Optimisation, Universite(cid:19) des Sciences et Technologies de Lille,
59655 Villeneuve d’Ascq Cedex, France
This volume is dedicated to two closely related subjects: interpolation and extrapolation. The
papers can be divided into three categories: historical papers, survey papers and papers presenting
new developments.
Interpolation is an old subject since, as noticed in the paper by M. Gasca and T. Sauer, the
term was coined by John Wallis in 1655. Interpolation was the (cid:12)rst technique for obtaining an
approximation of a function. Polynomial interpolation was then used in quadrature methods and
methods for the numerical solution of ordinary di(cid:11)erential equations.
Obviously, some applications need interpolation by functions more complicated than polynomials.
The case of rational functions with prescribed poles is treated in the paper by G. Mu(cid:127)hlbach. He
gives a survey of interpolation procedures using Cauchy{Vandermonde systems. The well-known
formulae of Lagrange, Newton and Neville{Aitken are generalized. The construction of rational
B-splines is discussed.
Trigonometric polynomials are used in the paper by T. Strohmer for the reconstruction of a signal
from non-uniformly spaced measurements. They lead to a well-posed problem that preserves some
important structural properties of the original in(cid:12)nite dimensional problem.
More recently, interpolation in several variables was studied. It has applications in (cid:12)nite di(cid:11)erences
and (cid:12)nite elements for solving partial di(cid:11)erential equations. Following the pioneer works of P. de
Casteljau and P. Be(cid:19)zier, another very important domain where multivariate interpolation plays a
fundamental role is computer-aided geometric design (CAGD) for the approximation of surfaces.
The history of multivariate polynomial interpolation is related in the paper by M. Gasca and T.
Sauer.
The paper by R.A. Lorentz is devoted to the historical development of multivariate Hermite
interpolation by algebraic polynomials.
In his paper, G. Walz treats the approximation of multivariate functions by multivariate Bernstein
polynomials. An asymptotic expansion of these polynomials is given and then used for building, by
extrapolation, a new approximation method which converges much faster.
E-mail address: [email protected] (C. Brezinski).
0377-0427/00/$-see front matter (cid:13)c 2000 Elsevier Science B.V. All rights reserved.
PII: S0377-0427(00)00352-6
x Preface/Journal of Computational and Applied Mathematics 122 (2000) ix{xi
Extrapolation is based on interpolation. In fact, extrapolation consists of interpolation at a point
outside the interval containing the interpolation points. Usually, this point is either zero or in(cid:12)nity.
Extrapolation is used in numerical analysis to improve the accuracy of a process depending of
a parameter or to accelerate the convergence of a sequence. The most well-known extrapolation
processes are certainly Romberg’s method for improving the convergence of the trapezoidal rule for
the computation of a de(cid:12)nite integral and Aitken’s (cid:1)2 process which can be found in any textbook
of numerical analysis.
An historical account of the development of the subject during the 20th century is given in the
paper by C. Brezinski.
The theory of extrapolation methods lays on the solution of the system of linear equations corre-
sponding to the interpolation conditions. In their paper, M. Gasca and G. Mu(cid:127)hlbach show, by using
elimination techniques, the connection between extrapolation, linear systems, totally positive matrices
and CAGD.
There exist many extrapolation algorithms. From a (cid:12)nite section S ;:::;S of the sequence (S ),
n n+k n
they built an improved approximation of its limit S. This approximation depends on n and k. When
at least one of these indexes goes to in(cid:12)nity, a new sequence is obtained with, possibly, a faster
convergence.
In his paper, H.H.H. Homeier studies scalar Levin-type acceleration methods. His approach is based
on the notion of remainder estimate which allows to use asymptotic information on the sequence to
built an e(cid:14)cient extrapolation process.
The most general extrapolation process known so far is the sequence transformation known under
the name of E-algorithm. It can be implemented by various recursive algorithms. In his paper,
N. Osada proved that the E-algorithm is mathematically equivalent to the Ford{Sidi algorithm. A
slightly more economical algorithm is also proposed.
When S depends on a parameter t, some applications need the evaluation of the derivative of S
with respect to t. A generalization of Richardson extrapolation process for treating this problem is
considered in the paper by A. Sidi.
Instead of being used for estimating the limit S of a sequence from S ;:::;S , extrapolation
n n+k
methods can also be used for predicting the next unknown terms S ;S ;:::. The prediction
n+k+1 n+k+2
properties of some extrapolation algorithms are analyzed in the paper by E.J. Weniger.
Quite often in numerical analysis, sequences of vectors have to be accelerated. This is, in particular,
the case in iterative methods for the solution of systems of linear and nonlinear equations.
Vector acceleration methods are discussed in the paper by K. Jbilou and H. Sadok. Using projec-
tors, they derive a di(cid:11)erent interpretation of these methods and give some theoretical results. Then,
various algorithms are compared when used for the solution of large systems of equations coming
out from the discretization of partial di(cid:11)erential equations.
Another point of view is taken in the paper by P.R. Graves-Morris, D.E. Roberts and A. Salam.
After reminding, in the scalar case, the connection between the (cid:15)-algorithm, Pade(cid:19) approximants and
continued fractions, these authors show that the vector (cid:15)-algorithm is the best all-purpose algorithm
for the acceleration of vector sequences.
There is a subject which can be related either to interpolation (more precisely, Hermite interpo-
lation by a rational function at the point zero) and to convergence acceleration: it is Pade(cid:19) approx-
imation. Pade(cid:19) approximation is strongly connected to continued fractions, one of the oldest subject
in mathematics since Euclid g.c.d. algorithm is an expansion into a terminating continued fraction.
Preface/Journal of Computational and Applied Mathematics 122 (2000) ix{xi xi
Although they were implicitly known before, Pade(cid:19) approximants were really introduced by Johan
Heinrich Lambert in 1758 and Joseph Louis Lagrange in 1776. Pade(cid:19) approximants have impor-
tant applications in many branches of applied sciences when the solution of a problem is obtained
as a power series expansion and some of its properties have to be guessed from its (cid:12)rst Taylor
coe(cid:14)cients. In this volume, two papers deal with nonclassical applications of Pade(cid:19) approximation.
M. Pre(cid:19)vost shows how Pade(cid:19) approximants can be used to obtain Diophantine approximations
of real and complex numbers and then proving irrationality. Pade(cid:19) approximation of the asymptotic
expansion of the remainder of a series also provides Diophantine approximations.
The solution of a discrete dynamical system can be related to matrix Hermite{Pade(cid:19) approximants,
an approach developed in the paper by V. Sorokin and J. van Iseghem. Spectral properties of the
band operator are investigated. The inverse spectral method is used for the solution of dynamical
systems de(cid:12)ned by a Lax pair.
Obviously, all aspects of interpolation and extrapolation have not been treated in this volume.
However, many important topics have been covered.
I would like to thank all authors for their e(cid:11)orts.
Journal of Computational and Applied Mathematics 122 (2000) 1{21
www.elsevier.nl/locate/cam
Convergence acceleration during the 20th century
C. Brezinski
Laboratoire d’Analyse Nume(cid:19)rique et d’Optimisation, UFR IEEA, Universite(cid:19) des Sciences et Technologies de Lille,
59655{Villeneuve d’Ascq cedex, France
Received 8 March 1999; received in revised form 12 October 1999
1. Introduction
In numerical analysis many methods produce sequences, for instance iterative methods for solving
systems of equations, methods involving series expansions, discretization methods (that is methods
depending on a parameter such that the approximate solution tends to the exact one when the pa-
rameter tends to zero), perturbation methods, etc. Sometimes, the convergence of these sequences
is slow and their e(cid:11)ective use is quite limited. Convergence acceleration methods consist of trans-
forming a slowly converging sequence (S ) into a new sequence (T ) converging to the same limit
n n
faster than the initial one. Among such sequence transformations, the most well known are certainly
Richardson’s extrapolation algorithm and Aitken’s (cid:1)2 process. All known methods are constructed
by extrapolation and they are often called extrapolation methods. The idea consists of interpolating
the terms S ;S ;:::;S of the sequence to be transformed by a sequence satisfying a certain
n n+1 n+k
relationship depending on parameters. This set of sequences is called kernel of the transformation
and every sequence of this set is transformed into a constant sequence by the transformation into
consideration. For example, as we will see below, the kernel of Aitken’s (cid:1)2 process is the set of
sequences satisfying 8n;a (S −S)+a (S −S)=0, where a and a are parameters such that
0 n 1 n+1 0 1
a +a 6= 0. If Aitken’s process is applied to such a sequence, then the constant sequence (T =S)
0 1 n
is obtained. The parameters involved in the de(cid:12)nition of the kernel are uniquely determined by the
interpolation conditions and then the limit of the interpolating sequence of the kernel is taken as an
approximation of the limit of the sequence to be accelerated. Since this limit depends on the index
n, it will be denoted by T . E(cid:11)ectively, the sequence (S ) has been transformed into a new sequence
n n
(T ).
n
This paper, which is based on [31], but includes new developments obtained since 1995, presents
my personal views on the historical development of this subject during the 20th century. I do not
pretend to be exhaustive nor even to quote every important contribution (if a reference does not
E-mail address: [email protected] (C. Brezinski)
0377-0427/00/$-see front matter (cid:13)c 2000 Elsevier Science B.V. All rights reserved.
PII: S0377-0427(00)00360-5
2 C. Brezinski/Journal of Computational and Applied Mathematics 122 (2000) 1{21
appear below, it does not mean that it is less valuable). I refer the interested reader to the literature
and, in particular to the recent books [55,146,33,144]. For an extensive bibliography, see [28].
I will begin with scalar sequences and then treat the case of vector ones. As we will see, a
sequence transformation able to accelerate the convergence of all scalar sequences cannot exist.
Thus, it is necessary to obtain many di(cid:11)erent convergence acceleration methods, each being suitable
for a particular class of sequences. Many authors have studied the properties of these procedures and
proved some important classes of sequences to be accelerable by a given algorithm. Scalar sequence
transformations have also been extensively studied from the theoretical point of view.
The situation is more complicated and more interesting for vector sequences. In the case of a
sequence of vectors, it is always possible to apply a scalar acceleration procedure componentwise.
However, such a strategy does not take into account connections which may exist between the
various components, as in the important case of sequences arising from iterative methods for solving
a system of linear or nonlinear equations.
2. Scalar sequences
Let (S ) be a scalar sequence converging to a limit S. As explained above, an extrapolation
n
method consists of transforming this sequence into a new one, (T ), by a sequence transformation
n
T :(S ) ! (T ). The transformation T is said to accelerate the convergence of the sequence (S ) if
n n n
and only if
T −S
lim n =0:
n!1 S −S
n
We can then say that (T ) converges (to S) faster than (S ).
n n
The (cid:12)rst methods to have been used were linear transformations
X1
T = a S; n=0;1;:::;
n ni i
i=0
where the numbers a are constants independent of the terms of the sequence (S ). Such a linear
ni n
transformation is usually called a summation process and its properties are completely determined
by the matrix A=(a ). For practical reasons, only a (cid:12)nite number of the coe(cid:14)cients a are di(cid:11)erent
ni ni
from zero for each n. Among such processes are those named after Euler, Cesaro and Ho(cid:127)lder. In
the case of linear methods, the convergence of the sequence (T ) to S for any converging sequence
n
(S ) is governed by the Toeplitz summability theorem; see [115] for a review. Examples of such
n
processes are
1 Xn
T = S
n n+1 i
i=0
or
1 Xn+k
T = S:
n k +1 i
i=n
In the second case, the sequence (T ) also depends on a second index, k, and the convergence has to
n
be studied either when k is (cid:12)xed and n tends to in(cid:12)nity, or when n is (cid:12)xed and k tends to in(cid:12)nity.
C. Brezinski/Journal of Computational and Applied Mathematics 122 (2000) 1{21 3
With respect to convergence acceleration, summation processes are usually only able to accelerate
the convergence of restricted classes of sequences and this is why the numerical analysts of the
20th century turned their e(cid:11)orts to nonlinear transformations. However, there is one exception:
Richardson’s extrapolation process.
2.1. Richardson’s process
It seems that the (cid:12)rst appearance of a particular case of what is now called the Richardson ex-
trapolation process is due to Christian Huygens (1629{1695). In 1903, Robert Moir Milne (1873)
applied the idea of Huygens for computing (cid:25) [101]. The same idea was exploited again by Karl
Kommerell (1871{1948) in his book of 1936 [78]. As explained in [143], Kommerell can be con-
sidered as the real discoverer of Romberg’s method although he used this scheme in the context of
approximating (cid:25).
Let us now come to the procedures used for improving the accuracy of the trapezoidal rule for
computing approximations to a de(cid:12)nite integral. In the case of a su(cid:14)ciently smooth function, the
error of this method is given by the Euler{Maclaurin expansion. In 1742, Colin Maclaurin (1698{
1746) [90] showed that its precision could be improved by forming linear combinations of the
results obtained with various stepsizes. His procedure can be interpreted as a preliminary version of
Romberg’s method; see [49] for a discussion.
In 1900, William Fleetwood Sheppard (1863{1936) used an elimination strategy in the Euler{
Maclaurin quadrature formula with h =r h and 1=r <r <r <(cid:1)(cid:1)(cid:1) to produce a better approxi-
n n 0 1 2
mation to the given integral [132].
In 1910, combining the results obtained with the stepsizes h and 2h, Lewis Fry Richardson (1881{
1953) eliminated the (cid:12)rst term in a discretization process using central di(cid:11)erences [119]. He called
this procedure the deferred approach to the limit or h2-extrapolation. The transformed sequence
(T ) is given by
n
h2 S(h )−h2S(h )
T = n+1 n n n+1 :
n h2 −h2
n+1 n
In a 1927 paper [120] he used the same technique to solve a 6th order di(cid:11)erential eigenvalue problem.
His process was called (h2;h4)-extrapolation. Richardson extrapolation consists of computing the
value at 0, denoted by T(n), of the interpolation polynomial of the degree at most k, which passes
k
through the points (x ;S );:::;(x ;S ). Using the Neville{Aitken scheme for these interpolation
n n n+k n+k
polynomials, we immediately obtain
x T(n)−x T(n+1)
T(n) = n+k+1 k n k
k+1 x −x
n+k+1 n
with T(n)=S .
0 n
Let us mention that Richardson referred to a 1926 paper by Nikolai Nikolaevich Bogolyubov
(born in 1909) and Nikolai Mitrofanovich Krylov (1879{1955) where the procedure (often called
the deferred approach to the limit) can already be found [11].
In 1955, Werner Romberg (born in 1909) was the (cid:12)rst to use repeatedly an elimination approach
for improving the accuracy of the trapezoidal rule [121]. He himself refers to the book of Lothar
Collatz (1910{1990) of 1951 [50]. The procedure became widely known after the rigorous error
4 C. Brezinski/Journal of Computational and Applied Mathematics 122 (2000) 1{21
analysis given in 1961 by Friedrich L. Bauer [3] and the work of Eduard L. Stiefel (1909{1978)
[138]. Romberg’s derivation of his process was heuristic. It was proved by Pierre-Jean Laurent in
1963 [81] that the process comes out from the Richardson process by choosing x =h2 and h =h =2n.
n n n 0
Laurent also gave conditions on the choice of the sequence (x ) in order that the sequences (T(n))
n k
tend to S either when k or n tends to in(cid:12)nity. Weaker conditions were given by Michel Crouzeix
and Alain L. Mignot in [52, pp. 52{55]. As we shall see below, extensions of Romberg’s method
to nonsmooth integrands leads to a method called the E-algorithm.
Applications of extrapolation to the numerical solution of ordinary di(cid:11)erential equations were stud-
ied by H.C. Bolton and H.I. Scoins in 1956 [12], Roland Bulirsch and Josef Stoer in 1964{1966 [47]
and William B. Gragg [65] in 1965. The case of di(cid:11)erence methods for partial di(cid:11)erential equations
was treated by Guri(cid:21)(cid:16) Ivanovich Marchuk and V.V. Shaidurov [91]. Sturm{Liouville problems are
discussed in [117]. Finally, we mention that Heinz Rutishauser (1918{1970) pointed out in 1963
[122] that Romberg’s idea can be applied to any sequence as long as the error has an asymptotic
expansion of a form similar to the Euler{Maclaurin’s.
For a detailed history of the Richardson method, its developments and applications, see [57,77,143].
2.2. Aitken’s process and the (cid:15)-algorithm
The most popular nonlinear acceleration method is certainly Aitken’s (cid:1)2 process which is given
by
S S −S2 (S −S )2
T = n n+2 n+1 = S − n+1 n ; n=0;1;:::
n S −2S +S n S −2S +S
n+2 n+1 n n+2 n+1 n
The method was stated by Alexander Craig Aitken (1895{1967) in 1926 [1], who used it to ac-
celerate the convergence of Bernoulli’s method for computing the dominant zero of a polynomial.
Aitken pointed out that the same method was obtained by Hans von Naegelsbach (1838) in 1876
in his study of Furstenau’s method for solving nonlinear equations [104]. The process was also
given by James Clerk Maxwell (1831{1879) in his Treatise on Electricity and Magnetism of 1873
[95]. However, neither Naegelsbach nor Maxwell used it for the purpose of acceleration. Maxwell
wanted to (cid:12)nd the equilibrium position of a pointer oscillating with an exponentially damped simple
harmonic motion from three experimental measurements. It is surprising that Aitken’s process was
known to Takakazu Seki (1642{1708), often considered the greatest Japanese mathematician. In his
book Katsuyo(cid:22) Sanpo(cid:22), Vol. IV, he used this process to compute the value of (cid:25), the length of a chord
and the volume of a sphere. This book was written around 1680 but only published in 1712 by his
disciple Murahide Araki. Parts of it can be found in [73]. Let us mention that the Japanese characters
corresponding to Takakazu have another pronounciation which is Ko(cid:22)wa. This is the reason why this
mathematician is often called, erroneously as in [29,31] Seki Ko(cid:22)wa.
What makes Aitken’s process so popular is that it accelerates the convergence of all linearly
converging sequences, that is sequences such that 9a 6= 1
S −S
lim n+1 =a:
n!1 S −S
n
It can even accelerate some logarithmic sequences (that is corresponding to a=1) which are those
with the slowest convergence and the most di(cid:14)cult to accelerate.
C. Brezinski/Journal of Computational and Applied Mathematics 122 (2000) 1{21 5
Aitken’s (cid:1)2 process is exact (which means that 8n;T =S) for sequences satisfying, a (S −S)+
n 0 n
a (S −S)=0; 8n; a a 6= 0; a +a 6= 0. Such sequences form the kernel of Aitken’s process.
1 n+1 0 1 0 1
The idea naturally arose of (cid:12)nding a transformation with the kernel
a (S −S)+(cid:1)(cid:1)(cid:1)+a (S −S)=0; 8n;
0 n k n+k
a a 6= 0; a + (cid:1)(cid:1)(cid:1) + a 6= 0. A particular case of k = 2 was already treated by Maxwell in his
0 k 0 k
book of 1873 and a particular case of an arbitrary value of k was studied by T.H. O’Beirne in 1947
[107]. This last work remains almost unknown since it was published only as an internal report. The
problem was handled in full generality by Daniel Shanks (1917{1996) in 1949 [130] and again in
1955 [131]. He obtained the sequence transformation de(cid:12)ned by
(cid:12) (cid:12)
(cid:12) S S (cid:1)(cid:1)(cid:1) S (cid:12)
(cid:12) n n+1 n+k (cid:12)
(cid:12)S S (cid:1)(cid:1)(cid:1) S (cid:12)
(cid:12) n+1 n+2 n+k+1(cid:12)
(cid:12) . . . (cid:12)
(cid:12) .. .. .. (cid:12)
(cid:12) (cid:12)
(cid:12) (cid:12)
S S (cid:1)(cid:1)(cid:1) S
T =e (S )= (cid:12) n+k n+k+1 n+2k (cid:12):
n k n (cid:12) (cid:1)2S (cid:1)(cid:1)(cid:1) (cid:1)2S (cid:12)
(cid:12) n n+k−1 (cid:12)
(cid:12) . . (cid:12)
(cid:12) .. .. (cid:12)
(cid:12) (cid:12)
(cid:12)(cid:1)2S (cid:1)(cid:1)(cid:1) (cid:1)2S (cid:12)
n+k−1 n+2k−2
When k = 1, Shanks transformation reduces to the Aitken’s (cid:1)2 process. It can be proved that
e (S )=S; 8n if and only if (S ) belongs to the kernel of the transformation given above. The same
k n n
ratios of determinants were obtained by R.J. Schmidt in 1941 [127] in his study of a method for
solving systems of linear equations.
The determinants involved in the de(cid:12)nition of e (S ) have a very special structure. They are
k n
called Hankel determinants and were studied by Hermann Hankel (1839{1873) in his thesis in
1861 [72]. Such determinants satisfy a (cid:12)ve-term recurrence relationship. This relation was used
by O’Beirne and Shanks to implement the transformation by computing separately the numerators
and the denominators of the e (S )’s. However, numerical analysts know it is di(cid:14)cult to compute
k n
determinants (too many arithmetical operations are needed and rounding errors due to the computer
often lead to a completely wrong result). A recursive procedure for computing the e (S )’s without
k n
computing the determinants involved in their de(cid:12)nition was needed. This algorithm was obtained in
1956 by Peter Wynn. It is called the (cid:15)-algorithm [147]. It is as follows. One starts with
(cid:15)(n) =0; (cid:15)(n)=S
−1 0 n
and then
1
(cid:15)(n) =(cid:15)(n+1)+ :
k+1 k−1 (cid:15)(n+1)−(cid:15)(n)
k k
Note that the numbers (cid:15)(n)’s (cid:12)ll out a two-dimensional array. The (cid:15)-algorithm is related to Shanks
k
transformation by
(cid:15)(n)=e (S ) and (cid:15)(n) =1=e ((cid:1)S ):
2k k n 2k+1 k n
Thus, the (cid:15)’s with an odd lower index are only auxiliary quantities. They can be eliminated from
the algorithm, thus leading to the so-called cross rule due to Wynn [153].
