Efficient algorithm for summation of some slowly convergent series

Applied Numerical Mathematics 60 (2010) 1442–1453

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

Efficient algorithm for summation of some slowly convergent series Paweł Wozny ´ ∗ Institute of Computer Science, University of Wrocław, ul. Joliot-Curie 15, 50-383 Wrocław, Poland

a r t i c l e

i n f o

Article history: Received 23 November 2009 Received in revised form 29 March 2010 Accepted 1 April 2010 Available online 8 April 2010 Keywords: Convergent acceleration Sequence transformations Orthogonal series

a b s t r a c t The Q transformation, introduced recently by Wozny ´ and Nowak, may serve as a good tool for summation of slowly convergence series. This approach can be easily applied to the case of generalized or basic hypergeometric series, as well as some orthogonal polynomial expansions. It is closely related to the famous Wynn’s epsilon algorithm, Weniger’s or ˇ Homeier’s transformations, and the method introduced by Cížek, Zamastil and Skála. However, it is difficult to use the algorithm proposed by Wozny ´ and Nowak in the general case, because of its high complexity, and some other restrictions. We propose another realization of the Q transformation, which results in obtaining a simpler and faster algorithm. Four illustrative numerical examples are given. © 2010 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction Problems of summation of slowly convergent series and convergent acceleration are of great importance from the practical point of view. One can find them in such areas as computational mathematics, physics or chemistry. So, these problems have been considered by many authors in recent decades (see, e.g., [15,5] and papers cited there). The family of so-called Levin-type methods is one of the most important among the other techniques of summation and convergent acceleration. This type of methods has been introduced by Davin Levin in [12], and can be used to improve the convergence of the wide class of sequences. See also [15,5,10,11,17,2,13]. Recently in [18], a new method has been proposed for the summation of slowly convergent series which can be easily and successfully applied to the case of generalized or basic hypergeometric series, as well as some orthogonal polynomial expansions. The efficiency of the method has been confirmed by some theoretical results, and many numerical experiments. The Q transformation, introduced in the mentioned paper, is closely related to the Wynn’s epsilon algorithm, Weniger’s ˇ [15] or Homeier’s [10] transformations, and the method proposed by Cížek, Zamastil and Skála [6,16]. However, the usage of the method presented in [18] for evaluating the transformation Q is difficult in the case of arbitrary series, because of its high complexity. The main aim of this paper is to propose a new, efficient algorithm of computing the Q transformation which numerical version can be applied to a general case. The paper is organized as follows. In Section 2, we shortly present the method of Wozny ´ and Nowak [18] for the summation of slowly convergent series using the transformation Q. We show also a disadvantage of the method which makes it difficult to apply (see Section 2.2). In Section 3, being the main part of the article, we propose a simpler, and faster algorithm of computing the transformation Q. Results of numerical experiments—done in Maple™ 8, using 32-digit arithmetic—are given in Section 4. In the sequel, we use the following notation. Let E be a shift operator: Exn = xn+1 , more generally Ek xn = xn+k (k ∈ Z), whereas I := E0 —an identity operator. The linear difference operator is the expression of the form

*

Fax: +48 71 3757801. E-mail address: [email protected].

0168-9274/$30.00 © 2010 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2010.04.001

P. Wo´zny / Applied Numerical Mathematics 60 (2010) 1442–1453

K=

v 



λ j (n)E j

1443

 λu (n) ≡ 0, λ v (n) ≡ 0 ,

j =u

while its order is a number ord K := v − u. For example, in the case of the well-known forward difference operator , we have  = E − I. The operations of addition, subtraction, and composition of linear difference operators, as well as their multiplication by scalars are defined as usual. Notice that all difference operators appearing in this paper always act on the variable n, and are linear, i.e., K(can ) = c K(an ), where c is an n-independent quantity. We often omit the words linear and difference, i.e., we write an operator instead of a linear difference operator. Let (c )n (c ∈ C, n ∈ N) denote the Pochhammer symbol defined by

