Applied Mathematics and Computation 215 (2009) 1622–1645
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Method of summation of some slowly convergent series Paweł Woz´ny, Rafał Nowak * 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
Keywords: Convergence acceleration Extrapolation Levin-type method Nonlinear sequence transformation Iterative methods Power series Orthogonal series Hypergeometric series Basic hypergeometric series
a b s t r a c t A new method of summation of slowly convergent series is proposed. It may be successfully applied to the summation of generalized and basic hypergeometric series, as well as some classical orthogonal polynomial series expansions. In some special cases, our algorithm is equivalent to Wynn’s epsilon algorithm, Weniger transformation [E.J. Weniger, Nonlinear sequence transformations for the acceleration of convergence and the summation of divergent series, Computer Physics Reports 10 (1989) 189–371] or the technique recently introˇ ízˇek et al. [J. C ˇ ízˇek, J. Zamastil, L. Skála, New summation technique for rapidly duced by C divergent perturbation series. Hydrogen atom in magnetic field, Journal of Mathematical Physics 44 (3) (2003) 962–968]. In the case of trigonometric series, our method is very similar to the Homeier’s H transformation, while in the case of orthogonal series — to the K transformation. Two iterated methods related to the proposed method are considered. Some theoretical results and several illustrative numerical examples are given. Ó 2009 Elsevier Inc. All rights reserved.
1. Introduction Many fields of numerical mathematics come up against the convergence problem. In the case of slow convergence, many methods of acceleration were proposed. But still the development of the new algorithms plays an important role. In this paper, a new method of summation of slowly convergent series is introduced, and its efficiency is examined by many numerical examples. The method is a Levin-type sequence transformation; its relations to the classic methods are analysed. For the purpose of this paper, let us consider an infinite series
s¼
1 X
an
n¼0
Pn1 P ð1Þ with terms an , partial sums sn ¼ j¼0 aj , and remainders rn ¼ 1 annihilate the j¼n aj . Let the linear difference operator L ð1Þ remainder r n , i.e. L ðrn Þ ¼ 0. Throughout the paper, we assume that all difference operators act on the variable n. Observe that
s¼
Lð1Þ ðsn Þ Lð1Þ ð1Þ
provided Lð1Þ ð1Þ–0. In general, it is difficult to obtain explicit form of the annihilator Lð1Þ . Let the linear difference operator LðmÞ ; m 2 N, be an approximation of Lð1Þ in the sense that
(
ðmÞ
LðmÞ ðr n Þ ¼ 0; LðmÞ ð1Þ – 0;
* Corresponding author. E-mail addresses:
[email protected] (P. Woz´ny),
[email protected] (R. Nowak). 0096-3003/$ - see front matter Ó 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2009.07.016
P. Woz´ny, R. Nowak / Applied Mathematics and Computation 215 (2009) 1622–1645
1623
ðmÞ
where r n :¼ rn r nþm . Then, we can expect that LðmÞ ðsn Þ=LðmÞ ð1Þ gives an approximation of s, of accuracy growing when m ! 1. In this paper, we consider the Q transformation, where QðmÞ : fsn g#fQðmÞ n g; m 2 N, is defined by
QðmÞ :¼ n
LðmÞ ðsn Þ : LðmÞ ð1Þ
Since LðmÞ is a linear difference operator, Q is a Levin-type sequence transformation; see, e.g., [1], [3, Section 2.7], [4] or [5]. Different Levin-type transformations, and hence different methods of convergence acceleration, were studied by Weniger [1], Homeier [4], Brezinski and Matos [6], and Matos [7]. In Section 2, we give a method of obtaining the annihilators LðmÞ , satisfying LðmÞ ¼ PðmÞ Lðm1Þ for some linear difference operators PðmÞ ; m 2 N. Thus we will consider
LðmÞ ¼ PðmÞ Pðm1Þ Pð1Þ and we will say that Q is determined by operators PðmÞ . Next, we will make use of the fact that numerators and denominators on the right side of (2.3) satisfy the recurrence relation
X nfðmÞg ¼ fPgfðmÞg ðX fðm—1Þg Þ n ðmÞ to propose an efficient algorithm of computing QðmÞ ; m 2 N, provided n . One could obtain explicit forms of the operators P the analytical form of terms an is known. In the sequel, we give some applications of the Q transformation. In Section 3, we give explicit forms of the operators PðmÞ in the case of generalized hypergeometric series, while in Section 4 we give analogous results for basic hypergeometric series. In Section 5, we show that our method may be successfully applied to the summation of some classical orthogonal polynomial series, too. Two iterated methods, related to the Q transformation, are considered in Section 6. In the sequel, we will use the following notation:
shift operator: Exn ¼ xnþ1 ; more generally Ek xn ¼ xnþk ; k 2 Z; identity operator: I :¼ E0 ; linear difference operator:
K¼
v X
kj ðnÞEj
ðku ðnÞX0; kv ðnÞX0Þ of order ord K :¼ v u;
ð1:1Þ
j¼u
forward difference operator: D :¼ E I. All computations and numerical tests are made in Maplee11 system, using 128-digit arithmetic. Following Paszkowski [8], we measure the accuracy of approximation r by
r accðrÞ :¼ log10 1; s i.e., by the number of exact significant decimal digits in approximation
r of the sum s.
2. The method Consider the series
s¼
1 X
an
ð2:1Þ
n¼0
Pn1 P with terms an , partial sums sn ¼ j¼0 aj , and remainders rn ¼ 1 j¼n aj . In Section 2.1, we propose a method of obtaining the ðmÞ annihilators L ; m 2 N, satisfying
LðmÞ ðrðmÞ n Þ ¼ 0;
ð2:2Þ
where
rðmÞ ¼ rn r nþm ¼ an þ anþ1 þ þ anþm1 : n We consider the transformation QðmÞ : fsn g#fQðmÞ n g, defined by
QðmÞ :¼ n
LðmÞ ðsn Þ ; LðmÞ ð1Þ
m 2 N:
In Section 2.2, we describe an algorithm of computing QðmÞ n .
ð2:3Þ
P. Woz´ny, R. Nowak / Applied Mathematics and Computation 215 (2009) 1622–1645
1624
2.1. Constructing the annihilators LðmÞ According to (2.2) with m ¼ 1, the operator Lð1Þ should annihilate the term an . A possible choice is the first-order operator
1 1 1 Lð1Þ :¼ D I ¼ E I: an anþ1 an
ð2:4Þ
The operators Lð2Þ ; Lð3Þ ; . . . are constructed recursively by
LðmÞ ¼ PðmÞ Lðm1Þ ;
m P 2:
ð2:5Þ ðmÞ
ðm1Þ
We give the details of a typical step of this construction. Since r n ¼ r n þ anþm1 and by the fact that Lðm1Þ annihilates the ðm1Þ , we are looking for a linear difference operator PðmÞ satisfying term rn
LðmÞ ðanþm1 Þ ¼ PðmÞ Lðm1Þ ðanþm1 Þ ¼ 0:
ð2:6Þ
Let us put
e L ðm1Þ :¼ Em1 Lð1Þ Emþ1
ð2:7Þ
and observe that e L ðm1Þ annihilates the term anþm1 . Now, assume that operators P
ðmÞ
ðmÞ
and R
are such that
L ðm1Þ : PðmÞ Lðm1Þ ¼ RðmÞ e
ð2:8Þ
We obtain ðmÞ ðm1Þ LðmÞ ðr ðmÞ ðrn Þ þ LðmÞ ðanþm1 Þ ¼ PðmÞ Lðm1Þ ðr nðm1Þ Þ þ RðmÞ e L ðm1Þ ðanþm1 Þ ¼ PðmÞ ð0Þ þ RðmÞ ð0Þ ¼ 0; n Þ ¼ L
and thus (2.2) holds. Note that, in view of (2.5), the operators LðmÞ are of the form (we assume Pð1Þ :¼ Lð1Þ )
LðmÞ ¼ PðmÞ Pðm1Þ Pð1Þ :
ð2:9Þ
Now we explain how to obtain operators P ator Lð1Þ is of the form
Lð1Þ ¼
‘ X
ð1Þ
ð1Þ
ki ðnÞEi ;
ðmÞ
ðmÞ
;R
satisfying (2.8). Without loss of generality we can assume that oper-
ð1Þ
k0 ðnÞX0; k‘ ðnÞX0:
ð2:10Þ
i¼0
Thus operator (2.7) can be written as follows:
e L ðm1Þ ¼
‘ X
ð1Þ
ki ðn þ m 1ÞEi ;
m P 2:
ð2:11Þ
i¼0
Assume that Lðm1Þ has the form
Lðm1Þ ¼
ðm1Þ‘ X
ðm1Þ
kj
ðnÞEj ;
m P 2:
ð2:12Þ
j¼0
Let us suppose that condition (2.8) can be satisfied by the operators of the form
PðmÞ ¼
‘ X
i pðmÞ RðmÞ ¼ i ðnÞE ;
ðm1Þ‘ X
i¼0
j qðmÞ j ðnÞE ;
ð2:13Þ
j¼0
where fpi ðnÞg‘i¼0 ; fqj ðnÞgm‘‘ j¼0 are unknown functions. By performing multiplication of the operators in (2.8), we obtain ðmÞ
‘ m‘‘ X X i¼0
ðmÞ
ðm1Þ
kj
ðmÞ
ðn þ iÞ pi ðnÞEiþj ¼
m‘‘ X
‘ X
j¼0
i¼0
j¼0
ð1Þ
ðmÞ
ki ðn þ m 1 þ jÞ qj ðnÞEiþj :
ð2:14Þ
Equating the variable coefficients of Ek on both sides of the equation, we obtain the following system of m‘ þ 1 equations:
X
ðm1Þ
kj
ðmÞ
ðn þ iÞ pi ðnÞ ¼
iþj¼k
X
ð1Þ
ðmÞ
ki ðn þ m 1 þ jÞ qj ðnÞ;
which are satisfied by a set of m‘ þ 2 unknown functions corresponding to k ¼ 0; 1; 2, are as follows: ðm1Þ
k ¼ 0; 1; . . . ; m‘
ðmÞ pðmÞ i ðnÞ and qj ðnÞ. The beginning three equations of the system,
ðmÞ ð1Þ ðmÞ ðnÞ 0 ðnÞ ¼ k0 ðn þ m 1Þ 0 ðnÞ; ðm1Þ ðmÞ ðm1Þ ðmÞ ð1Þ ðmÞ ð1Þ ðn þ 1Þ 1 ðnÞ þ k1 ðnÞ 0 ðnÞ ¼ k0 ðn þ mÞ 1 ðnÞ þ k1 ðn þ m k0 ðm1Þ ðmÞ ðm1Þ ðmÞ ðm1Þ ðmÞ ðn þ 2Þ 2 ðnÞ þ k1 ðn þ 1Þ 1 ðnÞ þ k2 ðnÞ 0 ðnÞ k0 ð1Þ ðmÞ ð1Þ ðmÞ ð1Þ ðmÞ ¼ k0 ðn þ m þ 1Þ 2 ðnÞ þ k1 ðn þ mÞ 1 ðnÞ þ k2 ðn þ m 1Þ 0 ðnÞ:
k0
ð2:15Þ
iþj¼k
p
q
p p
p
q
p
q
p
q
q
ðmÞ
1Þ q0 ðnÞ;
P. Woz´ny, R. Nowak / Applied Mathematics and Computation 215 (2009) 1622–1645
1625
Since the one unknown is superabundant, so it can be chosen arbitrarily. Then the next unknowns are determined uniquely. Observe that such calculations can be made using computer algebra system like Maplee or Mathematicaâ . Observe that left-hand side of (2.14) is actually an explicit form of LðmÞ (cf. (2.5)), which can be rearranged in the form (2.12) with m 1 replaced by m. Notice that such construction of operators LðmÞ depends on the choice of the first operator Lð1Þ , as well as of the operators ðmÞ and RðmÞ , satisfying (2.8) (cf. (2.13)). In most cases, solving (2.15) gives coefficients of the linear difference operators P PðmÞ ; RðmÞ of the form (2.13). Also, we have
ord PðmÞ ¼ ord e L ðm1Þ ¼ ord Lð1Þ ¼ ‘;
and ord RðmÞ ¼ ðm 1Þ‘:
ð2:16Þ
As we shall see later, the complexity and efficiency of our algorithm depend on the order of Lð1Þ , thus — on the orders of PðmÞ ; m 2 N. Choosing higher order operator Lð1Þ may result in a better acceleration of the convergence of series. Since the operators LðmÞ are of the form (2.9), we have
ord LðmÞ ¼ m‘:
ð2:17Þ ð1Þ
In Section 2.3, we will show that if ord L ¼ 1, then our algorithm is equivalent to Wynn’s epsilon algorithm. In the case of generalized hypergeometric series, we propose to use higher order operator Lð1Þ (see Section 3). Our method can then be ˇ ízˇek et al. (see [2], or [9]). In the case of basic identified as a special case of the G transformation recently proposed by C hypergeometric series, considered in Section 4, our transformation seems to be a new method of summation. In the case of orthogonal polynomial series with terms an ¼ bn W n ðxÞ, we propose to annihilate the term an by using the second-order linear recurrence equation, satisfied by polynomials W n ðxÞ (see Section 5 and (6.8)). Such idea is similar to the one proposed by Homeier in derivation of H; I, or K transformation; see [4,10–12]. Notice that annihilators LðmÞ of the form (2.9) were studied, among others, by Brezinski and Redivo-Zaglia [13], or Matos [7], in derivation of some Levin-type sequence transformations. Example 2.1. Consider the series 1 X
2 ln 2 ¼
an
with an :¼
n¼0
1 : ðn þ 1Þ 2n
ð2:18Þ
One can check that the operator
Lð1Þ :¼ ðn þ 1ÞI ð2n þ 4ÞE
ð2:19Þ
satisfies Lð1Þ ðan Þ ¼ 0. We show in details the construction of operators PðmÞ for m ¼ 2; 3. We get
e L ð1Þ ¼ ðn þ 2ÞI ð2n þ 6ÞE;
e L ð2Þ ¼ ðn þ 3ÞI ð2n þ 8ÞE;
respectively (cf. (2.11)). Since Lð1Þ has the form (2.10) with ‘ ¼ 1, the assumption (2.8) (for m ¼ 2) can be satisfied by operators ð2Þ
ð2Þ
ð2Þ
Pð2Þ ¼ p0 ðnÞI þ p1 ðnÞE;
ð2Þ
Rð2Þ ¼ q0 ðnÞI þ q1 ðnÞE ð2Þ
(cf. (2.13)). Remark that one of the unknowns can be chosen arbitrarily, so we put q0 ðnÞ ¼ n þ 1, and obtain ð2Þ
ð2Þ
ð2Þ
ð2Þ
Pð2Þ Lð1Þ ¼ ðn þ 1Þ p0 ðnÞ I þ ½ðn þ 2Þ p1 ðnÞ ð2n þ 4Þ p0 ðnÞ E ð2n þ 6Þ p1 ðnÞ E2 ; ð2Þ e ð1Þ
ð2Þ 1 ðnÞ
¼ ðn þ 2Þðn þ 1ÞI þ ½ðn þ 3Þ q
R L
ð2Þ 2 1 ðnÞ E :
ð2n þ 6Þðn þ 1ÞE ð2n þ 8Þ q
Equating the coefficients of I; E; E2 , we compute
pð2Þ pð2Þ qð2Þ 0 ðnÞ ¼ n þ 2; 1 ðnÞ ¼ ð2n þ 8Þ; 1 ðnÞ ¼ ð2n þ 6Þ: Hence
Pð2Þ ¼ ðn þ 2ÞI ð2n þ 8ÞE; and L
ð2Þ
ð2Þ ð1Þ
¼P L ð3Þ
P
ð2:20Þ 2
¼ ðn þ 1Þðn þ 2ÞI 4 ðn þ 2Þðn þ 3ÞE þ 4 ðn þ 3Þðn þ 4ÞE . Now we will find the operators
ð3Þ 0 ðnÞI
¼p
Rð2Þ ¼ ðn þ 1ÞI ð2n þ 6ÞE;
ð3Þ
ð3Þ
ð3Þ
ð3Þ
Rð3Þ ¼ q0 ðnÞI þ q1 ðnÞE þ q2 ðnÞE2 ;
þ p1 ðnÞE;
ð3Þ
satisfying (2.8) for m ¼ 3. We put q0 ðnÞ ¼ ðn þ 1Þðn þ 2Þ and obtain ð3Þ
ð3Þ
ð3Þ
Pð3Þ Lð2Þ ¼ ðn þ 1Þðn þ 2Þ p0 ðnÞI ðn þ 2Þðn þ 3Þð4 p0 ðnÞ p1 ðnÞÞE ð3Þ ð3Þ ð3Þ 4ðn þ 3Þðn þ 4Þ p1 ðnÞ p0 ðnÞ E2 þ 4ðn þ 4Þðn þ 5Þ p1 ðnÞE3 ;
ð3Þ ð3Þ ð3Þ ð3Þ L ð2Þ ¼ ðn þ 1Þðn þ 2Þðn þ 3ÞI ðn þ 4Þð2ðn þ 1Þðn þ 2Þ q1 ðnÞÞE þ ðn þ 5Þ q2 ðnÞ 2 q1 ðnÞ E2 2ðn þ 6Þ q2 ðnÞE3 : Rð3Þ e
P. Woz´ny, R. Nowak / Applied Mathematics and Computation 215 (2009) 1622–1645
1626
Using computer algebra system like Maplee or Mathematicaâ , one can compute
pð3Þ pð3Þ 0 ðnÞ ¼ n þ 3; 1 ¼ ð2n þ 12Þ; ð3Þ q1 ðnÞ ¼ ð2n þ 4Þð2n þ 8Þ; qð3Þ 2 ðnÞ ¼ ð2n þ 8Þð2n þ 10Þ: Hence we obtain
Pð3Þ ¼ ðn þ 3ÞI ð2n þ 12ÞE;
Rð3Þ ¼ ðn þ 1Þðn þ 2ÞI ð2n þ 4Þð2n þ 8ÞE þ ð2n þ 8Þð2n þ 10ÞE2 :
ð2:21Þ
2.2. Algorithm Now, we present an efficient algorithm of performing the Q transformation, determined by the operators PðmÞ , hence by the annihilators LðmÞ obtained by the method of Section 2.1. Let us denote the numerators and denominators of the right side ðmÞ :¼ LðmÞ ðsn Þ and DðmÞ :¼ LðmÞ ð1Þ (notice the dependence on n, since LðmÞ annihilates rn ), respectively. By (2.9), of (2.3) by N ðmÞ n n we have the following recurrence relations:
NðmÞ ¼ PðmÞ ðNðm1Þ Þ; n n
DðmÞ ¼ PðmÞ ðDðm1Þ Þ: n n
Thus, the algorithm of computing QðmÞ does not need the operators LðmÞ explicitly; only PðmÞ are needed. We assume the inin tialization Pð1Þ :¼ Lð1Þ . P1 ðmÞ Algorithm 1. Consider the series ; m 2 N, are such that j¼0 aj with partial sums sn . If the operators P ðmÞ ðmÞ ðmÞ ðmÞ ðm1Þ ð1Þ :¼ P P P annihilate r n , then array Qn can be computed in the following recurrent way: L
Nð0Þ n ¼ sn ;
Dð0Þ n ¼ 1;
N ðmÞ ¼ PðmÞ ðNðm1Þ Þ; n n QðmÞ :¼ n
NðmÞ n DðmÞ n
DðmÞ ¼ PðmÞ ðDnðm1Þ Þ; m P 1; n
:
In the following example, we are able to obtain the operators PðmÞ explicitly. Example 2.2. Consider the series (2.18) once more. We prove that continuing the construction described in Example 2.1 gives the operators PðmÞ of the form
PðmÞ :¼ ðn þ mÞI ð2n þ 4mÞE;
m 2 N;
ð2:22Þ
(cf. (2.19)–(2.21)). Letting LðmÞ :¼ PðmÞ Pðm1Þ Pð1Þ , one can check that
LðmÞ ðanþk Þ ¼
m! ðkÞm 2
nþk
ðn þ k þ 1Þmþ1
;
k 2 N0 ; m 2 N;
ð2:23Þ
where ðzÞk is the Pochhammer symbol defined for any complex number z by
ðzÞ0 :¼ 1;
ðzÞk :¼
k1 Y
ðz þ jÞ; k 2 N:
ð2:24Þ
j¼0
The proof of (2.23) is by induction on m 2 N. For any k 2 N0 we have
Lð1Þ ðanþk Þ ¼ ððn þ 1ÞI ð2n þ 4ÞEÞðanþk Þ ¼ ðn þ 1Þanþk ð2n þ 4Þanþkþ1 ¼
ðn þ 1Þðn þ k þ 2Þ ðn þ 2Þðn þ k þ 1Þ 2nþk ðn þ k þ 1Þ2
;
and (2.23) holds for m ¼ 1. Let us suppose that it holds for some natural number m 2 N. Then for any k 2 N0 we obtain
Lðmþ1Þ ðanþk Þ ¼ Pðmþ1Þ LðmÞ ðanþk Þ ¼ Pðmþ1Þ ¼
!
m! ðkÞm 2nþk ðn þ k þ 1Þmþ1
¼
m! ðkÞm ðn þ m þ 1Þ 2nþk ðn þ k þ 1Þmþ1
m! ðkÞm ðn þ 2m þ 2Þ 2nþk ðn þ k þ 2Þmþ1
m! ðkÞm
2
nþk
½ðn þ m þ 1Þðn þ k þ m þ 2Þ ðn þ 2m þ 2Þðn þ k þ 1ÞÞ : ðn þ k þ 1Þmþ2 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} ðmþ1ÞðmkÞ
This completes the proof. Since ðkÞm ¼ 0 for k < m, one can get that LðmÞ ðanþk Þ ¼ 0 for all k ¼ 0; 1; . . . ; m 1, and consequently that (2.2) holds. ðmÞ Hence, computing of Qn can be executed by Algorithm 1, i.e., in the following recurrent way:
P. Woz´ny, R. Nowak / Applied Mathematics and Computation 215 (2009) 1622–1645 n1 X
Nð0Þ n :¼ sn ¼
aj ;
1627
Dð0Þ n :¼ 1;
j¼0 ðm1Þ
:¼ PðmÞ ðNðm1Þ Þ ¼ ðn þ mÞNnðm1Þ ð2n þ 4mÞNnþ1 ; NðmÞ n n
m P 1;
ðm1Þ 4mÞDnþ1 ;
m P 1;
DðmÞ n
ðmÞ
:¼ P
:¼ QðmÞ n
ðDðm1Þ Þ n
NðmÞ n DðmÞ n
¼ ðn þ
mÞDðm1Þ n
ð2n þ
:
Given the partial sums s1 ; s2 ; . . . ; s5 , one can obtain the following array of quantities QðmÞ (underlined digits are exact): n ð0Þ
ð1Þ
ð2Þ
ð3Þ
ð0Þ
ð1Þ
ð2Þ
ð3Þ
ð0Þ
ð1Þ
ð2Þ
ð0Þ
ð1Þ
ð4Þ
Q1 ¼ 1:000 Q1 ¼ 1:37500 Q1 ¼ 1:38596491 Q1 ¼ 1:38628472 Q1 ¼ 1:38629408 Q2 ¼ 1:250 Q2 ¼ 1:38333 Q2 ¼ 1:38621795 Q2 ¼ 1:38629227 Q3 ¼ 1:333 Q3 ¼ 1:38542 Q3 ¼ 1:38627451 Q4 ¼ 1:365 Q4 ¼ 1:38601 ð0Þ Q5
¼ 1:377 ð4Þ
Let us observe that Q1 is the following combination of the first 5 partial sums: ð4Þ
Q1 ¼
5 40 280 2240 672 s1 s2 þ s3 s4 þ s5 : 501 167 167 501 167 ð0Þ
Thus our method gives five exact digits more than partial sum s5 does (cf. Q5 ). Using Wynn’s e-algorithm (see, e.g., [3, Section2.3]) and Aitken’s iterated D2 process (see, e.g., [1, (5.1–15)] or [5, pp. 149–152]) one obtains ð1Þ eð1Þ A2 ¼ 1:38611111: 4 ¼ 1:38596491;
2.3. Relation to the Wynn’s epsilon algorithm Wynn’s e-algorithm (see, e.g., [3, Section2.3] or [1, Section4]) is defined by
eðnÞ eðnÞ 1 ¼ 0; 0 ¼ sn ; ðnÞ ðnþ1Þ 1 ekþ1 ¼ ek1 þ ½ekðnþ1Þ eðnÞ n; k 2 N0 k ; and solves the system of equations
sr ¼ ^s þ
m 1 X
cj Dsrþj ;
r ¼ n; n þ 1; . . . ; n þ m;
ð2:25Þ
j¼0 ðnÞ
being a model of Shanks transformation (see [14], or [13, p. 78]), with unknown ^s and coefficients cj , in the sense that ^s ¼ e2m . Consider the series (2.1) and assume that Lð1Þ is of the form
Lð1Þ ¼ k0 ðnÞI þ k1 ðnÞE: ð2Þ
ð3Þ
ð2:26Þ ðmÞ
ðmÞ
Let us denote by L ; L ; . . . ; L the operators obtained by the method of Section 2.1. By (2.17), we get ord L ¼ m, and thus one can solve the system of Eqs. (2.25) by the following trick. Applying the operator LðmÞ to both sides of the first equation of (2.25), we obtain ðmÞ
L
ðsn Þ ¼ L
ðmÞ
ð^sÞ þ LðmÞ
m1 X
! cj Dsnþj :
j¼0
One can prove that the last term in the right-hand side vanishes as a consequence of (2.2). Thus, we get
^s ¼
LðmÞ ðsn Þ : LðmÞ ð1Þ
Thus we obtain the following relation of the proposed method to Wynn’s epsilon algorithm: ðnÞ
QðmÞ ¼ e2m : n
ð2:27Þ
Example 2.3. Let us consider the series (2.18). It is quite fast convergent, since:
accðs5 Þ ¼ 2:2;
accðs15 Þ ¼ 5:6;
accðs25 Þ ¼ 8:8;
accðs35 Þ ¼ 11:9; ð1Þ
accðs45 Þ ¼ 15:1: ð1Þ
However, let us use only the partial sums s1 ; s2 ; . . . ; s15 . Observe that L ¼ P is, in view of (2.22), of the form (2.26). The ð1Þ ð1Þ ð1Þ first row of the array obtained by e-algorithm contains the quantities e0 ; e1 ; . . . ; e14 . Algorithm 1 yields the quantities
P. Woz´ny, R. Nowak / Applied Mathematics and Computation 215 (2009) 1622–1645
1628
Table 1 Values of accðQðmÞ n Þ. Results obtained by e-algorithm (above the stepping line) are marked in red.
ð0Þ
ð1Þ
ð14Þ
ð14Þ
ð1Þ
Q1 ; Q1 ; . . . ; Q1 (cf. Table 1), which is more than e-algorithm gives, since Q1 ¼ e28 . The reason is that we have the formulae (2.22), which are available thanks to the knowledge of analytical form of the terms an , as described in Section 2.1. Wynn’s method does not need such information.
3. Generalized hypergeometric series Now, we present the application of the proposed method to the case of generalized hypergeometric series
!
pþ1 Fp
a0 ; a1 ; . . . ; ap
x b1 ; b2 ; . . . ; bp
:¼
1 X ða0 Þn ða1 Þn ða2 Þn ðap Þn n x ðb1 Þn ðb2 Þn ðbp Þn n! n¼0
(see, e.g., [15, Section2.1]), where ðzÞn is the Pochhammer symbol (see (2.24)). 3.1. Main results Without loss of generality we can assume that a0 :¼ 1, since one can write
!
pþ1 Fp
a0 ; a1 ; . . . ; ap
x b1 ; b2 ; . . . ; bp
¼ pþ2 Fpþ1
! 1; a0 ; a1 ; . . . ; ap x : 1; b1 ; b2 ; . . . ; bp
In Theorem 2, we present the explicit form of operators PðmÞ which are the most important input in Algorithm 1 of computing the proposed transformation in the case of considered series. We will show that our method of summation can be ˇ ízˇek et al. , recently introduced in [2], and analysed by Weniger identified as some special variant of the transformation of C in [9]. Numerical examples show that if we choose the second or higher order operator Pð1Þ , then our method is not equivalent to e-algorithm. Given a vector c ¼ ðc1 ; c2 ; . . . ; cp Þ, we will use the notation
Ac þ B :¼ ðAc1 þ B; Ac2 þ B; . . . ; Acp þ BÞ;
A; B 2 C;
and
hcin :¼
p Y
ðcj Þn ;
hci :¼ hci1 ;
n 2 N0 :
ð3:1Þ
j¼1
Hence, we write the considered series as
pþ1 Fp
! 1 X 1; a1 ; a2 ; . . . ; ap an x ¼ b1 ; b2 ; . . . ; bp n¼0
with an :¼
hain n ða1 Þn ða2 Þn ðap Þn n x ¼ x : hbin ðb1 Þn ðb2 Þn ðbp Þn
ð3:2Þ
P. Woz´ny, R. Nowak / Applied Mathematics and Computation 215 (2009) 1622–1645
1629
Theorem 2. Let us define
hbin n x I ; hain p X mp :¼ ðDj hb þ n þ mðp þ 1Þ j 2iÞDpj ; j j¼0
Pð1Þ :¼ Dp
PðmÞ
ð3:3aÞ m ¼ 2; 3 . . . :
ð3:3bÞ
Then the operators LðmÞ :¼ PðmÞ Pðm1Þ Pð1Þ can be written in the form
LðmÞ ¼ Dmp
hbinþm1 I ; hain xn
m 2 N;
ð3:4Þ
LðmÞ ðan þ anþ1 þ þ anþm1 Þ ¼ 0:
ð3:5Þ
and satisfy
Proof. By the following well-known property of operator D:
Dr ðfn IÞ ¼
r X r ðDj fnþrj ÞDrj ; j j¼0
r 2 N0 ;
ð3:6Þ
(cf. [16, Theorem 1.3.2]), we obtain
Dmp ðhb þ n þ m 2iIÞ ¼ PðmÞ Dðm1Þp ;
m ¼ 2; 3; . . . ;
PðmÞ being given by (3.3b). Hence, (3.4) follows by induction on m. To prove (3.5), notice that LðmÞ ðanþk Þ ¼ 0 for k ¼ 0; 1; . . . ; m 1, since
hbinþm1 anþk ¼ ha þ nik hb þ n þ kim1k xk hain xn is a polynomial in n of degree ðm 1Þp. h The following two corollaries give the form of operators PðmÞ with coefficients being polynomials in n. Corollary 3 (Gauss series). Given the series 2 F1
X 1 1; a x ¼ an b
with an :¼
n¼0
ðaÞn n x ; ðbÞn
let us define
b ðmÞ :¼ xðn þ m þ a 1ÞI ðn þ 2m þ b 2ÞE; P b ðmÞ
Then the operators L
b ðmÞ L
b ðmÞ
:¼ P
b ðm1Þ
P
b ð1Þ
P
m 2 N:
satisfy
ðan þ anþ1 þ þ anþm1 Þ ¼ 0:
ð3:7Þ
b ðmÞ in the sense that Proof. One can check that the operators PðmÞ , given in Theorem 2 (in the case of p ¼ 1), are related to P
b ðkÞ ¼ PðkÞ Qðk1Þ ; QðkÞ P
k 2 N;
where
Qð0Þ :¼ I;
QðkÞ :¼
ðbÞnþk1 I; k 2 N: ðaÞnþk xnþk
Thus LðmÞ ¼ QðmÞ b L ðmÞ , m 2 N, and consequently (3.7) holds. h Corollary 4 (Series 3 F 2 ). Given the series 3 F2
X 1 1; a1 ; a2 x ¼ an b ;b 1
2
n¼0
with an :¼
ða1 Þn ða2 Þn n x ; ðb1 Þn ðb2 Þn
let us define
b ðmÞ :¼ x2 hn þ 2m þ a 2i I 2xhn þ 2m þ a 1i½hn þ 2m þ b 2i P 2 mðm 1Þ E þ hn þ m þ b 1ihn þ 3m þ b 2i E2
ð3:8Þ
P. Woz´ny, R. Nowak / Applied Mathematics and Computation 215 (2009) 1622–1645
1630
b ðmÞ P b ð1Þ satisfy b ðm1Þ P with a ¼ ða1 ; a2 Þ; b ¼ ðb1 ; b2 Þ. Then the operators b L ðmÞ :¼ P
b L ðmÞ ðan þ anþ1 þ þ anþm1 Þ ¼ 0:
ð3:9Þ
Proof. The proof follows very closely the previous proof, since the operators PðmÞ , given in Theorem 2 (in the case of p ¼ 2), b ðmÞ in the sense that are related to P
b ðkÞ ¼ PðkÞ Qðk1Þ ; QðkÞ P
k 2 N;
where
Qð0Þ :¼ I;
QðkÞ :¼
ðb1 Þnþk1 ðb2 Þnþk1 I; ða1 Þnþ2k ða2 Þnþ2k xnþ2k
k 2 N:
Thus LðmÞ ¼ QðmÞ b L ðmÞ ; m 2 N, and hence (3.9) holds. h Example 3.1. Consider the series 3 F2
! 1 X 1; 12 ; 12 ð1Þn 1 ¼ 2 3 3 ; 2 2 n¼0 ð2n þ 1Þ
ð3:10Þ
summing up to the Catalan’s constant G 0:9159655941772190. We compared the results obtained by:
eðnÞ 2k – Wynn’s e-algorithm (see [1, (4.2-2)]), ðnÞ
Ak – Aitken’s iterated D2 process (see [1, (5.1-15)]), ðkÞ sn – The method proposed by Lewanowicz and Paszkowski in [17, Theorem 3.2], ðkÞ t u LðkÞ n ð1; fsn gÞ and Ln ð1; fsn gÞ – t and u variants of Levin transformation (see [1, (7.1-7), (7.3-5), (7.3-7)], Q – transformation (2.3) determined by operators (3.8). The input for all algorithms was the set of the partial sums s1 ; s2 ; . . . ; s15 of series (3.10). In Table 2, we present the values: ðkÞ
ð1Þ
ð1Þ
ð2kÞ
ð2kÞ
accðQ1 Þ; accðe2k Þ; accðAk Þ; accðt L1 ð1; fsn gÞÞ; accðu L1 ð1; fsn gÞÞ ðkÞ
for k ¼ 0; 1; . . . ; 7, respectively. In the case of method of Lewanowicz and Paszkowski, we present the antidiagonal skþ1 . Since the series (3.10) is alternating, e-algorithm and Aitken’s D2 process are quite efficient. Example 3.2. Consider the following positive series:
!
A ¼ 2 F1
1 1 ; 4 2 1 5 4
¼ 3 F2
! 1 X 1; 14 ; 12 an ; 1 ¼ 5 1; 4 n¼0
1 1 an :¼ 4 n 25 n ; ð1Þn 4 n
summing up to lemniscate constant A 1:31102877714605990523. In this case, e-algorithm is ineffective; we have, for inð1Þ stance, accðe14 Þ ¼ 1:61. The convergence of the series can be successfully accelerated by the u variant of Levin transformað14Þ tion, or the method of Lewanowicz and Paszkowski (see [17, Theorem 2.2]). Namely, we have accðu L1 ð1; fsn gÞÞ ¼ 11:50 ð7Þ and accðs8 Þ ¼ 13:88, respectively. The Q transformation, determined by operators (3.8), yields similar results; we have ð7Þ accðQ1 Þ ¼ 13:03. Notice that ord Lð1Þ ¼ 2, and much better results than e-algorithm (corresponding to the choice of first-order operator Lð1Þ ; see Section 2.3) are obtained. Thus the order of operator Lð1Þ , being of choice in the Section 2.1, has a great importance indeed.
Table 2 Accuracy of approximations to the sum of the series (3.10). Comparison of the proposed summation method, e-algorithm, Aitken’s iterated D2 process, method of Lewanowicz and Paszkowski [17], t and u variants of Levin transformation. k
0
1
2
ð1Þ accð 2k Þ ð1Þ accðAk Þ ðkÞ accðskþ1 Þ ð2kÞ accðt L1 ð1; fsn gÞÞ ð2kÞ accðu L1 ð1; fsn gÞÞ ðkÞ accðQ1 Þ
1.04
2.59
4.14
e
3
4 5.67
5
7.21
8.74
6
7
10.27
11.81
1.04
2.59
4.40
6.40
8.60
13.06
11.81
14.46
1.04
3.00
4.89
6.76
8.62
10.46
12.30
14.14
1.04
4.91
5.76
11.24
12.11
14.23
16.91
20.15
1.04
2.70
5.10
8.78
10.95
13.95
16.61
18.78
1.04
3.44
6.20
9.55
12.27
14.95
18.15
21.19
P. Woz´ny, R. Nowak / Applied Mathematics and Computation 215 (2009) 1622–1645
1631
3.2. Some theoretical results Theorem 5. The Q transformation, applied to the generalized hypergeometric series (3.2) with x–1, is regular, i.e.,
! 1; a1 ; a2 ; . . . ; ap for all m 2 N: x b1 ; b2 ; . . . ; bp
lim QðmÞ n ðsn Þ ¼ pþ1 Fp
n!1
Proof. Let LðmÞ ; m 2 N, denote the operators (3.4). We have
QðmÞ n
LðmÞ ðsn Þ ¼ ¼ ðmÞ L ð1Þ
hainþmp xnþmp hbinþm1
LðmÞ ðsn Þ
hainþmp xnþmp hbinþm1
ðmÞ
L
ð1Þ
Pmp
ðmpÞ
j¼0 kj Pmp
¼
ðn; xÞ snþj
ðmpÞ
j¼0 kj
ðn; xÞ
ð3:11Þ
;
where
hainþmp xnþmp ðmpÞ kj ðn; xÞ :¼ hbinþm1
" mpj
ð1Þ
# mp mp hbinþjþm1 ¼ ðxÞmpj ha þ n þ jimpj hb þ n þ m 1ij ; j ¼ 0;1; .. .; mp: nþj hainþj x j j ðmpÞ
Denote the denominator in the last leg of (3.11) by Mðn; xÞ. Since kj ðn; xÞ are polynomials in n of degree mp2 with leading mp coefficient ðxÞmpj , under assumption x–1, the function Mðn; xÞ is also a polynomial in n of degree mp2 ; more j precisely 2
Mðn; xÞ ¼ ð1 xÞmp nmp þ Oðnmp
2 1
n ! 1:
Þ;
Let sn denote the partial sums of the series (3.2) and let limn!1 sn ¼ s. Then
s
QðmÞ n
LðmÞ ðsÞ LðmÞ ðsn Þ ¼ ¼ ðmÞ L ð1Þ LðmÞ ð1Þ
hainþmp xnþmp hbinþm1
LðmÞ ðs sn Þ
xnþmp
hainþmp hbinþm1
LðmÞ ð1Þ
Pmp ¼
ðmpÞ ðn; xÞ ðs snþj Þ j¼0 kj Pmp ðmpÞ ðn; xÞ j¼0 kj
¼
ðmpÞ mp X kj ðn; xÞ
Mðn; xÞ
j¼0
ðs snþj Þ ! 0; n!1
because ðmpÞ
kj
ðn; xÞ
Mðn; xÞ
¼ Oð1Þ;
j ¼ 0; 1; . . . ; mp; n ! 1:
In the following theorem, we present an efficiency of our transformation in the case of generalized hypergeometric series. Theorem 6. pþ1 Fp
! 1; a1 ; a2 ; . . . ; ap nþmðpþ1Þ Þ; x QðmÞ n ðsn Þ ¼ Oðx b1 ; b2 ; . . . ; bp ðmpÞ
Proof. Let us use the functions kj satisfying
Mðn; xÞ ¼ Oð1Þ;
x ! 0:
ðn; xÞ and Mðn; xÞ introduced in the proof of Theorem 5. Now, Mðn; xÞ is a function of x
x ! 0:
Similarly, we obtain ðmÞ
s
QðmÞ n
LðmÞ ðs sn Þ L ¼ ¼ LðmÞ ð1Þ ¼
ðmpÞ mp X kj ðn; xÞ j¼0
Mðn; xÞ
P
m1 k¼0 anþk ðmÞ
L
þ rnþm
ð1Þ
ðmpÞ
ðmpÞ
ðn; xÞ ¼ Oðxmpj Þ;
LðmÞ ðr nþm Þ ¼ ðmÞ ¼ L ð1Þ
hainþmp xnþmp hbinþm1
LðmÞ ðrnþm Þ
hainþmp xnþmp hbinþm1
LðmÞ ð1Þ
r nþjþm :
This completes the proof, since functions kj
kj
r nþmþj ¼
1 X i¼0
ðn; xÞ and r nþmþj satisfy
hainþmþiþj nþmþiþj x ¼ Oðxnþmþj Þ; x ! 0: hbinþmþiþj
Pmp ¼
ðmpÞ ðn; xÞ r nþjþm j¼0 kj
Mðn; xÞ
P. Woz´ny, R. Nowak / Applied Mathematics and Computation 215 (2009) 1622–1645
1632
Remark 7. The authors were unable to show that Q transformation accelerates the convergence in the sense that
lim
n!1
Qðmþ1Þ ðsn Þ s n ¼ 0: ðmÞ Qn ðsn Þ s
This seems to be a very difficult problem, in general. 3.3. Relation to the Cízˇek et al. transformation By (2.3) and (3.4), the Q transformation in the case of generalized hypergeometric series can be written in the following form:
QðmÞ n
i h sn Dmp ðn þ b1 Þm1 ðn þ b2 Þm1 ðn þ bp Þm1 h i ¼ h ¼ hbi Dmp hainþm1 Dmp ðn þ b1 Þm1 ðn þ b2 Þm1 ðn þ bp Þm1 xn
Dmp
h
hbinþm1 hain xn n
sn an 1 an
i i:
ð3:12Þ
Remind that the G transformation of Cízˇek et al. is defined by m1 GðmÞ n ðfqk gk¼1 ; fsn g; f
xn gÞ :¼
Dm Dm
hQ
m1 k¼1 ðn
þ qk Þ xsnn
m1 k¼1 ðn
þ qk Þ x1n
hQ
i i;
ð3:13Þ
where fqk gm1 k¼1 is a set of parameters, and xn is a sequence of remainder estimates (see [2] or [1, (2.13)]). In special cases qk 1, or qk ¼ k (for k ¼ 1; 2; . . . ; m 1) the Levin and Weniger [1, Section8.2] transformations are obtained, respectively. Now, it is easy to observe that the transformation (3.12) is equivalent to the transformation (3.13) with some particular choice of parameters and remainder estimates. Namely, we have the following. Proposition 8. Equation
m 1 QðmÞ ¼ Gnðm Þ ðfqk gk¼1 ; fsn g; fxn gÞ n
holds, where
m :¼ mp; n
x :¼ ðn þ 1Þ
m 1 fqk gk¼1
ð3:14aÞ p1
an ; ð3:14bÞ 9 8 > > = < :¼ 1; 1; . . . ; 1; b1 ; b2 ; . . . ; bp ; b1 þ 1; b2 þ 1; . . . ; bp þ 1; . . . ; b1 þ m 2; b2 þ m 2; . . . ; bp þ m 2 : ð3:14cÞ |fflfflfflfflfflffl ffl {zfflfflfflfflfflffl ffl } > > ; : p1 times
Example 3.3. For the series
s ¼ 5 F4
1; 1; 2; 2; 7 99 1:0376328566238592296948; 3; 5; 6; 9 100
we obtain ð5Þ
Q1 1:03763285662385922975208; ð5Þ
ð5Þ
ð20Þ
ð20Þ
which yields accðQ1 Þ ¼ 19:26. Notice that Q1 ¼ G1 ðfqk g; fsn g; fxn gÞ; see Proposition 8 ðp ¼ 4Þ. We computed G1 and compared the choices of parameters fqk g and remainder estimates fxn g given in (3.14) with the classic ones. As shown in Table 3, the best results are obtained with the parameters chosen as in (3.14), i.e., by using Q transformation (cf. (3.12)).
Table 3 ð20Þ Values of accðG1 ðfqk g19 k¼1 ; fsn g; fxn gÞÞ. qk
(3.14c) 1 (Levin L [1, (7.1-7), b ¼ 1]) k (Weniger S [1, (8.2-7), b ¼ 1]) 2
k (cf. [2, Tables 1 and 2])
xn (3.14b)
an1 (t-variant)
ðn þ 1Þan1 (u-variant)
an aann1 an1 (v-variant)
19.26 18.44 14.20 1.51
15.46 14.46 16.04 10.60
16.78 15.83 17.33 8.39
17.67 16.69 18.24 9.73
P. Woz´ny, R. Nowak / Applied Mathematics and Computation 215 (2009) 1622–1645
1633
4. Basic hypergeometric series The Q transformation may also be successfully applied to the case of the basic hypergeometric series, or q-hypergeometric series,
pþ1 /p
!
a0 ; a1 ; . . . ; ap
q; x b1 ; b2 ; . . . ; bp
1 X ða0 ; qÞn ða1 ; qÞn ða2 ; qÞn ðap ; qÞn n x ðb1 ; qÞn ðb2 ; qÞn ðbp ; qÞn ðq; qÞn n¼0
:¼
(see, e.g., [18, (1.2.22)]), where ðz; qÞk is the q-Pochhammer symbol, defined for any complex number z by
8 > <
ðz; qÞk :¼
1; k1 Q
> :
k ¼ 0;
ð1 zqj Þ; k > 0:
j¼0
Without loss of generality we may assume that a0 :¼ q, since one can write
pþ1 /p
a0 ; a1 ; . . . ; ap
!
q; x b1 ; b2 ; . . . ; bp
¼ pþ2 /pþ1
! q; a0 ; a1 ; . . . ; ap q; x : q; b1 ; b2 ; . . . ; bp
We append the notation (3.1) by
hc; qin :¼
p Y ðcj ; qÞn ; j¼1
and write the considered series shortly as pþ1 /p
! 1 X q; a1 ; . . . ; ap an q; x ¼ b1 ; b2 ; . . . ; bp n¼0
with an :¼
ha; qin n x : hb; qin
ð4:1Þ
Before we propose a q-analogue of Theorem 2, let us introduce the operator Dq and its powers by
D0q :¼ I; Drq :¼ ðE qr1 IÞDr1 q ;
r 2 N:
One can check, that q-analogue of (3.6) is
Drq ðfn IÞ ¼
r X r j¼0
where
j
ðDjq fnþrj ÞDrj q ;
r 2 N0 ;
ð4:2Þ
q
r is the Gauss binomial coefficient, defined by j q
r j
:¼
q
ðq; qÞr ; ðq; qÞj ðq; qÞrj
r; j 2 N0 ;
j 6 r:
Theorem 9. Let us define
hb; qin n x I ; ha; qin p
X mp ~ ðmpp;pjÞ ; :¼ ðDjq h1 b qnþmðpþ1Þj2 iÞD q j q j¼0
Pð1Þ :¼ Dpq PðmÞ
m ¼ 2; 3 . . . ;
where
~ ði;0Þ :¼ I; D q
~ ði;rÞ :¼ D ~ ðiþ1;r1Þ ðE qi IÞ; D q q
Then the operators L
LðmÞ
ðmÞ
ðmÞ
:¼ P
P
ðm1Þ
hb; qinþm1 ¼ Dmp I ; q ha; qin xn
ð1Þ
P
m 2 N;
i 2 Z; r 2 N: are of the form
ð4:3Þ
and satisfy
LðmÞ ðan þ anþ1 þ þ anþm1 Þ ¼ 0:
ð4:4Þ
P. Woz´ny, R. Nowak / Applied Mathematics and Computation 215 (2009) 1622–1645
1634
Table 4 Values of accðQðmÞ n Þ for the approximations of the sum of the series (4.5), with a ¼ 9=10; b ¼ 5=8; c ¼ 1=2 and q ¼ 1=5. nnm
0
1
2
3
4
5
6
1 2 3 4 5 6 7 8 9 10 11 12 13
0.42 0.48 0.53 0.59 0.64 0.69 0.74 0.79 0.84 0.89 0.94 1.00 1.05
3.18 4.67 6.13 7.58 9.03 10.48 11.93 13.38 14.83 16.28 17.73
8.75 11.64 14.50 17.35 20.19 23.04 25.89 28.73 31.58
14.77 18.36 21.92 25.47 29.01 32.56 36.11
22.17 26.46 30.72 34.97 39.21
30.97 35.96 40.91
41.15
Proof. The proof follows very closely the proof of Theorem 2. Namely, using the the property (4.2), we obtain nþk2 Dkp iIÞ ¼ PðkÞ Dqðk1Þp ; q ðh1 b q
k ¼ 2; 3; . . . ; m 1;
and thus the proof of (4.3) follows by induction on m. Next, observe that
hb; qinþm1 anþk ¼ haqn ; qik hbqnþk ; qim1k xk ; ha; qin xn
k ¼ 0; 1; . . . ; m 1;
is a polynomial in qn of degree ðm 1Þp. Since operator Drq annihilates all the polynomials in qn of degree at most r 1, we have that LðmÞ ðanþk Þ ¼ 0 for all k ¼ 0; 1; . . . ; m 1. This completes the proof of (4.4). h Summarizing, in the case of q-hypergeometric series (4.1), we can compute the Q transformation by Algorithm 1 with operators PðmÞ given in Theorem 9. Such a transformation does not seem to be related to any classic method of convergence acceleration. Example 4.1. Let us consider the following q-hypergeometric series:
2 /1
a; b q; a; b q; c=ab ¼ 3 /2 q; c=ab : c q; c
ð4:5Þ
The Heine’s q-analogue of Gauss summation formula is given by
2 /1
a; b ðc; =a; qÞ1 ðc=b; qÞ1 q; c=ab ¼ ðc; qÞ ðc=ab; qÞ c 1
1
(see, e.g., [18, (1.5.1)]), where ðz; qÞ1 is the generalized q-Pochhammer symbol, defined by
ðz; qÞ1 ¼
1 Y ð1 zqj Þ: j¼0
Let us put a ¼ 9=10; b ¼ 5=8; c ¼ 1=2 and q ¼ 1=5. Then the sum s of the series equals approximately 1.6172358319254889075. Let us notice the very slow convergence of the considered series (cf. the second column of Table 4. We used only the partial sums s1 ; s2 ; . . . ; s13 ; and the operators PðmÞ ; m ¼ 1; 2; . . . ; 6; given in Theorem 9. The results obtained are presented in Table 4.
5. Orthogonal polynomial series Let us consider the series 1 X
an ;
an :¼ bn f n ;
ð5:1Þ
n¼0
where fn satisfies a linear difference equation of the second order
jn f n þ kn f nþ1 þ ln f nþ2 ¼ 0; n 2 N0 :
ð5:2Þ
It is known that Eq. (5.2) is satisfied, if fn ¼ W n ðxÞ, where W n ðxÞ are orthogonal polynomials in an interval I with respect to a weight function xðxÞ. More specifically, if W n is
P. Woz´ny, R. Nowak / Applied Mathematics and Computation 215 (2009) 1622–1645
1635
(I) the nth Chebyshev polynomial of the first kind T n ðI ¼ ½1; 1; xðxÞ ¼ ð1 x2 Þ1=2 Þ or of the second kind U n ðI ¼ ½1; 1; xðxÞ ¼ ð1 x2 Þ1=2 Þ, then
jn ¼ 1; kn ¼ 2x; ln ¼ 1;
ð5:3Þ
(II) the nth Legendre polynomial P n ðI ¼ ½1; 1; xðxÞ 1Þ, then
jn ¼ n þ 1; kn ¼ ð2n þ 3Þx; ln ¼ n þ 2; (III) the nth Laguerre polynomial
LnðaÞ ðI
ð5:4Þ
¼ ½0; 1Þ; xðxÞ ¼
xa ex ; a
> 1), then
jn ¼ n þ 1 þ a; kn ¼ ð2n þ 3 þ a xÞ; ln ¼ n þ 2; x2
(IV) the nth Hermite polynomial Hn ðI ¼ ð1; 1Þ; xðxÞ ¼ e
ð5:5Þ Þ, then.
jn ¼ n þ 1; kn ¼ x; ln ¼ 1=2;
ð5:6Þ
see, e.g., [8, p. 36], or [19, Sections1.13, 1.11, 1.8.3, 1.8.2]. The case of Fourier series
s ¼ A0 =2 þ
1 X ðAn cos na þ Bn sin naÞ
ð5:7Þ
n¼1
corresponds to the case I, since the functions fn ¼ sin na; cos na satisfy the Eq. (5.2) with
jn ¼ 1; kn ¼ 2 cos a; ln ¼ 1:
ð5:8Þ
The summation of such series was studied, among others, by Homeier [10]. He proposed the H transformation, defined as follows: 1 Nð0Þ sn =xn ; n ¼ ðn þ bÞ
1 Dð0Þ n ¼ ðn þ bÞ =xn ; ðk1Þ
ðk1Þ
ðk1Þ NðkÞ 2 cos a ðn þ k þ bÞNnþ1 þ ðn þ 2k þ bÞN nþ2 ; n ¼ ðn þ bÞN n
DðkÞ n
¼ ðn þ
bÞDðk1Þ n
2 cos a ðn þ k þ
ðk1Þ bÞDnþ1
þ ðn þ 2k þ
ðk1Þ bÞDnþ2 ;
ð5:9Þ
ðkÞ ðkÞ HðkÞ n ða; b; fsn g; fxn gÞ ¼ N n =Dn ;
see, e.g., [4, (155)]. The transformation of the Fourier series (5.7) to accelerate their convergence, especially near the singularities, was proposed by Oleksy [20]. The hierarchically consistent iteration of Hð1Þ transformation was studied by Homeier [11] in derivation of his I transformation. More generally, in the case of orthogonal series (5.1), hierarchically consistent iteration of the transformation
s0n :¼
jn xsnn þ kn xsnþ1 þ ln xsnþ2 nþ1 nþ2 1 1 jn x1n þ kn xnþ1 þ ln xnþ2
ð5:10Þ
yields K transformation, defined as follows:
Nð0Þ Dð0Þ n ¼ sn =xn ; n ¼ 1=xn ; 1 ðk1Þ ðk1Þ ðk1Þ þ knþk1 Nnþ1 þ lnþk1 Nnþ2 Þ; N ðkÞ n ¼ ðkÞ ðjnþk1 N n dn 1 ðk1Þ ðk1Þ ðk1Þ þ knþk1 Dnþ1 þ lnþk1 Dnþ2 Þ; DðkÞ n ¼ ðkÞ ðjnþk1 Dn dn
ð5:11Þ
ðkÞ ðkÞ ðkÞ KðkÞ n ðfsn g; fxn g; fdn g; fjn ; kn ; ln gÞ ¼ N n =Dn ;
see, e.g., [12], or [4, (174)]. The transformation (5.10) is exact for the sequences of the form
sn ¼ s þ xn ðc Pn þ d Q n Þ; where c and d are arbitrary constants, while P n and Q n are two linearly independent solution of the Eq. (5.2). For instance, it is exact in the case
sn þ bn fn ¼ s: Let us notice that QðmÞ transformation is exact for any sequence sn satisfying
sn þ rðmÞ ¼ sn þ bn fn þ bnþ1 fnþ1 þ þ bnþm1 fnþm1 ¼ s: n
ð5:12Þ
Since Homeier’s transformations are Levin-type sequence transformations, they can be performed by the Algorithm 1 with some special operators PðmÞ . Notice that such operators do not correspond to the model (5.12), since usually the ðmÞ operators LðmÞ :¼ PðmÞ Pðm1Þ Pð1Þ do not annihilate rn ; cf. (2.2). Indeed, in the case of orthogonal series (5.1), we usually ðmÞ ðmÞ are not able to obtain explicit forms of operators P , and so the operators LðmÞ ; m 2 N, annihilating r n . In the sequel, we
P. Woz´ny, R. Nowak / Applied Mathematics and Computation 215 (2009) 1622–1645
1636
consider some special cases of the series (5.1), where we found operators PðmÞ explicitly for m 2 N. Some numerical computations, comparing Q transformation with Homeier’s methods and the methods recently introduced by Paszkowski [8], are given. In the following theorem, we present an explicit form of operators PðmÞ in the case of some Hermite series. Theorem 10. Let z 2 C, and
an :¼
1 Hn ðxÞ; zn n!
n 2 N0 :
Let us define the operators
PðmÞ :¼ I zxE þ z2 ðn þ 3m 1Þ=2E2 ;
m 2 N:
Then the operators LðmÞ :¼ PðmÞ Pðm1Þ Pð1Þ satisfy
LðmÞ ðan þ anþ1 þ þ anþm1 Þ ¼ 0:
ð5:13Þ
Proof. We proceed by induction on m 2 N. One can check that
Lð1Þ ðan Þ ¼ Pð1Þ
1
zn n!
Hn
¼
1 ððn þ 1ÞHn ðxÞ xHnþ1 ðxÞ þ Hnþ2 ðxÞ=2Þ ¼ 0 zn ðn þ 1Þ!
(cf. (5.6)). Suppose that (5.13) holds for some m 2 N. Observe that it is enough to show that Lðmþ1Þ ðanþm Þ ¼ 0. Using the following property:
Pðkþ1Þ Ej ¼ Ej Pð1Þ þ 2ð3k jÞEjþ2 ;
j; k 2 N0 ;
we obtain
We can do the cancellations shown above, since Pð1Þ ðan Þ ¼ 0.
h
In Theorem 12, preceded by a lemma, we present the forms of operators PðmÞ ; m 2 N, in the case of the series (5.1) with some classes of sequences bn ; fn . Lemma 11. Let us define the function
rðnÞ as follows
ðA þ 1Þn rðnÞ :¼ 1; ðA þ 1 dÞn
ð5:14Þ
where A; d; 1 2 C. Further, in the case d ¼ 0 (resp. d–0), let us define the operators
b ðmÞ :¼ u ðn þ m þ AÞI þ w ðn þ 2m þ BÞE þ v ðn þ 3m þ CÞE2 ; P
ð5:15Þ
where u; w; v; B; C 2 C (resp. A ¼ B ¼ C) are chosen so that
b ð1Þ ðrðnÞ an Þ ¼ 0: P
ð5:16Þ
Then
b ðmÞ P b ð1Þ ½rðnÞðan þ anþ1 þ þ anþm1 Þ ¼ 0; b ð2Þ P P
m 2 N:
ð5:17Þ
b ðkÞ defined in (5.15) satisfy Proof. One can check that operators P
b ðkþ1Þ E2k P b ðkÞ Ek1 ¼ E1k P b ðkþ1Þ Ek P b ðkÞ ; P Let us define the operator
D :¼ ð1 þ dÞðuI þ wE þ vE2 Þ:
k 2 N:
ðaÞ
P. Woz´ny, R. Nowak / Applied Mathematics and Computation 215 (2009) 1622–1645
1637
Observe that
b ð1Þ ðrðnÞEÞ ¼ ðE P b ð1Þ DÞðrðnÞIÞ: P
ðbÞ
One can also check the following property:
b ðkþ1Þ D ¼ D P b ðkÞ ; P
k 2 N:
ðcÞ
The proof of (5.17) is by induction on m 2 N. By (5.16), it holds for m ¼ 1. Suppose that (5.17) holds for some m 2 N. Observing
b ð1Þ ½rðnÞðan þ anþ1 þ þ anþm Þ ¼ P b ðmþ1Þ P b ð2Þ P b ðmþ1Þ P b ðmÞ P b ðmÞ P b ð1Þ ðrðnÞ anþm Þ; P we will show that right side of above formula vanishes. Namely, using the properties (a), (b), (c), we obtain
b ðmþ1Þ P b ðmÞ P b ð1Þ ðrðnÞ anþm Þ ¼ ½ P b ð1Þ ðrðnÞEÞðanþm1 Þ b ð2Þ P b ðmþ1Þ P b ð2Þ P P ðbÞ
b ðmþ1Þ P b ð2Þ ðE P b ð1Þ DÞðrðnÞanþm1 Þ ¼ ½P b ðmþ1Þ P b ð3Þ E0 P b ð2Þ E P b ð2Þ DÞðrðnÞanþm1 Þ b ð4Þ P b ð1Þ Þ ð P b ðmþ1Þ P b ð3Þ P ¼ ½ð P
ða;cÞ
b ð3Þ E2 P b ð2Þ P b ð1Þ Þ ð P b ðmþ1Þ P b ð4Þ E1 P b ðmþ1Þ P b ð3Þ D P b ð1Þ ÞðrðnÞanþm1 Þ ¼ ½ð P
ða;cÞ
ða;cÞ b ð2Þ P b ðmÞ P b ð1Þ ÞðrðnÞanþm1 Þ b ðmþ1Þ Em P b ðmÞ P b ð1Þ Þ ðD P ¼ . . . ¼ ½ðE1m P b ðmþ1Þ Em DÞ P b ð1Þ ðanþm1 Þ: b ðmÞ P b ð2Þ P ¼ ðE1m P
Now, observe that
b ð2Þ P b ðmÞ P b ð2Þ P b ðmÞ P b ð1Þ ðanþm1 Þ ¼ P b ð1Þ ða1 þ a2 þ þ anþm1 Þ ¼ 0: P This completes the proof. h ðmÞ Theorem 12. Let the sequences bn ; W n ðxÞ, and the functions AnðmÞ ; BðmÞ n , and C n ; m 2 N, be defined as follows:
Case
I
II
III
bn
1 zn ðnþtÞv
1 zn
1 zn
W n ðxÞ
T n ðxÞ; U n ðxÞ
Pn ðxÞ
AðmÞ n BðmÞ n C ðmÞ n
nþmþtþv 2
nþm
LnðaÞ ðxÞ nþmþa
Here x; z; t 2 C and
2zxðn þ 2m þ t þ v 2Þ
2zxðn þ 2m 1=2Þ
2zðn þ 2m þ ða x 1Þ=2Þ
z2 ðn þ 3m þ t þ v 2Þ
z2 ðn þ 3m 1Þ
z2 ðn þ 3m 1Þ
v 2 N0 . Further, define the second-order operators
b ðmÞ :¼ AðmÞ I þ BðmÞ E þ C ðmÞ E2 ; P n n n
m 2 N:
ð5:18Þ
b ð1Þ ððn þ tÞ IÞ; in the case I; b ðm1Þ P b ðmÞ P P v 1 ðmÞ b ðm1Þ b ð1Þ ; b P in the cases II; III; P P
ð5:19Þ
Then the operators
( b L ðmÞ :¼ satisfy
b L ðmÞ ðan þ anþ1 þ þ anþm1 Þ ¼ 0;
m 2 N;
ð5:20Þ
with an ¼ bn W n ðxÞ. Proof. In all the cases, the proof follows by Lemma 11. Let us put in (5.14)
ðd; 1Þ :¼
ðv 1; ðtÞv 1 Þ; in the case I; ð0; 1Þ;
in the cases II; III;
where the number A, and consequently the numbers u; w; v; B; C, are such that definition (5.15) agrees with (5.18). Observe that
rðnÞ :¼
ðn þ tÞv 1 ; in the case I; 1;
in the cases II; III:
P. Woz´ny, R. Nowak / Applied Mathematics and Computation 215 (2009) 1622–1645
1638
Since
b ð1Þ ¼ ðjn I þ kn E þ l E2 Þðzn nðnÞ IÞ with nðnÞ :¼ zn P n
n þ t þ v 1; in the case I; 1;
in the cases II; III
(cf. (5.3)–(5.5)), one can see that assumption (5.16) is satisfied, too. h Example 5.1. Consider the series
2
p pffiffiffiffiffiffiffiffiffiffiffiffiffiffi
2 þ 2x ¼
2
1 X
bn T n ðxÞ;
bn :¼
n¼0
ð1Þn n2 1=4
ð5:21Þ
(cf. [21, (2), p. 143]). Depending on x 2 ½1; 1, its partial sums may behave very irregularly, so the classic methods of convergence acceleration usually fail. By Theorem 12 (cf. case I with z ¼ 1; t ¼ 1=2; v ¼ 2), we obtain the operators
b ðmÞ :¼ ðn þ m 1=2ÞI þ 2xðn þ 2m 1=2ÞE þ ðn þ 3m 1=2ÞE2 ; P satisfying
b ðmÞ P b ð1Þ ½ðn 1=2Þðbn T n ðxÞ þ bnþ1 T nþ1 ðxÞ þ þ bnþm1 T nþm1 ðxÞÞ ¼ 0: b ð2Þ P P
ð5:22Þ
Let us define the operators
b ð1Þ ½ðn 1=2ÞI; Pð1Þ :¼ P
b ðmÞ ; PðmÞ :¼ P
m P 2:
ð5:23Þ
For the summation of the series (5.21) one can use Wynn’s e-algorithm or, among others, the K transformation (5.11) with the parameters dðkÞ n :¼ 1=ðn þ 1Þ and the remainder estimates xn :¼ bn . However, the Algorithm 1 of performing the Q transformation, determined by operators (5.23), looks similar to the H transformation with parameters b ¼ 1=2; cos a ¼ x, and the remainder estimates xn :¼ ðn 1=2Þ2 ; see (5.9). In the sequel we test these methods, as well as some techniques given by Paszkowski [8]. A. Let x ¼ 1. Since T n ð1Þ ¼ ð1Þn , the partial sums of the series (5.21) are as follows
sn ¼
n1 X
bk T k ð1Þ ¼
n1 X
1
k¼0
k 1=4
2
k¼0
¼ 2
2 ; 2n 1
and hence s ¼ limn!1 sn ¼ 2. One can check that the Qð1Þ transformation, applied to the series (5.21), is exact, i.e.,
Qnð1Þ ¼ 2 for all n 2 N: The H and K transformations, as well as the variants u; v of Levin transformation, are exact, too. However, the series is convergent logarithmically:
accðs10 Þ ¼ 1:28;
accðs100 Þ ¼ 2:30;
accðs600 Þ ¼ 3:08;
accðs1000 Þ ¼ 3:30;
ð1Þ 30 Þ
and e-algorithm yields only accðe ¼ 2:70. B. Let x ¼ 0. The convergence of the series (5.21) is slow again, since we have
accðs10 Þ ¼ 2:85;
accðs100 Þ ¼ 4:92;
accðs400 Þ ¼ 6:13;
accðs1000 Þ ¼ 6:93:
Using Q transformation, one can compute some approximation of the sum of the series. Namely, using the operators (5.23) in ð15Þ ð15Þ Algorithm 1 yields accðQ1 Þ ¼ 23:18. Using Homeier’s transformations we obtain accðK1 Þ ¼ 17:07 and ð15Þ accðH1 Þ ¼ 20:37. Since T 2nþ1 ð0Þ ¼ 0 for n 2 N0 , the e-algorithm is unable to accelerate the convergence of the series; ð1Þ namely, we obtain only accðe30 Þ ¼ 3:28. But application of e-algorithm to the simplified form of the series (5.21), i.e.,
pffiffiffi
2
p 2 2
¼s¼
1 X
ð1Þn b2n ¼
n¼0
1 X n¼0
ð1Þn ; 4n2 1=4 ð1Þ
results in better acceleration of the convergence; namely we obtain accðe14 Þ ¼ 12:06. C. Now, let x ¼ 0:5. Then we obtain
accðs10 Þ ¼ 2:95;
accðs100 Þ ¼ 4:97;
accðs600 Þ ¼ 6:53;
accðs800 Þ ¼ 9:86;
accðs200 Þ ¼ 8:05;
accðs400 Þ ¼ 6:18;
accðs1000 Þ ¼ 6:97: ðkÞ
Using Algorithm 1 with operators (5.23), we obtain the following values of accðQ1 Þ for k ¼ 1; 3; . . . ; 15:
3:10;
6:38;
9:63;
13:18;
17:10;
20:64;
23:76;
27:26:
ð15Þ Q1
ð5:24Þ
Notice that computing requires the knowledge of partial sums s1 ; s2 ; . . . ; s31 . By Wynn’s e-algorithm the following values ð1Þ of quantities accðek Þ; k ¼ 2; 6; . . . ; 30, were obtained:
1:62;
4:10;
6:87;
8:31;
10:07;
12:17;
14:26;
16:01:
P. Woz´ny, R. Nowak / Applied Mathematics and Computation 215 (2009) 1622–1645
1639
Much worse results are obtained by using the variants t; u; v of Levin transformation, e.g., ð30Þ
accðt L1 Þ ¼ 4:05;
ð30Þ
accðu L1 Þ ¼ 3:94;
ð30Þ
accðv L1 Þ ¼ 2:14:
We also tried the method Wðr; xj Þ recently proposed by Paszkowski in [8]. We set the parameters r ¼ 0; xj ¼ bjþ1 , and obðkÞ tained the following accuracy of quantities s1 ; k ¼ 2; 6; . . . ; 30:
4:37;
9:00;
14:15;
17:84;
21:44;
24:90;
28:27;
31:58:
ð30Þ accðs1 Þ
ð5:25Þ ð15Þ accðQ1 Þ
D. Let x ¼ 1. Then using first 31 partial sums, we obtained ¼ 32:77 and ¼ 33:89. Thus our method is better than Paszkowski’s method. If we use first 50 (resp. 100) partial sums, then our method becomes better for x P 0:75 (resp. x P 0:45). We made numerical computations for x ¼ 0:95; 0:9; . . . ; 1:0 and compared the following methods of summation of the series (5.21), using first 50 partial sums (after colon, we give the notation of quantities tested): Usual summation: s50 . ð1Þ Wynn’s e-algorithm: e48 . ð24Þ Homeier’s Hða; b; fsn g; fxn gÞ transformation with cos a ¼ x; b ¼ 1=2 and xn ¼ ðn 1=2Þ2 : H1 . ð24Þ ðkÞ ðkÞ Homeier’s Kðfsn g; fxn g; fdn g; f1; 2x; 1gÞ transformation with dn ¼ 1=ðn þ 1Þ and xn ¼ bn : K1 . ð24Þ Q transformation, determined by operators (5.23): Q1 . ð50Þ Paszkowski’s method Wðr; xj Þ with r ¼ 0; xj ¼ bjþ1 : s1 . On Fig. 1, we present the accuracies obtained. We observed that e-algorithm fails for x ¼ 0, exceptionally. Computing of ð24Þ ð1Þ Q1 was performed by Algorithm 1 and lasted about 3:5 times longer than in the case of computing e24 . Let us notice that computations using Paszkowski’s method lasted about 30 times longer than e-algorithm. Remark 13. In Example 5.1 we considered the Chebyshev series with coefficients bn of the form
bn ¼
1 ; zn ðn þ tÞv
ð5:26Þ
where z ¼ 1; t ¼ 1=2; v ¼ 2. If we choose bigger value of parameter t, then Q transformation (determined by the operators corresponding to the case I of Theorem 12) can yield much better results than the method Wð0; bjþ1 Þ of Paszkowski [8]. For instance, using first 50 partial sums of the series
Fig. 1. Accuracy of the summation of the series (5.21) for x ¼ 0:95; 0:9; . . . ; 1.
P. Woz´ny, R. Nowak / Applied Mathematics and Computation 215 (2009) 1622–1645
1640 1 X n¼0
ð1Þn T n ðxÞ ðthe choice z ¼ 1; t ¼ 80; v ¼ 2Þ; ðn þ 80Þðn þ 81Þ
ð5:27Þ
our transformation is more efficient for x J 0:55; see Fig. 2. The choice t ¼ 130 makes our method better for all x ¼ 0:95; 0:9; . . . ; 1. Results presented on Fig. 2 correspond to the following methods:
Usual summation: s50 . ð1Þ Wynn’s e-algorithm: e48 . ð24Þ Homeier’s Hða; b; fsn g; fxn gÞ transformation with cos a ¼ x; b ¼ t and xn ¼ ðn þ tÞ2 : H1 . ð24Þ ðkÞ ðkÞ Homeier’s Kðfsn g; fxn g; fdn g; f1; 2x; 1gÞ transformation with dn ¼ 1=ðn þ 1Þ, and xn ¼ bn : K1 . ð24Þ Q transformation, determined by operators given by Theorem 12: Q1 . ð50Þ Paszkowski’s method Wðr; xj Þ with r ¼ 0, xj ¼ bjþ1 : s1 . Analysing Fig. 2, one can also observe that Hð24Þ and Qð24Þ transformations behave very similarly.
Example 5.2. Consider the following trigonometric series: 1 X n¼0
1 f zn ðn þ tÞv n
with f n ¼ sin nx; cos nx
ð5:28Þ
(cf. case I of Theorem 12). Combining (5.3) with (5.8), in the case of series (5.28), one can apply the Q transformation deterhas cos x instead of x. For mined by the operators corresponding to the case I of Theorem 12, where the definition of BðmÞ n instance, let us consider the series 1 X n¼0
cos nx : nþbþ1
ð5:29Þ
Observing that (5.18) are of the form
b ðmÞ ¼ ðn þ m þ bÞI 2 cos xðn þ 2m þ bÞE þ ðn þ 3m þ bÞE2 ; P
Fig. 2. Accuracy of the summation of the series (5.27) for x ¼ 0:95; 0:9; . . . ; 1.
P. Woz´ny, R. Nowak / Applied Mathematics and Computation 215 (2009) 1622–1645
1641
Table 5 Accuracy of the summation of the series (5.29) with b ¼ 10. 0:1p
x ð10Þ accðQ1 Þ ð10Þ accðH1 Þ ð21Þ accðs1 Þ ð21Þ accð~s1 Þ
0:2p
0:3p
0:4p
0:5p
0:6p
0:7p
0:8p
0:9p
5.81
10.28
13.37
15.41
18.15
19.44
22.63
23.71
26.39
p 27.43
4.91
9.44
11.88
14.29
16.50
18.65
20.93
24.44
26.65
27.48
12.00
18.18
21.74
24.88
26.96
29.25
30.02
32.43
31.36
31.43
12.13
19.16
22.18
24.81
27.76
28.64
30.65
32.35
30.92
30.84
we obtain the following form of Algorithm 1:
Nð0Þ n ¼ sn ;
Dð0Þ n ¼ 1; ðk1Þ
ðk1Þ
ðk1Þ N ðkÞ 2 cos x ðn þ 2k þ bÞNnþ1 þ ðn þ 3k þ bÞNnþ2 ; n ¼ ðn þ k þ bÞN n ðk1Þ DðkÞ 2 cos x ðn þ 2k þ n ¼ ðn þ k þ bÞDn
ðk1Þ bÞDnþ1
þ ðn þ 3k þ
ðk1Þ bÞDnþ2 ;
ð5:30Þ
ðkÞ ðkÞ QðkÞ n :¼ N n =Dn ;
which is similar to the algorithm (5.9) of performing HðkÞ n ðx; b; fsn g; fxn gÞ. We applied the following methods to the case of summation of the series (5.29): (a) Q transformation given by (5.30), (b) transformation HðkÞ n ðx; b; fsn g; fxn gÞ of Homeier [10], (c) method Wðr; xj Þ of Paszkowski [8]. In the case of Homeier’s transformation, considering classic variants of remainder estimates xn , does not accelerate the convergence, at all. We choose remainder estimates:
xn :¼
1 : nþbþ1
We made a numerical computations for b ¼ 10 and b ¼ 100. In the case of Paszkowski’s method, we put r ¼ 0, xj ¼ 1=ðj þ b þ 1Þ (cf. [8, Case (1), p. 41]) and r ¼ 1; xj 1 (cf. [8, Case (2), p. 41]). The arrays obtained we denote by sðkÞ n ðkÞ and ~sn , respectively. Putting x ¼ 0:1 p; 0:2 p; . . . ; p, we obtain the accuracies presented in the Table 5 ðb ¼ 10Þ and Table 6 ðb ¼ 100Þ. The input for all algorithms was the set of the partial sums s1 ; s2 ; . . . ; s21 of the series (5.29). One can observe that Paszkowski’s method results in the best convergence acceleration, in case b ¼ 10. In each case, the transformation (5.30) gives accuracies growing with x and permits us to obtain similar accuracy as Homeier’s does. In the case b ¼ 100, both methods yield much higher accuracy than Paszkowski’s Wðr; xj Þ method. 6. Iterated methods Knowing the explicit form of operators PðmÞ satisfying (2.2), allows us to apply the transformation
QðmÞ n ðsn Þ ¼
LðmÞ ðsn Þ ; LðmÞ ð1Þ
LðmÞ ¼ PðmÞ Pðm1Þ Pð1Þ ;
ð6:1Þ
in a very efficient way, using the Algorithm 1. Usually, it is difficult to obtain the operators PðmÞ ; m 2 N, explicitly. In most cases, it is enough to obtain the operators Pð1Þ ; Pð2Þ ; . . . ; Pðm0 Þ , m0 being a small number. But if it is still hard task, then we propose to use one of the following two iterated methods. I. Let us consider the classic idea to iterate the Qð1Þ transformation. We define the U transformation by
Uð0Þ n :¼ sn ;
UðmÞ :¼ n
Pð1Þ ðUðm1Þ Þ n ; ð1Þ P ð1Þ
m 2 N:
ð6:2Þ
Notice that such classic methods as iterated D2 Aitken’s process, or Brezinski’s #-algorithm [1, (10.1-9)], were obtained using such idea. II. We propose the following simplification: establish Pð1Þ and put
PðmÞ :¼ Pð1Þ ;
m ¼ 2; 3 . . . :
Then the transformation (6.1) simplifies to the V transformation, defined by
VðmÞ :¼ n
½Pð1Þ m ðsn Þ : ½Pð1Þ m ð1Þ
ð6:3Þ
Notice that both iterated methods require only the operator Pð1Þ . Since Pð1Þ ¼ Lð1Þ , it may be obtained as in the Section 2.1.
P. Woz´ny, R. Nowak / Applied Mathematics and Computation 215 (2009) 1622–1645
1642
Table 6 Accuracy of the summation of the series (5.29) with b ¼ 100. x
0:1p
0:2p
0:3p
0:4p
0:5p
0:6p
0:7p
0:8p
0:9p
ð10Þ accðQ1 Þ ð10Þ accðH1 Þ ð21Þ accðs1 Þ ð21Þ accð~s1 Þ
15.95
21.18
24.66
28.54
29.67
32.07
33.73
36.31
39.23
p 46.15
14.65
21.10
25.26
27.63
29.91
32.09
34.33
36.92
40.77
47.69
10.02
13.78
16.26
18.10
19.38
20.35
21.21
22.52
22.11
22.02
7.40
10.62
12.63
14.22
15.36
16.24
17.04
18.51
17.80
17.71
The following example shows that V transformation can yield not much worse results than original Q transformation (with proper operators PðmÞ ). Of course, such behavior is quite unpredictable and, in general, it is hard to say if iterated transformation can accelerate the convergence of considered series. The following example shows that the V transformation performs better in some cases, while U – in other cases. Next we give a remark considering the case of summation of some orthogonal polynomial series given by Paszkowski [8]. They do not fit to any case of Theorem 12. In such a case, numerical tests show that the Paszkowski’s method is much more efficient than the iterated methods U and V. In the Example 6.2, we present the series which is neither orthogonal nor hypergeometric. We make a special choice of operator Pð1Þ , and obtain much better results by using the second iterated method V. Example 6.1. Consider the orthogonal polynomial series (5.21) and suppose we do not know the operators (5.23). Since
ðI 2xE þ E2 ÞðT n ðxÞÞ ¼ 0 (cf. (5.3)), one can put
b ð1Þ :¼ ðI 2xE þ E2 Þ P
1 I ; bn
ð6:4Þ
b ð1Þ ðbn T n ðxÞÞ ¼ 0. Now, one can use the both iterated methods with Pð1Þ ¼ P b ð1Þ . Observe that one can also consider and obtain P ð1Þ the operator P of the form
b ð1Þ Pð1Þ ¼ Q P
ð6:5Þ
with arbitrary operator Q. Let us put x ¼ 0:5 (cf. part C of Example 5.1). We put
Q¼
ð1Þn I ð2n 1Þð2n þ 1Þ
in (6.5), and obtained
Pð1Þ ¼ I þ
2n þ 3 ð2n þ 3Þð2n þ 5Þ 2 Eþ E : 2n 1 ð2n 1Þð2n þ 1Þ
ð6:6aÞ ð15Þ
ð15Þ
We computed the transformations (6.2) and (6.3), obtaining accðU1 Þ ¼ 15:03; accðV1 Þ ¼ 12:55. One can also consider
Q¼
ð1Þn1 E1 ; ð2n 3Þð2n 1Þ
Pð1Þ ¼ E1 þ
2n þ 1 ð2n þ 1Þð2n þ 3Þ Iþ E; 2n 3 ð2n 3Þð2n 1Þ
ð6:6bÞ
Q¼
ð1Þn E2 ; ð2n 5Þð2n 3Þ
Pð1Þ ¼ E2 þ
2n 1 1 ð2n 1Þð2n þ 1Þ E þ I: 2n 5 ð2n 5Þð2n 3Þ
ð6:6cÞ
or
However, the best results are obtained for the choice (6.6a); see Table 7.
Table 7 Accuracy of the summation of the series (5.21) for x ¼ 0:5, obtained by U; V transformations, determined by operators Pð1Þ of the form (6.6c) and (6.7c). Pð1Þ
n
accðUn
(6.6a) (6.6b) (6.6c)
1 15 31
15.03 13.87 11.48
12.55 9.89 9.62
(6.7a) (6.7b) (6.7c)
1 15 31
15.03 13.87 11.48
21.86 25.26 26.87
ð15Þ
Þ
ð15Þ
accðVn
Þ
P. Woz´ny, R. Nowak / Applied Mathematics and Computation 215 (2009) 1622–1645
1643
We considered also the operators Q of the form
Q 2 fð1Þn I; ð1Þn1 E1 ; ð1Þn E2 g; and get the following variants of the operator Pð1Þ with polynomial coefficients:
Pð1Þ ¼ ð2n 1Þð2n þ 1ÞI þ ð2n þ 1Þð2n þ 3ÞE þ ð2n þ 3Þð2n þ 5ÞE2 ;
ð6:7aÞ
Pð1Þ ¼ ð2n 3Þð2n 1ÞE1 þ ð2n 1Þð2n þ 1ÞI þ ð2n þ 1Þð2n þ 3ÞE;
ð6:7bÞ
Pð1Þ ¼ ð2n 5Þð2n 3ÞE2 þ ð2n 3Þð2n 1ÞE1 þ ð2n 1Þð2n þ 1ÞI;
ð6:7cÞ
respectively. The obtained accuracies are shown in Table 7. One can observe that, in the case of using operators (6.6), U transformation yields better results than V. But when operators (6.7) are used, the V transformation gives essentially better results. The normalization of operators (6.6) into (6.7) does matter in the case of V transformation, and does not in the case of U transformation. Remark 14. In the case of orthogonal series, the methods given by Paszkowski [8] are more general than our transformation, and can be successfully applied in the case of many important orthogonal series of the form (5.1). We are not able to obtain explicit forms of operators PðmÞ in the case of Legendre series: 1 X
1 P n ðxÞ; nþ1
1 X
1 P n ðxÞ; ð2n 1Þð2n þ 3Þ n¼1 n¼1 1 1 X X logð1 þ 2n1 ÞPn ðxÞ; ðDÞ ¼ ð2n þ 1Þzn Pn ðxÞ ðCÞ ¼
ðAÞ ¼
ðBÞ ¼
n¼1
n¼1
(see [8, (10–13)]). However, we can apply the iterated methods, since for any orthogonal series (5.1), one can choose the following operator Pð1Þ :
1 Pð1Þ :¼ ðjn I þ kn E þ ln E2 Þ I : bn
ð6:8Þ
We observed that better results can be obtained by using the transformation
RðmÞ :¼ n
LðmÞ ðsn Þ ; LðmÞ ð1Þ
LðmÞ :¼ PðmÞ Pðm1Þ Pð1Þ ;
ð6:9Þ
determined by operators
b ð1Þ 1 I ; Pð1Þ :¼ P bn
b ðmÞ ; m P 2; PðmÞ :¼ P
ð6:10Þ
b ðmÞ ; m P 1, are given by (5.15) and u; w; v; A; B; C are such that Pð1Þ fits to (6.8). Notice that this is possible where operators P in each case (5.3)–(5.6) since jn ; kn ; ln are linear polynomials in n. We introduced the notation RðmÞ since the operators LðmÞ , given in (6.9), do not satisfy (2.2), as it is required in Q transformation. One can check that in the case of Legendre polynomial series, the R transformation is determined by operators (6.10), where
b ðmÞ ¼ ðn þ mÞI 2xðn þ 2m 1=2ÞE þ ðn þ 3m 1ÞE2 P (cf. (5.4)). For instance, let x ¼ 0:4 and z ¼ 0:8. We tested the summation of the series (A)–(D) by both iterated methods, determined by operator (6.8), and by the R transformation. We also tested Homeier’s K transformation. Results are presented in the Table 8. The last row corresponds to the case Wð0; bjþ1 Þ of Paszkowski’s method [8], where bn are the coefficients of P n ðxÞ in (A)–(D). In the case of K transformation, we put the parameters dðkÞ n ¼ 1=ðn þ 1Þ; xn ¼ bn (cf. (5.11)). Both iterated methods and R transformation perform much worse than Paszkowski’s method. It may be caused by wrong choice of operators PðmÞ , in the sense that the operators (6.10) do not satisfy the assumption in Algorithm 1, in the case of the series (A)–(D). More precisely, (2.2) holds only for m ¼ 1. In the case (D), we observed that Homeier’s transformation performs better than R transformation if z is closer and closer 1. In contrast with our transformation, Paszkowski’s method can be successfully applied also for x close to 1.
Table 8 ð20Þ ð20Þ ð20Þ ð41Þ Accuracies of U1 ; V1 ; R1 and s1 , in the case of summation of the series (A)–(D) with x ¼ 0:4 and z ¼ 0:8. (A) ð20Þ accðU1 Þ ð20Þ accðV1 Þ ð20Þ accðR1 Þ ð20Þ accðK1 Þ ð41Þ accðs1 Þ
(B)
(C)
(D)
8.68
9.96
7.39
3.30
20.12
11.51
18.86
6.78
22.18
20.95
18.79
19.52
21.13
21.74
19.89
24.42
47.82
34.38
32.54
33.06
P. Woz´ny, R. Nowak / Applied Mathematics and Computation 215 (2009) 1622–1645
1644
Example 6.2. Consider the series 1 X
an ;
an :¼
n¼0
nn pffiffiffiffiffiffiffiffiffiffiffiffi ðn þ 2Þ n þ 1
ðn ¼ 1Þ:
ð6:11Þ
A. Let n ¼ 1. Then (6.11) is the alternating series, and its convergence can be accelerated by classic method as Wynn’s ealgorithm, or Aitken’s iterated D2 process, very well. For instance, using first 21 partial sums, we obtain ð1Þ
ð1Þ
accðe20 Þ ¼ 16:0;
accðA10 Þ ¼ 20:1:
Consider the operator
1 Pð1Þ ¼ Lð1Þ :¼ D I an
ð6:12Þ
(cf. (2.4)). We are not able to find explicit forms of operators PðmÞ given by the method described in Section 2.1. But even if we do, results obtained by Q transformation will be the same as the ones given by e-algorithm: ð20Þ
ð1Þ
accðQ1 Þ ¼ accðe40 Þ ¼ 31:30; see Section 2.3. Using the iterated methods determined by operator (6.12) yields ð20Þ
accðU1 Þ ¼ 13:1;
ð20Þ
accðV1 Þ ¼ 5:3:
The efficiency of the U transformation is higher probably because it is based on an idea similar to the one used in the Aitken’s process. If we use D2 instead of D in (6.12), then we obtain the second-order operator Pð1Þ which results in a little bit better ð10Þ ð10Þ efficiency of both iterated methods: accðU1 Þ ¼ 13:9, accðV1 Þ ¼ 8:3. We observed that much better results one can obtain by using V transformation, determined by the following second-order operator Pð1Þ :
pffiffiffiffiffiffiffiffiffiffiffiffi Pð1Þ :¼ Dððn þ 2Þðn þ 3ÞDÞ n þ 1I ;
ð6:13Þ
where D :¼ E nI. One can check that results obtained by the first iterated method remain the same. Namely, we obtain ð10Þ
accðU1 Þ ¼ 13:9;
ð10Þ
accðV1 Þ ¼ 22:6:
B. Now, let n ¼ 1. Then the series is convergent logarithmically, and thus very slowly. Wynn’s and Aitken’s methods do not accelerate its convergence, at all: ð1Þ
accðe20 Þ ¼ 1:1;
ð1Þ
accðA10 Þ ¼ 1:7:
Using the iterated methods determined by operator Pð1Þ of the form (6.12), or with D2 instead of D, does not result in better ð10Þ accuracies. However, the second iterated method determined by operator (6.13) yields accðV1 Þ ¼ 8:3. Let us notice that usual summation of the first million terms of the series (6.11) yields less than three exact significant decimal digits of its limit. Acknowledgements The authors wish to express their thanks to Prof. E.J. Weniger for his interest in this research, and valuable comments and suggestions. They also wish to thank their supervisor Prof. S. Lewanowicz for his pertinent comments and helpful discussion. Finally, they are grateful to the anonymous referees for their suggestions and remarks which helped to improve the paper. References [1] E.J. Weniger, Nonlinear sequence transformations for the acceleration of convergence and the summation of divergent series, Comput. Phys. Rep. 10 (1989) 189–371. [2] J. Cˇízˇek, J. Zamastil, L. Skála, New summation technique for rapidly divergent perturbation series. Hydrogen atom in magnetic field, J. Math. Phys. 44 (3) (2003) 962–968. [3] C. Brezinski, M.R. Zaglia (Eds.), Extrapolation methods: theory and practice, Studies in Computational Mathematics, vol. 2, North-Holland, 1991. [4] H.H. Homeier, Scalar Levin-type sequence transformations, J. Comput. Appl. Math. 122 (1–2) (2000) 81–147. [5] J. Wimp, Sequence transformations and their applications, Mathematics in Science and Engineering, vol. 154, Academic Press, New York – London, 1981. [6] C. Brezinski, A.C. Matos, A derivation of extrapolation algorithms based on error estimates, J. Comput. Appl. Math. 66 (1–2) (1996) 5–26. [7] A.C. Matos, Linear difference operators and acceleration methods, IMA J. Numer. Anal. 20 (3) (2000) 359–388. [8] S. Paszkowski, Convergence acceleration of orthogonal series, Numer. Algorithms 47 (1) (2008) 35–62. ˇ ízˇek, Zamastil, Skála. I. Algebraic theory, J. Math. [9] E.J. Weniger, Mathematical properties of a new Levin-type sequence transformation introduced by C Phys. 45 (3) (2004) 1209–1246. [10] H.H.H. Homeier, A Levin-type algorithm for accelerating the convergence of Fourier series, Numer. Algorithms 3 (1992) 245–254. [11] H.H.H. Homeier, Some applications of nonlinear convergence accelerators, Int. J. Quantum Chem. 45 (6) (1993) 545–562. [12] H.H.H. Homeier, Nonlinear convergence acceleration for orthogonal series, in: R. Gruber, M. Tomassini (Eds.), Proceedings of the 6th Joint EPS–APS International Conference on Physics Computing, Physics Computing’94, European Physical Society, Boite Postale 69, CH-1213 Petit-Lancy, Geneva, Switzerland, 1994, pp. 47–50. [13] C. Brezinski, M. Redivo-Zaglia, On the kernel of sequence transformations, Appl. Numer. Math. 16 (1994) 239–244.
P. Woz´ny, R. Nowak / Applied Mathematics and Computation 215 (2009) 1622–1645 [14] [15] [16] [17] [18] [19] [20] [21]
1645
D. Shanks, Nonlinear transformations of divergent and slowly convergent sequences, J. Math. Phys. 34 (1955) 1–42. G.E. Andrews, R. Askey, R. Roy, Special Functions, Cambridge University Press, Cambridge, 1999. G.M. Phillips, Interpolation and approximation by polynomials, CMS Books in Mathematics, vol. 14, Springer, New York, 2003. S. Lewanowicz, S. Paszkowski, An analytic method for convergence acceleration of certain hypergeometric series, Math. Comput. 64 (210) (1995) 691– 713. G. Gasper, M. Rahman, Basic hypergeometric series, Encyclopedia of Mathematics and its Applications, second ed., vol. 96, Cambridge University Press, 2004. R. Koekoek, R.F. Swarttouw, The Askey-scheme of Hypergeometric Orthogonal Polynomials and its q-Analogue, Technical Report, Delft University of Technology, 1998. C. Oleksy, A convergence acceleration method of fourier series, Comput. Phys. Commun. 96 (1) (1996) 17–26. S. Paszkowski, Numerical applications of Chebyshev polynomials and series (Zastosowania numeryczne wielomianow i szeregow Czebyszewa.), PWN, Warszawa, 1975.