# Convergence acceleration of series

## 1  Introduction

Numerous mathematical constants are calculated as limit of series and many of those series are very slow to converge requiring therefore methods to accelerate their convergence. In this section, we recall some definitions and elementary results on well known series.

For any series ĺkak, we denote by sn its partial sums, that is
ę
ę
ę
ę
ę
ę
 s0=a0,
 s1=a0+a1,
 sn=a0+a1+...+an= nĺ k=0 ak.

We are interested in the behavior of sn as n tends to infinite (that is infinite series). The series whose terms have alternative signs are said to be alternating series, we write them as
 sn=a0-a1+...+(-1)nan= nĺ k=0 (-1)kak,
where the (ak) are positive numbers.

Definition 1 An infinite series ĺ ak is said to be convergent if its partial sums (sn) tends to a limit S (called the sum of the series) as n tends to infinite. Otherwise the series is said to be divergent.

The following examples are elementary and occur frequently and it's usual to compare a given series to one of those.

Example 2 (Geometric series) For the geometric series which is defined by ak=xk, the partial sum sn is then given by
 sn=1+x+x2+...+xn= 1-xn+1 1-x ,
hence it converges for | x| < 1 to S=1/(1-x) and diverges otherwise.

Example 3 (Harmonic series) The harmonic series is defined for k > 0 by ak=1/k, and the partial sum sn (sometimes denoted Hn and called Harmonic number) is
 sn=1+ 1 2 + 1 3 ...+ 1 n .
We have for n=2p
 sn
 =
 ćč 1+ 1 2 öř + ćč 1 3 + 1 4 öř + ćč 1 5 +..+ 1 8 öř +...+ ćč 1 2p-1+1 +...+ 1 2p öř
 sn
 >
 1 1 2 2 1 4 4 1 8 +...+2p-1. 1 2p = p 2 ,
therefore the partial sums are unbounded and the harmonic series is divergent (this important result was known to Jakob Bernoulli (1654-1705)  and much before to the french bishop Nicolas Oresme (1323-1382)).

Example 4 (Riemann Zeta series) A natural generalization of the harmonic series is to consider the partial sum
 sn=1+ 1 2s + 1 3s +...+ 1 ns
which converges to the Riemann Zeta function z (s) for any complex numbers s such as Â (s) > 1. The case s=1 is the Harmonic series and the limits for s=2p are well known since Euler who established for example that z(2)=p2/6, z(4)=p4/90,...

## 2  Kummer's acceleration method

Ernst Kummer's (1810-1893) method is very natural to understand and may be used to accelerate the convergence of many series. It goes back to 1837 and the idea is to subtract from a given convergent series ĺak another equivalent series ĺbk whose sum C=ĺk ł 0bk is well known. More precisely suppose that
 lim k®Ą ćč ak bk öř =r ą 0,
then the transformation is given by
 Ąĺ k=0 ak=r Ąĺ k=0 bk+ Ąĺ k=0 ćč 1-r bk ak öř ak=rC+ Ąĺ k=0 ćč 1-r bk ak öř ak.
(1)

The convergence of the latest series is faster because 1-rbk/ak tends to 0 as k tends to infinity .

Example 5 Consider
 S
 =
 Ąĺ k=1 1 k2 ,
 C
 =
 Ąĺ k=1 1 k(k+1) = Ąĺ k=1 ćč 1 k - 1 k+1 öř =1,
we have r = C=1 and Kummer's transformation (1) becomes
 S=1+ Ąĺ k=1 ćč 1- k2 k(k+1) öř 1 k2 =1+ Ąĺ k=1 1 k2(k+1) .
The process can be repeated with this time
 S*
 =
 Ąĺ k=1 1 k2(k+1)
 C
 =
 Ąĺ k=1 1 k(k+1)(k+2) = 1 4 ,
giving
 S= 5 4 +2 Ąĺ k=1 1 k2(k+1)(k+2) .
After N applications of the transformation
 S= Nĺ k=1 1 k2 +N! Ąĺ k=1 1 k2(k+1)(k+2)...(k+N) ,
whose convergence may be rapid enough for suitable values of N (this example is due to Knopp ). Note, that we have used the following family of series to improve the convergence
 C
 =
 Ąĺ k=1 bk= Ąĺ k=1 1 k(k+1)...(k+N) = 1 N.N! ,
 bk
 ~
 1 kN+1 when k® Ą.