(c )0 := 1,

(c )n := c (c + 1) · · · (c + n − 1) (n = 1, 2, . . .).

We measure the accuracy of approximation

   σ   acc(σ ) := − log10 1 − . s

σ of a nonzero number s by computing the quantity (1.1)

Hence, acc(σ ) is the number of exact significant decimal digits in the approximation

σ of the number s.

2. The transformation Q In this section, we briefly describe a method of summation of slowly convergent series which has been recently proposed in [18]. We also point out a disadvantage of the method which makes it difficult to apply in a general case (see Section 2.2). 2.1. The method Let us consider an infinite series

s :=

∞ 

ak .

(2.1)

k =0

Let sn denote the nth partial sum of the series (2.1), while rn —its nth remainder, i.e.,

sn :=

n −1 

rn :=

ak ,

k =0

∞ 

an+k .

k =0

Suppose we know the linear difference operator L(∞) with the property



L(∞) (rn ) = 0, L(∞) (1) = 0.

Using the relation s = sn + rn , it is easy to justify that

s=

L(∞) (sn ) . L(∞) (1)

In general, it is very difficult, or even impossible, to find the operator L(∞) . So, let us consider the linear operator L(m) being an approximation of the operator L(∞) in the following sense:

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩

L

(m)

 m −1 



an+k

= 0,

k =0

L(m) (1) = 0.

Thus, one can expect that for sufficiently large m,

s≈

L(m) (sn ) . L(m) (1) (m)

In [18], the transformation Q(m) : {sn } → {Qn } defined by (m)

Qn (sn ) :=

L(m) (sn ) L(m) (1)

has been considered. The authors gave a recurrent algorithm of constructing operators L(m) related by the formula

(2.2)

1444

P. Wo´zny / Applied Numerical Mathematics 60 (2010) 1442–1453

L(m) = P(m) L(m−1) (m = 2, 3, . . .)

(2.3)

for some operators P(m) . Algorithm 1. (See [18].) Let L(1) be the linear difference operator such that

L(1) (an ) = 0. Assume that for m  2, we know the operator L(m−1) with the property

L

(m−1)

 m −2 



an+k

= 0.

k =0

The operator L(m) is constructed in the following way. 1. Find the operators P(m) , and R(m) such that

L(m) , P(m) L(m−1) = R(m) where

L(m−1) := Em−1 L(1) E1−m (for details, see [18, §2.1]). 2. Set L(m) := P(m) L(m−1) . From the given construction, it follows that

L

(m)

 m −1 



= 0,

an+k

k =0

as well as that we have

L(m) = P(m) P(m−1) · · · P(1) , where P(1)

(2.4)

:= L(1) . Using the form (2.4) of the operator

L(m) , one can compute the transformation

Q efficiently.

Algorithm 2. (See [18].) (0)

(0)

1. Set N n := sn , D n := 1. (m) (m) 2. Compute quantities N n , and D n in the following recurrent way: (m)

Nn

(m)

3. Then Qn

 (m−1)  , := P(m) Nn

=

(m)

Dn

 (m−1)  (m = 1, 2, . . .). := P(m) D n

(m)

Nn

(m)

Dn

.

As it turned out, the described method of summation of series is closely related to the famous Wynn’s epsilon algorithm, ˇ Weniger’s [15] and Homeier’s [10] transformations, as well as the method by Cížek, Zamastil and Skála [6,16]. For the details ˇ of the mentioned connections, see §2.3 (Wynn’s epsilon algorithm), §3.3 (Weniger’s, and Cížek et al. transformations) and §5 (Homeier’s transformation) in [18]. Notice that this approach seems to be similar to a general construction principle for sequence transformations proposed in [15, §3.2], and developed later by Brezinski and Redivo-Zaglia [3,4], and Brezinski and Matos [2]. The Q transformation can be successfully applied to the case of many types of series, e.g.,

• generalized hypergeometric series (see, e.g., [1, §2.1]),



r +1 Fr

α1 , α2 , . . . , αr +1 β1 , β2 , . . . , βr



x

:=

∞  (α1 )k (α2 )k · · · (αr +1 )k k =0

(β1 )k (β2 )k · · · (βr )k

·

xk k!

,

• basic hypergeometric series (see, e.g., [8, (1.2.22)]),



r +1

φr

α1 , α2 , . . . , αr +1

where (a; q)k :=

β1 , β2 , . . . , βr k−1

j =0 (1 − aq

j



q; x

:=

∞  (α1 ; q)k (α2 ; q)k · · · (αr +1 ; q)k k =0

(β1 ; q)k (β2 ; q)k · · · (βr ; q)k

) denotes the q-Pochhammer symbol,

·

xk

(q; q)k

,

P. Wo´zny / Applied Numerical Mathematics 60 (2010) 1442–1453

1445

• and also some expansions of the type f (x) =

∞ 

ck P k (x),

k =0

where { P k } is a sequence of orthogonal polynomials. 2.2. A disadvantage of the method One can use the Q transformation for summation of slowly convergent series in many cases. The transformation is particularly attractive when the explicit expressions for operators P(m) are known. We have this situation, for example, in the case of the series (see [14, (2), p. 143])

−2 −

π√ 2

2 + 2x =

∞  (−1)k k =0

k2

−

1 4

T k (x)





x ∈ [−1, 1] ,

(2.5)

where T k is the kth Chebyshev polynomial of the first kind. Using the recurrence relation for polynomials T n , we find that

  L(1) = (2n − 1)(2n − 3)E−1 + 2x 4n2 − 1 I + (2n + 1)(2n + 3)E. One can prove (cf. (2.4), [18, Example 5.1]) that

P(m) = (2n − 1)E−1 + 2x(2n + 2m − 1)I + (2n + 4m − 1)E. (m)

Using Algorithm 2, we can easily compute values Qn which—as examined—give good results (see [18, Fig. 1]). Now, let us consider the series of the form (see [14, (8), p. 144])

cos(q arccos x) +

sin(qπ ) qπ

=−

∞ 2q sin qπ  (−1)k

π

k =0

k2 − q2

T k (x)





x ∈ [−1, 1]; q ∈ /Z ,

(2.6)

which for q = 12 reduces to the expansion (2.5). For q = 12 , we are not able to find the explicit forms for operators P(m) . Of course, the use of Q transformation is still possible. For this purpose, we have to derive (symbolically or numerically) the operators P(2) , P(3) , . . . one after another, using Algorithm 1. However, this approach requires solving many systems of linear equations whose orders are increasing with m, hence, it is very expensive. We describe this below in a general situation. Suppose that ord L(1) = r. From the relation (2.3), and Algorithm 1, it follows that

ord P(m) = r ,

ord L(m−1) = ord R(m) = r (m − 1)

(m = 2, 3, . . .).

Assume that for l = 2, 3, . . . , m − 1 the operators L(l) have been determined. So, to find the operator

L(m) , we have to solve symbolically a system of m · r + 1 linear equations. The order of the linear equation, which has to be solved, depends on m, and increases in each step. Consequently, Algorithm 1 is not efficient in a general case. This is the main reason, why we need another, i.e., faster, method of constructing operators L(m) . 3. A new algorithm Below, we give a new, efficient algorithm of constructing operators L(m) which requires solving symbolically a system of linear equations of a fixed order being independent of m. Its numerical realization in the case ord P(1) = 2 is described in details in Section 3.2, while Section 3.3 deals with the general situation, i.e., when ord P(1) = r (r ∈ N). 3.1. The method Let us consider two sequences of operators {P(m) }m , and {Q(m) }m which are constructed recursively using the following algorithm. Algorithm 3. 1. Suppose we know the operator P(1) such that