## 3  Richardson Extrapolation method

Suppose that the partial sum verifies the relation
 S=sn+ a np +O ćč 1 np+1 öř ,
(2)
where S is the limit of the sum, p a known integer and a a constant which we don't need to determine. Then, following the ideas of Lewis Fry Richardson (1881-1953) , it's natural to define the transformation
 sn(1)= 2ps2n-sn 2p-1 =s2n+ 1 2p-1 (s2n-sn)
in order to eliminate the term in a from relation (2). The process may be repeated on the new series
 sn(2)
 =
 2p+1s2n(1)-sn(1) 2p+1-1 ,
 sn(3)
 =
 2p+2s2n(2)-sn(2) 2p+2-1 ,
 ...

Example 6 Suppose that you have computed, by geometric means and like Archimedes did, the perimeters of the circumscribed polygon with respectively 6,12,24,48 and 96 sides on a unit circle then in this case S=p, p=2 and Richardson's process becomes here
 sn(1)
 =
 s2n+ 1 3 ( s2n-sn) ,
 sn(2)
 =
 s2n(1)+ 1 15 (s2n(1)-sn(1)) ,
 sn(3)
 =
 s2n(2)+ 1 63 (s2n(2)-sn(2)) ,
 ...
It produces the approximations (the number of sides is 6n)
 n
 sn
 sn(1)
 sn(2)
 sn(3)
 1
 3.4641016151...
 3.1324865405...
 3.1416562605...
 3.1415925429...
 2
 3.2153903091...
 3.1410831530...
 3.1415935385...
 3.1415926532...
 4
 3.1596599420...
 3.1415616394...
 3.1415926670...
 8
 3.1460862151...
 3.1415907278...
 16
 3.1427145996...
and finally we find p » s1(4) =3.141592653(637...) which is correct to nearly 10 digits and therefore much more accurate than the 2 digits value estimated by Archimedes !

## 4  Aitken's acceleration and related methods

### 4.1  Aitken's d2-process

In 1926, Alexander Aitken (1895-1967) found a new procedure to accelerate the rate of convergence of a series, the Aitken d2-process . It consists in constructing a new series S(1), whose partial sums sn(1) are given by
 sn(1)=sn+1- (sn+1-sn)2 sn+1-2sn+sn-1 = sn-1sn+1-sn2 sn+1-2sn+sn-1 ,
(3)
where (sn-1,sn,sn+1) are three successive partial sums of the original series. Of course, it's possible to repeat the process on the new series S(1) to obtain S(2)... This process was in fact known to the Japanese mathematician Takakazu Seki Kowa (1642-1708) who used it to compute 10 correct digits of the constant p by accelerating the convergence of the polygon algorithm.

For example for the converging geometric series with 0 < x < 1:
 sn
 =
 1+x+x2+...+xn= 1-xn+1 1-x and the limit is
 S
 =
 lim n® Ą sn= 1 1-x ,
hence
 sn=(1-xn+1)S=S-xn+1S,
and if we replace (sn-1,sn,sn+1) by this expression in (3) we have sn(1)=S ..., the exact limit is given just by three evaluations of the original series. This indicates that for a series almost geometric, Aitken's process may be very efficient.

Example 7 If we consider the extremely slow (logarithmic rate) converging series
 S= lim n® Ą sn= lim n® Ą ćč nĺ k=0 (-1)k k+1 öř =log(2),
The following array illustrates the result of three consecutive applications of Aitken's process to the 7 first computed terms of this series :

 n
 sn
 sn(1)
 sn(2)
 sn(3)
 0
 1.0000000000...
 1
 0.5000000000...
 0.7000000000...
 2
 0.8333333333...
 0.6904761904...
 0.6932773109...
 3
 0.5833333333...
 0.6944444444...
 0.6931057563...
 0.6931488693...
 4
 0.7833333333...
 0.6924242424...
 0.6931633407...
 5
 0.6166666666...
 0.6935897435...
 6
 0.7595238095...
A formal expression may be obtained in this case ; the expressions of (sn-1,sn,sn+1) are
 sn=sn-1+ (-1)n n+1 ,sn+1=sn+ (-1)n+1 n+2