P(1) (an ) = 0. Set Q(1) := P(1) .

1446

P. Wo´zny / Applied Numerical Mathematics 60 (2010) 1442–1453

2. Solving a system of linear equations (see Section 3.2, and Section 3.3), we find recursively the operators P(m) , and Q(m) (m = 2, 3, . . .) such that

P(m) , P(m) Q(m−1) = Q(m)

(3.1)

P(m) := EP(m−1) E−1 .

(3.2)

where

Let us define

L(m) := P(m) L(m−1) (m = 2, 3, . . .).

L(1) := P(1) , Then

L(m) = P(m) P(m−1) · · · P(1) (m = 1, 2, . . .). Theorem 4. Let L(m) = P(m) P(m−1) · · · P(1) , where operators P(m) have been obtained using Algorithm 3. Then

L

(m)

 m −1 

= 0.

an+k

(3.3)

k =0

Proof. Clearly, the property (3.3) holds for m = 1 (see Step 1 of Algorithm 3). Assume that for any l < m the operator L(l) = P(l) P(l−1) · · · P(1) has the property

L

(l)

 l −1 

an+k

= 0.

(3.4)

k =0

Then we have also

L(l) (an+l−1 ) = 0 (l = 1, 2, . . . , m − 1).

(3.5)

Using assumption (3.4) for l := m − 1, the way of constructing operators P(m) and Q(m) (see (3.1)), definition (3.2) of the operator P(m) , as well as property (3.5) for l := m − 1, we obtain

L

(m)

 m −1 

an+k

(m) (m−1)

=P

L

 m −2 

k =0

an+k + an+m−1

= P(m) L(m−1) (an+m−1 )

k =0

(m)

=P

(3 ) (2 ) (1 )

···P

P

P

(an+m−1 ) = P(m) · · · P(3) P(2) Q(1) (an+m−1 )

= P(m) · · · P(3) Q(2) P(1) (an+m−1 ) = P(m) · · · Q(3) P(2) P(1) (an+m−1 ) = · · · = Q(m) P(m−1) · · · P(1) (an+m−1 ) = Q(m) EP(m−1) · · · P(2) P(1) (an+m−2 ) P(2) = Q(m) EL(m−1) (an+m−2 ) = 0.

2

Thus, the operators P(m) constructed by Algorithm 3 allow us to compute the transformation Q (see (2.2), and Algorithm 2). Most importantly, their calculation is simpler and much less complex than the method presented in Algorithm 1. More precisely, let ord P(1) = r. Then

P(m) = ord Q(m) = r ord P(m) = ord

(m = 1, 2, . . .).

Suppose that the operators P(l) (l = 2, 3, . . . , m − 1) have been constructed. Taking under account that Q(1) = P(1) , one can obtain the operator P(m) by solving the system of 2r linear equations (for details, see Sections 3.2 and 3.3). Recall that in the previous approach (see Section 2.1), the order of the system of linear equations, which has to be solved to obtain the same operator, depends on m and is equal to m · r + 1.

P. Wo´zny / Applied Numerical Mathematics 60 (2010) 1442–1453

1447

3.2. The case ord P(1) = 2 In the case ord P(1) = 1 (cf. Step 1 of Algorithm 3), the presented method does not add anything new. One can show the following relation to the Wynn’s epsilon algorithm: n Qn(m) = ε2m

(see [18, Section 2.3]). So, let P(1) be a second order operator. Without loss of generality, one can assume the following form of the operator ( 1) P :

P(1) := E−1 + p 0 (n)I + p 1 (n)E (1 )

(1 )



(1 )



p 1 (n) ≡ 0 .

Recalling that Q(1) = P(1) , and taking into account the way of constructing operators P(m) , Q(m) (see Step 2 of Algorithm 3), we conclude that also the mentioned operators must be of the form

P(m) = E−1 + p 0 (n)I + p 1 (n)E, (m)

(m)

Q(m) = E−1 + q0 (n)I + q1 (n)E (m = 2, 3, . . .) (m)

(m)

for some functions p i Using the identity

(m)

(m)

, qi

(i = 0, 1) of the variable n.

P(m) , P(m) Q(m−1) = Q(m) P(m) , and definition (3.2) of the operator (m) (m)

P(m) = EP(m−1) E−1 = E−1 + p 0 (n + 1)I + p 1 (n + 1)E,

we obtain



(m−1)

0 = q0

(m)

(m−1)

(n − 1) + p 0 (n) − p 0

 (m) (n) − q0 (n) E−1

 (m−1)  (m) (m−1) (m) (m−1) (m) (m−1) (m) + q1 (n − 1) + p 0 (n)q0 (n) + p 1 (n) − p 1 (n) − q0 (n) p 0 (n + 1) − q1 (n) I  (m)  (m−1) (m) (m−1) (m) (m−1) (m) (m−1) + p 0 (n)q1 (n) + p 1 (n)q0 (n + 1) − q0 (n) p 1 (n + 1) − q1 (n) p 0 (n + 2) E  (m)  (m−1) (m) (m−1) + p 1 (n)q1 (n + 1) − q1 (n) p 1 (n + 2) E2 . (m−1)

Thus, the coefficients of E−1 , I, E, E2 have to be equal to 0. Assuming the functions p i and solving the system of 4 linear equations (cf. (3.6)), which has the structure

⎡ ⎤ ⎡     ⎤ p (m) (n) ⎡⎤ 0 ⎥     ⎢ (m)  ⎥ ⎢ q0 (n) ⎥ ⎢ ⎢ ⎥ ⎥⎢ ⎢  ⎥=⎢⎥ ⎦ ⎢ p (m) (n) ⎥ ⎣ ⎣ ⎦ 0 0 ⎦   ⎣ 1 0 0 0 (m) q1 (n)



(m−1)

, qi

(3.6)

(i = 0, 1) are known,



— nonzero elements ,

(m)

(m)

(m)

(m)

one can find the expressions for the coefficients of the operators P(m) , Q(m) , i.e., functions p 0 , p 1 , and q0 , q1 . Admittedly, the system of linear equations, which has to be solved, is simple. However, we have to remember that the (m−1) (m−1) coefficients appearing in this linear system are functions of the variable n (p i , qi (i = 0, 1)) which may be very complicated for larger values of m. Therefore, it is worth to consider a numerical version of the sketched procedure. (1) (1) (1) (1) Assume we know the values p 0 (n), p 1 (n), q0 (n), q1 (n) for n = 2, 3, . . . , N − 1, where 3| N. Proceeding as described (m)

(m)

earlier, i.e., solving the appropriate systems of linear systems (for fixed n, m), we can compute the values p 0 (n), p 1 (n), (m)

(m)

(m)

q0 (n), q1 (n) for m = 2, 3, . . . , N /3, and n = m + 1, m + 2, . . . , N − 2m + 1. It allows us to calculate the values Qn for m = 2, 3, . . . , N /3, and n = m + 1, m + 2, . . . , N − 2m + 1. So, for example, if N = 12, one can compute the following table:

1448

P. Wo´zny / Applied Numerical Mathematics 60 (2010) 1442–1453

s1

= Q1(0)

s2

= Q2(0)

Q2(1)

s3

= Q3(0)

Q3(1)

s4

= Q4

( 0)

(1 )

Q4

( 0)

Q5

( 0)

Q6

( 0)

Q7

( 0)

Q8

s5

= Q5

s6

= Q6

s7

= Q7

s8

= Q8

Q3(2) (2 )

Q4

(2 )

Q5

(2 )

Q6

(2 )

Q7

Q4

(1 )

Q5

(1 )

Q6

(1 )

Q7

(1 )

Q8

(1 )

Q9

( 0)

Q9

s9

= Q9

s10

= Q10

( 0)

Q10