so that
 sn(1)=sn+1+ (-1)n(n+1) (n+2)(2n+3)

For n=1000, it becomes
 s1000
 =
 0.69(264...),
 s1000(1)
 =
 0.693147180(435...).

### 4.2  The e-algorithm

In 1956, Peter Wynn obtained a generalization of Aitken's algorithm which is called the e-algorithm . It is defined by the transformation
 en(k+1)=en+1(k-1)+ 1 en+1(k)-en(k)
(4)

starting with en(-1)=0 and en(0)=sn where sn is the series to be accelerated.

Example 8 Again we take the series
 lim n®Ą ćč nĺ k=0 (-1)k k+1 öř =log(2),
the following array illustrates the result of consecutive applications of Wynn's process (4) to the same as in the previous example 7 first computed terms of the series :

 n
 en(-1)
 en(0)=sn
 en(1)
 en(2)
 en(3)
 en(4)
 en(5)
 0
 0
 1.0000000000...
 -2
 0.7000000000...
 -102
 0.6933333333...
 -3852
 1
 0
 0.5000000000...
 3
 0.6904761904...
 248
 0.6930894308...
 12015
 2
 0
 0.8333333333...
 -4
 0.6944444444...
 -490
 0.6931693989...
 3
 0
 0.5833333333...
 5
 0.6924242424...
 852
 4
 0
 0.7833333333...
 -6
 0.6935897435...
 5
 0
 0.6166666666...
 7
 6
 0
 0.7595238095...
and finally en(6)=0.6931524547... Note that the odd values of the index are intermediate quantities.

## 5  Euler-Maclaurin summation formula

Suppose that the (ak) may be written as the image of a given differentiable function f, that is
 ak=f(k)
and so
 sn= nĺ k=0 ak= nĺ k=0 f(k)
the following famous result due to Leonhard Euler (1732) and Colin Maclaurin (1742, ) states that
 sn= óő n 0 f(t)dt+ 1 2 ( f(0)+f(n)) + mĺ k=1 B2k (2k)! ( f2k-1(n)-f2k-1(0)) +em,n
where em,n is a remainder given by
 em,n= 1 (2m+1)! óő n 0 B2m+1(x- ëx ű )f2m+1(x)dx