s11

= Q11

( 0)

Q11

s12

( 0)

(3 ) (3 )

(4 )

Q5

(3 ) (3 )

(2 ) (2 )

(1 ) (1 )

= Q12

∞

Remark 5. Note that the presented numerical algorithm can be applied to the arbitrary series k=0 ak . Indeed, let us assume we know the numerical values of the terms a0 , a1 , . . . , a N . Now, all we need to use the method is the operator P(1) . If we do not have any additional information about the sequence {ak } then we can define the operator P(1) , for example, by

P(1) = E−1 + p 0 (n)I + p 1 (n)E, (1 )

(1 )

where (1 )

p 0 (n) := −2

an−1 an

an−1

(1 )

p 1 (n) :=

,

(n = 1, 2, . . . , N − 1).

an+1

Then we have

P(1) (an ) = 0 for n = 1, 2, . . . , N − 1, (m)

what allows us to compute the transformation Qn

for m = 2, 3, . . . , N /3, and n = m + 1, m + 2, . . . , N − 2m + 1.

3.3. The case ord P(1) = r Now, we briefly describe the general situation. Let ord P(1) = r, where r ∈ N. Without loss of generality, we can assume that the operator P(1) has the following form:

P(1) =

v 

(1 )

p j (n)E j

(u < v ; v − u = r )

j =u

(1)

for some functions p j

(1)

(1)

( j = u , u + 1, . . . , v ) with p u (n) ≡ 1, and p v (n) ≡ 0. From Algorithm 3 of constructing operators

P(m) , and Q(m) , it follows that these operators are of the form

P(m) =

v 

(m)

p j (n)E j ,

j =u

Q(m) =

v 

(m)

q j (n)E j

(m  1),

j =u

where (m)

p u (n) ≡ 1,

(m)

qu (n) ≡ 1 (m = 1, 2, . . .). (m−1)

Suppose that operators P(m−1) , Q(m−1) have been constructed, i.e., that the coefficients p j (m)

(m)

1, . . . , v ) are known. The coefficients p j , and q j 2r linear equations implying by the relation

P(m) P(m) Q(m−1) = Q(m) (cf. Step 2 of Algorithm 3, and Section 3.2).

of the operators P(m) ,

(m−1)

, and q j

( j = u, u +

Q(m) can be computed by solving the system of

P. Wo´zny / Applied Numerical Mathematics 60 (2010) 1442–1453

1449

(11)

Fig. 1. The quantities acc(s33 ) (++), and acc(Q12 ) (− ◦− −) ◦ for q = 3.33, and x = −1, −0.96, . . . , 1.

4. Examples The numerical version of the algorithm of constructing operators L(m) , presented in Section 3, has been implemented in the computer algebra system Maple™ by R. Nowak (the packages are available upon request; please write to the following e-mail address: [email protected]). Results of the experiments given below have been obtained using Maple™ 8 and 32-digit arithmetic. 4.1. First example Let us back to the expansion (2.6),

cos(q arccos x) +

sin(qπ ) qπ

=−

∞ 2q sin qπ  (−1)k

π

k =0

k2 − q2

T k (x)





x ∈ [−1, 1]; q ∈ /Z .

Using the well-known recurrence relation for Chebyshev polynomials,

T n (x) = 2xT n−1 (x) − T n−2 (x), we find the operator P(1) :

      P(1) = (n − 1)2 − q2 E−1 + 2x n2 − q2 I + (n + 1)2 − q2 E with the property

 (1 )

P

 (−1)n T n (x) = 0. n2 − q 2

Set q = 3.33. The partial sum s33 of the series (2.6) (x ∈ [−1, 1]) gives no more than 5 exact significant decimal digits. We (11) obtain better results computing the transformation Q12 being the linear combination of the partial sums s1 , s2 , . . . , s33 . In Fig. 1, the results of numerical tests made for x = −1, −0.96, . . . , 1 are presented. 4.2. Second example Now, let us consider the following expansion in terms of Jacobi polynomials [7, §10.20, (3)]:

1450

P. Wo´zny / Applied Numerical Mathematics 60 (2010) 1442–1453

(11)

Fig. 2. The quantities acc(s33 ) (++), and acc(Q12 ) (− ◦− −) ◦ for

ρ = −3.33, α = 5.3, β = 1.3, and x = −0.96, −0.92, . . . , 1.

∞

(1 − x)ρ = 2ρ

Γ (α + ρ + 1)Γ (α + β + 1)  (α ,β) ck P k (x), Γ (α + 1)Γ (α + β + ρ + 2)

(4.1)

k =0

where

ck :=

(2k + α + β + 1)(−ρ )k (α + β + 1)k , (α + 1)k (α + β + ρ + 2)k

and x ∈ (−1, 1),

α , β > −1, −ρ < min(α + 1, α2 + 34 ).

The operator P(1) has to be chosen so that

  (α ,β) P(1) cn P n (x) = 0. In this case, we have

P(1) =

2(n + β)(2n + α + β + 2) −1 E 2n + α + β − 1

 (n + α + β + ρ + 1)  I − (2n + α + β)(2n + α + β + 2)x + α 2 − β 2 (n + α + β)(n − ρ − 1) 2(n + 1)(2n + α + β)(n + α + 1)(n + α + β + ρ + 1)2 + E, (2n + α + β + 3)(n − ρ − 1)2 (n + α + β)

what follows from the recurrence relation for Jacobi polynomials (see, e.g., [1, Section 2, Exercise 22]). (11) In Fig. 2, the accuracies given by the partial sum s33 of the series (4.1), and the transformation Q12 for α = 5.3, β = 1.3, and x = −0.96, −0.92, . . . , 1 are shown. 4.3. Third example In the case of the series

0.0185355891429332393 . . . =

∞  k =0

(−1)k , √ (k + 33) k + 1

one can define

P(1) =

√

√

n(n + 32)E−1 + 2 n + 1(n + 33)I +

√

n + 2(n + 34)E.

ρ = −3.33,

P. Wo´zny / Applied Numerical Mathematics 60 (2010) 1442–1453

1451

(11)

Fig. 3. The quantities acc(s33 ) (++), and acc(Q12 ) (− ◦− −) ◦ for xk (k = 1, 2, . . . , 20).

Then

 P(1)

(−1)n √ (n + 33) n + 1

 = 0.

The considered series is slowly convergent. We have, for example,

acc(s100 ) = 1.7,

acc(s500 ) = 2.65,

acc(s1000 ) = 3.08,

acc(s10000 ) = 4.57

(11)

(see (1.1)), while the transformation Q12 , which can be computed using the partial sums s1 , s2 , . . . , s33 , gives more than (11)

28 exact significant decimal digits. More precisely, acc(Q12 ) = 28.65. 4.4. Fourth example Let the following expansion be given [9, Eq. (1.516.1)]:

1 2

2

ln (1 + x)

=

 ∞  (−1)k+1 xk+1 k =1

k+1

1+

1 2

+ ··· +

1



k





x2 < 1 .

Setting

P(1) = x2n2 E−1 + x(n + 1)(2n + 1)I + (n + 1)(n + 2)E, one can check that

 P(1)

  (−1)n+1 xn+1 1 1 1 + + ··· + = 0. n+1 2 n

Let us define

xk = −1 +

2k 101

(k = 1, 2, . . . , 100). (11)

In Fig. 3, the accuracies given by the partial sum s33 , and the transformation Q12 In Fig. 4, one can see the same quantities for xk and k = 80, 81, . . . , 100.

for xk (k = 1, 2, . . . , 20) are presented.

1452

P. Wo´zny / Applied Numerical Mathematics 60 (2010) 1442–1453

(11)

Fig. 4. The quantities acc(s33 ) (++), and acc(Q12 ) (− ◦− −) ◦ for xk (k = 80, 81, . . . , 100).

5. Conclusions ´ and Nowak in [18]. Recall that the We propose a new realization of the Q transformation introduced recently by Wozny algorithm presented in the mentioned paper requires solving of the systems of linear equations whose orders are increasing in each step, what can be seen as a disadvantage of this approach. In this article, we propose the faster and more efficient method. It requires solving of the systems of linear equations of the same order in each step of the new algorithm. Thus, the Q transformation can be successfully applied to the wide class of slowly convergent series. Acknowledgements The author would like to thank to Prof. E.J. Weniger, and to Prof. S. Lewanowicz for their valuable comments and suggestions concerning this research. He also wish to thank to R. Nowak for the implementation of the presented method in the computer algebra system Maple™. Finally, the author would like to express his gratitude to the referees for the remarks and suggestions which were helpful in preparing the revised version of the manuscript. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]