where the B2k are Bernoulli's Numbers and Bi(x) being Bernoulli's polynomials (see Bernoulli's number essay for more details or  for a proof).

Example 9 Suppose that s > 1 and
 f(x)= 1 (x+1)s
then
 sn-1
 =
 1+ 12s + 13s +...+ 1ns
 =
 1s-1 ćç č 1- 1ns-1 ö÷ ř + 12 ćç č 1+ 1ns ö÷ ř + B22 ćç č s 1 ö÷ ř ćç č 1- 1ns+1 ö÷ ř +...
 + B2m2m ćç č s+2m-2 2m-1 ö÷ ř ćç č 1- 1ns+2m-1 ö÷ ř +em,n
if n® Ą  in this identity, we find the relation
 z(s) = Ą ĺ k = 1 1ks = 1s-1 + 12 + B22 ćç č s 1 ö÷ ř +...+ B2m2m ćç č s+2m-2 2m-1 ö÷ ř +em,Ą
and by computing z (s)-sn-1 we find the useful algorithm to estimate z (s)
 z(s) = 1+ 12s +...+ 1ns + 1s-1 1ns-1 - 12ns + m ĺ i = 1 B2i2i ćç č s+2i-2 2i-1 ö÷ ř 1ns+2i-1 +( em,Ą-em,n) ,
and for any given integer m, the remainder tends to zero (Exercise : check this with care!).

## 6  Convergence improvement of alternating series

Very efficient methods exist to accelerate the convergence of an alternating series
 S= Ąĺ k=0 (-1)kak,
one of the first is due to Euler. Usually the transformed series converges geometrically with a rate of convergence depending on the method.

Because the speed of converge may be dramatically improved, there is a trick due to Van Wijngaarden which permits to transform any series ĺbk with positive terms bk into an alternating series. It may be written
 Ąĺ k=0 bk= Ąĺ k=0 (-1)kak
with
 ak= Ąĺ j=0 2j b2j(k+1)-1=bk+2b2k+1+4b4k+3+8b8k+7+...
and because 2j(k+1)-1 increases very rapidly with j, only a few terms of this sum are required.

Sometimes a closed form exists for the ak, for example if
 bk= 1 (k+1)2 ,
then
 ak
 =
 Ąĺ j=0 2jb2j(k+1)-1= Ąĺ j=0 2j ( 2j(k+1)) 2 = 1 (k+1)2 Ąĺ j=0 1 2j
 =
 2 (k+1)2 ,
and
 Ąĺ k=0 1 (k+1)2 =2 Ąĺ k=0 (-1)k (k+1)2 .

### 6.1  Euler's method

In 1755, Euler gave a first version of what is now called Euler's transformation. To establish it, observe that the sum of the alternating series may be written as
 2S=a0+(a0-a1)-(a1-a2)+(a2-a3)-...=a0+S1,
with
 S1=D1a0-D1a1+D1a2-...+(-1)kD1ak+...,
and the first difference operator is denoted by
 D1ak=ak-ak+1.

The new series S1 is also alternating and the process may be applied again
 2S1=D1a0+(D1a0-D1a1)-(D1a1-D1a2)+(D1a2-D1a3)-...=D1a0+T1,
with this time
 T1=D2a0-D2a1+D2a2-...+(-1)kD2ak+...,
and the second difference operator is denoted by
 D2ak=D1ak-D1ak+1,
the following theorem is then easily deduced by repeating the process (we omit some details given in ).

Theorem 10 (Euler) Let S= ĺk=0n(-1)kak be a convergent alternating series, then
 S= lim n®Ą sn= lim n®Ą ćč 1 2 nĺ k=0 (-1)k 2k Dka0 öř ,
(5)
and D is the forward difference operator defined by
 Dka0 = Dk-1a0-Dk-1a1 = (-1)k k ĺ j = 0 (-1)j ćç č k j ö÷ ř aj.

The first values of Dka0 are
 D0a0
 =
 a0,
 D1a0
 =
 a0-a1,
 D2a0
 =
 D1a0-D1a1=a0-2a1+a2,
 ...,
 Dka0
 =
 a0-ka1+ k(k-1) 2! a2- k(k-1)(k-2) 3! a3+...+(-1)kak.

Example 11 Sometimes it's possible to find a closed form for the transformation (5), take
 S
 =
 p 4 =1- 1 3 + 1 5 - 1 7 +...= Ąĺ k=0 (-1)kak with
 ak
 =
 1 2k+1 ,
the rate of convergence of this famous series is extremely slow, but Euler's method gives
 D0a0
 =
 1,
 D1a0
 =
 1- 1 3 = 1.2 1.3 ,
 D2a0
 =
 D1a0-D1a1= 1.2 1.3 - 1.2 3.5 = 1.2.4 1.3.5 ,
 D3a0
 =
 D2a0-D2a1= 1.2.4 1.3.5 - 1.2.4 3.5.7 = 1.2.4.6 1.3.5.7 ,
 ...
this suggest the expression (prove it by induction ...)
 Dna0= 1.2.4...(2n) 1.3.5...(2n+1) .
and
 p 4
 =
 1- 1 3 + 1 5 - 1 7 +...= 1 2 ćč 1+ 1.2 1.3 1 2 + 1.2.4 1.3.5 1 22 +... öř
 p 4
 =
 1 2 ćč 1+ 1 3 + 1.2 3.5 + 1.2.3 3.5.7 ... öř = 1 2 Ąĺ k=0 2kk!2 (2k+1)! ,
relation also given by Euler.

### 6.2  Improvement

In 1991, a generalization of Euler's method was given in  (this work in fact generalizes a technique that was used to obtain irrationality measures of some classical constants like log(2) for example). This algorithm is very general, powerful and easy to understand. For an alternating series ĺ(-1)k ak, we assume that there exist a positive function w(x) such that
 ak= óő 1 0 xkw(x)dx.

This relation permits to write the sum of the series as
 Ąĺ k=0 (-1)kak = óő 1 0 ćč Ąĺ k=0 (-1)kxk w(x)dx öř = óő 1 0 w(x) 1+x dx.

Now, for any sequence of polynomials Pn(x) of degree n with Pn(-1) ą 0, we denote by Sn the number
 Sn= 1 Pn(-1) óő 1 0 Pn(-1)-Pn(x) 1+x w(x)dx.

Notice that Sn is a linear combination of the number (ak)0 Ł k < n since if we write Pn(x) = ĺk=0n pk (-x)k, we have easily
 Sn = 1 Pn(-1) n-1ĺ k=0 ck (-1)k ak       with        ck= nĺ j = k+1 pj
(6)

The number Sn satisfies
 Sn = óő 1 0 w(x) 1+x dx - óő 1 0 Pn(x)w(x) Pn(-1)(1+x) dx = S- óő 1 0 Pn(x)w(x) Pn(-1)(1+x) dx.

Therefore we deduce
 | Sn-S| Ł 1 | Pn(-1)| óő 1 0 | Pn(x)| w(x) (1+x) dx Ł Mn |Pn(-1)| | S| ,       Mn= sup x Î [0,1] |pn(x)|

This inequality suggests to choose polynomials with small value of Mn/|Pn(-1)|.

#### 6.2.1  Choice of a family of polynomials

1. A first possible choice is
 Pn(x) = (1-x)n = n ĺ k = 0 ćç č n k ö÷ ř (-x)k,       for which    Pn(-1) = 2n,Mn = 1.

 | Sn-S| Ł | S| 2n

where
 Sn = 12n n-1 ĺ k = 0 (-1)kckak,      ck = n ĺ j = k+1 ćç č n j ö÷ ř .

This choice corresponds in fact to Euler's method.

2. Another choice is
 Pn(x) = (1-2x)n = n ĺ k = 0 2k ćç č n k ö÷ ř (-x)k,    for which    Pn(-1) = 3n,Mn = 1.
 | Sn-S| Ł | S| 3n
where
 Sn = 13n n-1 ĺ k = 0 (-1)kckak,      ck = n ĺ j = k+1 2j ćç č n j ö÷ ř .

3. Chebyshev's polynomials (see ) shifted to [0,1] provide a more efficient acceleration. They satisfy the relation
 Pn(sin2t)=cos(2nt)
(7)
and explicitly writes in the form
 Pn(x) = n ĺ j = 0 nn+j ćç č n+j 2j ö÷ ř 4j(-x)j.
The relation (7) show that Mn=1 and |Pn(-1)| ł [ 1/2](3+Ö8)n > 5.828n/2. This family leads to the following acceleration process :
 | Sn-S| Ł 2| S| 5.828n
where
 Sn = 1Pn(-1) n-1 ĺ k = 0 (-1)kckak,      ck = n ĺ j = k+1 nn+j ćç č n+j 2j ö÷ ř 4j.

4. Other families of orthogonal polynomials such as Legendre's polynomials, Niven's polynomial may give interesting accelerations. More details and results may be found in the very interesting paper .

#### 6.2.2  Examples

Once the choice of a sequence of polynomials is made, it can be applied to compute the value of many alternating series such as
 log(2)= Ąĺ k=0 (-1)k k+1
 ak= 1 k+1 = óő 1 0 xk dx
 p 4 = Ąĺ k=0 (-1)k 2k+1
 ak= 1 2k+1 = óő 1 0 xk x-1/2 2 dx
 (1-2-s)z(s)= Ąĺ k=0 (-1)k (k+1)s ,
 ak= 1 (k+1)s = 1 G(s) óő 1 0 xk | log(x)| s-1dx.
Notice that this latest method is very efficient and may be used to compute the value of the zeta function at values of s with Â(s) > 0 (see ). Another beautiful alternating series whose convergence can be accelerated in this way is
 log(2) ćč g- 1 2 log(2) öř = log(2) 2 - log(3) 3 + log(4) 4 - log(5) 5 +Ľ
where g is the Euler constant (see  for a proof of this formula).

## 7  Acceleration with the Zeta function

If you know how to evaluate the Zeta function at integers values, there is an easy and convenient way to transform your original series in a geometric converging one (based on ). Suppose that you want to estimate a series which has the form
 S=A+ Ąĺ k=2 f ćč 1 k öř
(8)
where A is a constant and where the analytic function f may be written
 f(z)= Ąĺ n=2 anzn
then
 S=A+ Ąĺ k=2 ćč Ąĺ n=2 an kn öř =A+ Ąĺ n=2 an ćč Ąĺ k=2 1 kn öř
hence (8) may be transformed to
 S=A+ Ąĺ n=2 an( z(n)-1) .
(9)

Observe that
 z(n)-1 ~ 1 2n ,
therefore the transformed series (9) has a geometric rate of convergence. This rate may be improved if a few terms of the original series are computed, this time the limit is given by
 S=A+ Mĺ k=2 f ćč 1 k öř + Ąĺ k=M+1 f ćč 1 k öř =A+ Mĺ k=2 f ćč 1 k öř + Ąĺ n=2 an ćč Ąĺ k=M+1 1 kn öř
and again
 S=A+f ćč 1 2 öř +...+f ćč 1 M öř + Ąĺ n=2 an ćč z(n)-1- 1 2n -...- 1 Mn öř
but this time
 z(n)-1- 1 2n -...- 1 Mn =z(n,M+1) ~ 1 (M+1)n
and z(s,a) is the Hurwitz Zeta Function. The rate of convergence remains geometric but this rate may be taken as large as desired by taking a large enough value for M.

### 7.1  Examples

1. The first natural example comes with the (almost) definition series for Euler's constant
 S=g = 1+ Ąĺ k=2 ćč 1 k +log ćč 1- 1 k öř öř
here
 f(z)=z+log(1-z)=- Ąĺ n=2 zn n
and the transformed series is
 g = 1- Ąĺ n=2 ( z(n)-1) n .

2. Another interesting series is Mercator's relation :
 S=log(2)=1- 1 2 + 1 3 - 1 4 +...= 1 2 + Ąĺ k=2 1 (2k-1)2k
and this time
 f(z)= z2 2(2-z) = Ąĺ n=2 zn 2n
giving
 log(2)= 1 2 + Ąĺ n=2 ( z(n)-1) 2n
and if we compute two more terms
 log(2)= 37 60 + Ąĺ n=2 (z(n)-1-1/2n-1/3n) 2n .

3. A check relation

Observe that if we take for the function f ,an=1 for every n
 f(z)= Ąĺ n=2 zn= 1 1-z -1-z= z2 1-z

so that for k > 1
 f ćč 1 k öř = 1 k(k-1) = 1 k-1 - 1 k

and because clearly
 S= Ąĺ k=2 f ćč 1 k öř = Ąĺ k=2 ćč 1 k-1 - 1 k öř =1,

we find the relation
 ĺ n ł 2 ( z(n)-1) = 1.

With the same method come the two other relations
 ĺ n ł 1 ( z(2n)-1)
 =
 3 4
 ĺ n ł 1 ( z(2n+1)-1)
 =
 1 4
which may be used to check the evaluations of the Zeta function on consecutive integer values.

## References


M. Abramowitz and I. Stegun, Handbook of Mathematical Functions, Dover, New York, (1964)


A.C. Aitken, On Bernoulli's numerical solution of algebraic equations, Proc. Roy. Soc. Edinburgh, (1926), vol. 46, p. 289-305


P. Borwein, An efficient algorithm for the Riemann Zeta function, (1995)


H. Cohen, F. Rodriguez Villegas, D. Zagier, Convergence acceleration of alternating series, Bonn, (1991)


P. Flajolet and I. Vardi, Zeta Function Expansions of Classical Constants, (1996)


X. Gourdon and P. Sebah, Numbers, Constants and Computation, World Wide Web site at the adress http://numbers.computation.free.fr/Constants/constants.html, (1999)


R.L. Graham, D.E. Knuth and O. Patashnik, Concrete Mathematics, Addison-Wesley, (1994)


K. Knopp, Theory and application of infinite series, Blackie & Son, London, (1951)


C. Maclaurin, A Treatise of fluxions, Edinburgh, (1742)


L.W. Richardson, The deferred Approach to the Limit, Philosophical Transactions of the Royal Society of London, (1927), serie A, vol. 226


P. Wynn, On a device for computing the em(Sn) transformation, MTAC, (1956), vol. 10, p. 91-96

Back to Numbers, Constants and Computation

File translated from TEX by TTH, version 3.01.
On 8 Jan 2002, 16:57.