G.E. Andrews, R. Askey, R. Roy, Special Functions, Cambridge University Press, Cambridge, 1999. C. Brezinski, A.C. Matos, A derivation of extrapolation algorithms based on error estimates, J. Comput. Appl. Math. 66 (1996) 5–26. C. Brezinski, M. Redivo-Zaglia, A general extrapolation procedure revisited, Adv. Comput. Math. 2 (1994) 461–477. C. Brezinski, M. Redivo-Zaglia, On the kernel of sequence transformations, Appl. Numer. Math. 16 (1994) 239–244. C. Brezinski, M.R. Zaglia, Extrapolation Methods: Theory and Practice, North-Holland, Amsterdam, 1991. ˇ J. Cížek, J. Zamastil, L. Skála, New summation technique for rapidly divergent perturbation series. Hydrogen atom in magnetic field, J. Math. Phys. 44 (2003) 962–968. A. Erdélyi, W. Magnus, F. Oberhettinger, F.G. Tricomi, Higher Transcendental Functions, vol. II, McGraw–Hill Book Company, Inc., New York, 1953, based, in part, on notes left by Harry Bateman. G. Gasper, M. Rahman, Basic Hypergeometric Series, 2nd ed., Encyclopedia of Mathematics and its Applications, vol. 96, Cambridge University Press, Cambridge, 2004. I.S. Gradshteyn, I.M. Ryzhik, Table of Integrals, Series and Products, 7th ed., Academic Press Inc., New York, 2007. H.H. Homeier, A Levin-type algorithm for accelerating the convergence of Fourier series, Numer. Algorithms 3 (1992) 245–254. H.H. Homeier, Scalar Levin-type sequence transformations, J. Comput. Appl. Math. 122 (2000) 81–147. D. Levin, Development of non-linear transformations of improving convergence of sequences, Int. J. Comput. Math. B 3 (1973) 371–388. A.C. Matos, Linear difference operators and acceleration methods, IMA J. Numer. Anal. 20 (2000) 359–388. S. Paszkowski, Numerical Applications of Chebyshev Polynomials and Series (Zastosowania Numeryczne Wielomianów i Szeregów Czebyszewa), PWN, Warszawa, 1975. E.J. Weniger, Nonlinear sequence transformations for the acceleration of convergence and the summation of divergent series, Comput. Phys. Rep. 10 (1989) 189–371.

P. Wo´zny / Applied Numerical Mathematics 60 (2010) 1442–1453

1453

ˇ [16] E.J. Weniger, Mathematical properties of a new Levin-type sequence transformation introduced by Cížek, Zamastil, and Skála. I. Algebraic theory, J. Math. Phys. 45 (2004) 1209–1246. [17] J. Wimp, Sequence Transformations and Their Applications, Mathematics in Science and Engineering, vol. 154, Academic Press, New York, 1981. [18] P. Wozny, ´ R. Nowak, Method of summation of some slowly convergent series, Appl. Math. Comput. 215 (2009) 1622–1645.