- Email: [email protected]
A generalization of the G-transformation and the related algorithms
Applied Numerical Mathematics 60 (2010) 1221–1230
Contents lists available at ScienceDirect
Applied Numerical Mathematics www.elsevier.com/locate/apnum
A generalization of the G-transformation and the related algorithms Claude Brezinski a , Yi He b,c , Xing-Biao Hu b,∗ , Jian-Qing Sun b,c a
Laboratoire Paul Painlevé, UMR CNRS 8524, UFR de Mathématiques Pures et Appliquées, Université des Sciences et Technologies de Lille, 59655-Villeneuve d’Ascq cedex, France LSEC, Institute of Computational Mathematics and Scientific, Engineering Computing, AMSS, Chinese Academy of Sciences, P.O. Box 2719, Beijing 100190, PR China c Graduate School of the Chinese Academy of Sciences, Beijing, PR China b
a r t i c l e
i n f o
Article history: Received 30 November 2009 Received in revised form 4 March 2010 Accepted 10 March 2010 Available online 18 March 2010 Keywords: G-transformation G-algorithm rs-algorithm qd-algorithm
a b s t r a c t In this paper, a generalization of the G-transformation is proposed together with the corresponding auxiliary rs-algorithm for implementation. We show that the related generalization of the rs-algorithm is equivalent to a generalized qd-algorithm. Applications of the generalization of the G-transformation to the computation of infinite integrals are also given. © 2010 IMACS. Published by Elsevier B.V. All rights reserved.
1. Introduction When a sequence ( S n ) is slowly converging to its limit S, it can be transformed, without altering its terms, into a new sequence ( T n ) which, under some assumptions, converges faster to the same limit, that is such that
lim
n→∞
Tn − S Sn − S
= 0.
Sequence transformations are based on the idea of extrapolation. There exist many such sequence transformations since a general transformation able to accelerate the convergence of all sequences cannot exist as proved in [3]. Each sequence transformation can only accelerate the convergence of some classes of sequences. In many cases, the terms of the sequence ( T n ), which also often depend on a second index denoted by k, are expressed as ratios of determinants, thus needing the use of a recursive algorithm for the implementation of the transformation. On sequence transformations and the related algorithms, consult [2]. In [4], the authors proposed a nonlinear sequence transformation called the G-transformation, and applied it to the computation of improper integrals. In [5], a higher order G-transformation was introduced, which was later proved to be a particular case of the E-transformation [1], the most general sequence transformation known so far. These transformations are implemented by the G-algorithm and the E-algorithm, respectively. The former was proposed, together with the rs-algorithm as its auxiliary algorithm, for implementing this higher transformation [8], while the latter was presented with the help of another set of recursion relations. As far as the G-transformation is concerned, the auxiliary rs-algorithm is less costly than the auxiliary algorithm for the E-algorithm.
*
Corresponding author. E-mail addresses: [email protected] (C. Brezinski), [email protected] (Y. He), [email protected] (X.-B. Hu), [email protected] (J.-Q. Sun).
0168-9274/$30.00 © 2010 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2010.03.008
1222
C. Brezinski et al. / Applied Numerical Mathematics 60 (2010) 1221–1230
It has been shown that the rs-algorithm is equivalent to the qd-algorithm [9,2]. The ε -algorithm [12] is a well-known extrapolation method, and two generalizations of this algorithm were proposed and studied in [2,11]. They all are particular cases of the G-transformation. Similarly to the generalizations of the ε -algorithm, we will now study a generalization of the G-transformation. By using determinantal identities, we will obtain a generalized G-algorithm and the related rs-algorithm for implementing the transformation. The generalized rs-algorithm reduces to the original rs-algorithm for a particular choice of the auxiliary sequence involved. The generalized rs-algorithm is equivalent to a generalized qd-algorithm. Imposing boundary value conditions on the generalized qd-algorithm, we will obtain a generalized qd-equation with a Lax representation. There are many connections between extrapolation algorithms and discrete integrable systems which could be considered as discretization of partial differential equations. An important property of such systems is their integrability, that is the possibility of obtaining their solution into a close form such as, for example, ratios of determinants. This is, in particular, the case of extrapolation algorithms such as the ε -algorithm and its generalizations, and the G-transformation which have a determinantal representation. Thus, the results presented into this paper are motivated by this connection. This paper is organized as follows. In Section 2, we give a generalization of the G-transformation, and the related G-algorithm and rs-algorithm. In Section 3, we obtain a generalized qd-algorithm, and show that it is equivalent to the generalized rs-algorithm related to the generalization of the G-algorithm. In Section 4, applications of the generalization of the G-transformation to the computation of improper integrals are given. Conclusions and a discussion are presented in Section 5. 2. A generalization of the G -transformation In this section, we give a generalization of the G-transformation. Then, we derive the generalized G-algorithm to implement it on the basis of a generalized rs-algorithm. This generalized G-transformation has the form
un v n .. . k −1 R vn (n) Gk = yn v n .. . k −1 R
vn
u n +1 v n +1
.. .
R k −1 v n +1 y n +1 v n +1
.. .
R k −1 v n +1
.. . k −1 ··· R v n+k , ··· yn+k ··· v n+k .. . k −1 ··· R v n+k ··· ···
un+k v n+k
where (xn ) and ( yn ) are auxiliary sequences, un = S n yn , v n = xn yn , and the operator R is defined by R v n = ( v n / yn ) and R k+1 v n = ( R k v n / yn ). Obviously, the condition yn = 0 for all n has to be satisfied. When ∀n, yn = 1, the general Gtransformation reduces to the original one. It must be noticed that this transformation enters into the general framework of the E-algorithm [1] by dividing the first columns in the numerator and the denominator by yn , the second ones by yn+1 , and so on. Set
vn Rv n .. . k −1 R vn (n) rk = yn v n .. . k −2 R
(n)
sk =
vn
yn vn R
vn Rv n .. . k −1 R
.. .
R k −1 v n +1 y n +1 v n +1
.. .
R k −2 v n +1 y n +1 v n +1
.. .
k −1
v n +1 R v n +1
.. .
vn
vn
k −1
R v n +1 v n +1 R v n +1
.. .
R k −1 v n +1
··· ···
v n+k−1 R v n+k−1
k −1 ··· R v n+k−1 , ··· yn+k−1 ··· v n+k−1 .. . · · · R k−2 v n+k−1 ··· yn+k ··· v n+k .. . k −1 ··· R v n+k . ··· v n+k−1 ··· R v n+k−1 .. . k −1 ··· R v n+k−1
Then, a generalized rs-algorithm holds
.. .
C. Brezinski et al. / Applied Numerical Mathematics 60 (2010) 1221–1230
1223
Theorem 1. (n) (n) The quantities rk and sk defined above can be recursively computed by the algorithm (n)
(n)
s 0 = 1,
r1 = xn ,
(n)
(n+1)
(n) r k +1
(n+1)
sk+1 = sk
= rk
(n+1) r k +1 (n) r k +1
n = 0, 1, . . . ,
−1 ,
k = 0, 1 , . . . , n = 0, 1 , . . . ,
(n+1)
sk
1
−
yn+k+1 s(n) k
1 yn
,
k = 1 , 2 , . . . , n = 0, 1 , . . . .
(n)
Proof. From the expressions of the quantities rk
(n)
sk+1
(n+1)
sk
y n vn . .. k R vn = vn R vn . .. k
··· ···
y n +1 v n +1
R v n +1
.. . R k v n +1
R vn
(n)
v n+k+1
.. .
R k v n +1 v n +1
(2)
and sk , we have
yn+k+1
.. .
(1)
· · · R k v n+k+1 · · · v n+k · · · R v n+k .. . k · · · R v n+k
v n +1 R v n +1 .. . k −1
··· ···
v n +2 R v n +2
v n+k R v n+k
.. .
.. .
R v n +1 y n +1
R k −1 v n +2 · · · y n +2 ···
v n +1
v n +2
.. .
.. .
R k −1 v n +1
R k −1 v n +2
R k−1 v n+k yn+k+1
···
v n+k+1
.. . · · · R k−1 v n+k+1
Applying Sylvester’s determinantal identity, we get
y n +1 ··· yn+k+1 v n+1 v n +2 ··· yn v v n +1 · · · v n+k+1 R v n+1 R v n +2 ··· n . .. .. .. .. . . . . . . k R v n R k v n+1 · · · R k v n+k+1 R k−1 v n+1 R k−1 v n+2 · · · v n +2 · · · v n+k+1 yn y n +1 v n +1 Rv R v n+2 · · · R v n+k+1 v n v n +1 n +1 = .. .. .. .. .. . . . . . k R v n+1 R k v n+2 · · · R k v n+k+1 R k−1 v n R k−1 v n+1 v n +1 · · · v n+k yn+1 y n +2 vn Rv R v n+1 · · · R v n+k v n+1 v n +2 n − . .. .. .. .. .. . . . . k R v n R k v n+1 · · · R k v n+k R k−1 v n+1 R k−1 v n+2
v n+k R v n+k
.. .
R k−1 v n+k ··· yn+k ··· v n+k
.. .
· · · R k−1 v n+k ··· yn+k+1 ··· v n+k+1 . .. . k −1 ··· R v n+k+1
Then, the relation (1) is proved. Let us now turn to the proof of the second relation. We have
(n) r k +1 (n+1)
rk
=
vn Rv n . . . k
R vn yn vn
.. .
R k −1 v n
v n +1 R v n +1
··· ···
R k v n +1 y n +1 v n +1
· · · R k v n+k R k−2 v n+1 R k−2 v n+2 · · · ··· yn+k v n+1 v n +2 ··· v n+k R v n+1 R v n +2 .. .. .. . . . · · · R k−1 v n+k R k−1 v n+1 R k−1 v n+2
.. .
.. .
R k −1 v n +1
v n+k R v n+k
.. .
y n +1 v n +1
.. .
y n +2 v n +2
.. .
··· ···
yn+k v n+k
.. .
R k−2 v n+k ··· v n+k ··· R v n+k
.. .
· · · R k−1 v n+k
.
By definition of the operator R, it holds
y n +1 v n +1
y n +2 v n +2
R k −2 v n +1
R k −2 v n +2
.. .
.. .
··· ···
yn+k v n+k
.. .
· · · R k−2 v n+k
··· R v n+k−1 R v n +1 .. .. = yn+1 yn+2 · · · yn+k . . . k −1 k −1 R v n +1 · · · R v n+k−1
.
1224
C. Brezinski et al. / Applied Numerical Mathematics 60 (2010) 1221–1230
Then, Sylvester’s identity gives
vn Rv n . . . k
··· ···
v n +1 R v n +1
.. .
v n+k R v n+k
.. .
y n +1 v n +1
.. .
y n +2 v n +2
··· ···
.. .
yn+k v n+k
.. .
· · · R k v n+k R k−2 v n+1 R k−2 v n+2 · · · R k−2 v n+k v n +1 · · · v n+k vn ··· R v n+k−1 Rv R v n +1 R v n+1 · · · R v n+k n .. .. = yn+1 yn+2 · · · yn+k . . . . . .. .. .. k −1 k −1 k R v n +1 · · · R v n+k−1 k k R v n R v n+1 · · · R v n+k ⎛ v n +1 ··· v n+k−1 vn R v n+1 · · · R v n+k R v R v n +1 ··· R v n+k−1 ⎜ n .. .. = yn+1 yn+2 · · · yn+k ⎜ .. .. .. . . ⎝ k . . . k R v n+1 · · · R v n+k k−1 R v n R k−1 v n+1 · · · R k−1 v n+k−1 ⎞ v n +2 ··· v n+k v n +1 R v n · · · R v n+k−1 R v R v n +2 · · · R v n+k ⎟ . n + 1 . ⎟ − .. .. .. .. ⎠ .. k . . . k R v n · · · R v n+k−1 k−1 R v n+1 R k−1 v n+2 · · · R k−1 v n+k y n +2 ··· yn+k+1 v n v n +1 ··· v n+k−1 y n +1 v v n +2 ··· v n+k+1 R v n R v n +1 ··· R v n+k−1 n +1 1 = .. .. .. .. .. .. yn+k+1 . . . . . . k −1 k −1 k −1 k −1 k −1 k −1 R v n +1 R v n +2 · · · R v n+k+1 R vn R v n +1 · · · R v n+k−1 y n +1 ··· yn+k v n+1 v n +2 ··· v n+k yn v v n +1 ··· v n+k R v n+1 R v n +2 ··· R v n+k n 1 . − . . . . . . .. .. .. .. .. .. yn k −1 k −1 k −1 k −1 k −1 k −1 R vn R v n +1 · · · R v n+k R v n +1 R v n +2 · · · R v n+k
R vn
R k v n +1
From the above relations, we obtain (2), which ends the proof of the theorem.
2
(n)
The quantities G k+1 can be recursively computed by (n)
G 0 = Sn, (n)
(n)
G k = G k −1 −
(3) (n+1) (n) G k−1 − G k−1 (n) r . (n+1) (n) k
rk
− rk
(4)
The kernel of the transformation is given by the following theorem Theorem 2. (n) A necessary and sufficient condition that ∀n, G k = S is that, ∀n,
( S n − S ) y n = b 0 v n + b 1 R v n + · · · + b k −1 R k −1 v n . 3. A generalized qd-algorithm In this section, we present a generalized qd-algorithm related to the generalized rs-algorithm (1)–(2). Set (n)
qk = (n)
ek =
(n+1) (n+1) sk−1 , (n) (n) rk sk−1 (n) (n) rk+1 sk . (n+1) (n+1) rk sk−1
rk
Then, the generalized qd-algorithm is as follows
(5)
(6)
C. Brezinski et al. / Applied Numerical Mathematics 60 (2010) 1221–1230
1225
Theorem 3. (n) (n) The qk and ek defined above can be computed by the generalized qd-algorithm (n)
(n)
e 0 = 0, 1
(n)
ek +
q1 = xn+1 /xn ,
yn+k
(n) (n)
1
(n)
qk −
yn
(7)
(n+1) = e k −1
1
+
(n+1)
yn+k+1
1
−
qk
y n +1
(8)
,
(n−1) (n−1) qk+1 .
qk e k = e k
(9)
Proof. Substituting (5) and (6) into both sides of (9) immediately shows that this formula holds. (n) (n) In order to prove formula (8), we first rewrite qk and ek using (1)–(2) as (n)
qk = (n)
ek =
(n+1) (n+1) sk−1 (n) (n) rk sk−1 (n) rk+1 sk(n) (n+1) (n+1) rk sk−1
rk
(n+1)
=
rk
r k +1
(n+1)
r k −1
yn
(n+1)
(n+1)
r yn+k rk yn+k = yn+k k(n+1) + , yn r (n) r k −1 k
(n) (n) r r 1 − 1 = k(+n)1 − (nk+ +1 )
(n+1)
rk
(n)
rk
1
+
(n+1)
(n)
rk
(n)
=
(n)
rk
rk
rk
rk
and (n)
qk = (n)
ek =
(n+1) (n+1) sk−1 (n) (n) rk sk−1 (n) rk+1 sk(n) (n+1) (n+1) rk sk−1
rk
=
(n)
sk
(n+1)
+1
1
sk
sk−1
=
(n+1)
sk−1 (n)
sk−1
(n+1)
1
−
yn+k+1 s(n) k
(n+1)
(n)
sk
=
(n)
sk−1
+
,
(n)
sk−1
(n)
sk
(n+1)
yn
sk−1
sk−1
=
(n+1)
sk
1
yn+k+1 s(n+1) k −1
−
1
(n)
sk
. yn s(n+1) k −1
Then, for the left-hand side of (8), we get
LHS of (8) =
1
yn+k
(n+1)
rk
yn+k (n+1) + · (n) yn r k −1 rk
(n+1)
=
rk
(n+1)
r k −1
+
1
rk
(n+1) r k −1
+
(n+1)
sk
yn+k+1 s(n+1) k −1
(n+1)
=
(n+1)
rk
yn+k
1
(n+1)
sk
1
+
+
rk
(n+1)
sk
1
1
(n)
sk
− yn+k+1 s(n+1) yn s(n+1) k −1 k −1
(n+1)
rk
(n)
yn
yn+k+1 s(n+1) k −1
(n)
−1−
sk
−
1 yn
(n+1)
sk−1
by (1) .
For the right-hand side of (8), we have (n+1)
(n+1)
r r 1 RHS of (8) = k(n+1) − k(n+2) + y n + k +1 r k −1 r k −1 (n+1)
=
rk
(n+1)
r k −1
+
rk
(n+1) r k −1
+
(n+1)
1
sk
1
sk
+ yn+k+1 s(n+1) k −1
(n+1)
=
(n+1)
yn+k+1 s(n+1) k −1
(n+2)
(n+1)
sk
(n+1) sk−1
+
sk−1
−
1
(n+1) y n +1 sk−1 (n+2) (n+1) sk−1 rk 1 1 − − ( n + 1 ) ( n + 2 ) yn+k+1 s y n +1 r k −1 k −1
by (2)
= LHS of (8), which completes our proof.
2 (n)
(n)
If we impose the boundary condition e 0 ≡ em ≡ 0 on the generalized qd-algorithm (8) and (9), where m is a positive integer, then we obtain a Lax pair for this generalized qd-equation
⎧ (n) (n) ⎪ e 0 ≡ em ≡ 0 ⎪ ⎪ ⎪ ⎨ 1 (n) 1 1 1 (n) (n+1) (n+1) qk − = e k −1 + qk − , ek + y y y y ⎪ n n +1 n+k n+k+1 ⎪ ⎪ ⎪ ⎩ q(n) e (n) = e (n−1) q(n−1) , k = 1, . . . , m. k k k k +1
k = 1, . . . , m − 1,
1226
C. Brezinski et al. / Applied Numerical Mathematics 60 (2010) 1221–1230
In fact, setting
⎛
(n)
q1 y n +1
⎜ ⎜ ⎜ ⎜ ⎜ ⎜ Ln = ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝
⎞
(n) (n)
q1 e 1 y n +1
1 yn
−
(n)
q2 y n +2
1
(n)
+ e1 −
1 yn
.. .. . . .. .
..
. ..
..
. (n)
(n)
⎜ ⎜ ⎜ Rn = ⎜ ⎜ ⎜ ⎝
1
en1
..
.
0
⎞
(n)
q m −1 e m −1 yn+m−1
.
1
⎛
⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
qm yn+m
(n) + em −1 −
1 yn
,
m×m
⎟ ⎟ ⎟ ⎟ .. (n) ⎟ . e m −1 ⎟ ⎠
..
.
1
m×m
then, we can write down its Lax form as follows Theorem 4. The generalized qd-equation has the Lax representation
L n +1 R n = R n L n . Proof. The result can be easily verified by definition of L n , R n , and the qd-equation.
2
As shown in [6,7,10], the original qd-algorithm can be used to compute the eigenvalues of a given tridiagonal matrix. Similarly, from Theorem 4, we see that the generalized qd-algorithm can also be used to compute the eigenvalues of a matrix having the same form as L n . 4. Numerical experiments In this section, we apply the generalized G-transformation based on the generalized rs-algorithm to compute infinite integrals of the form
∞ S=
f (x) dx. a
a
We set S n = a n f (x) dx, where a = a0 < a1 < · · · < an < · · · satisfying an → ∞, n → ∞, and xn = f (an ). 0 ∞In our algorithm, { yn } is an auxiliary sequence which can influence the convergence rate. It is obvious that S − S n = an f (x) dx. From Theorem 2 and the definition of the integral, we guess that a proper choice of y n may be yn = O (1/an ). The numerous numerical examples performed confirm that this choice leads to a much higher precision. Then, in the following numerical experiments, this is the choice which is made. Example 1. Consider the sequence
n3 Sn =
2 x3
dx
1
+∞
which converges to S = 1 (2/x3 ) dx = 1. Set yn = 1/an = 1/(3n2 + 3n + 1). The corresponding numerical results are presented in the following tables.
C. Brezinski et al. / Applied Numerical Mathematics 60 (2010) 1221–1230
(n)
1227
(n)
|G 2 − 1|
|G 3 − 1|
n
yn = 1
yn = 1/an
yn = 1
yn = 1/an
1 2 3 4 5 6 7
0.000260 0.0000643 0.0000164 0.0000052 0.0000019 0.0000008 0.00000042
0.0005331 0.0000667 0.0000140 0.0000040 0.0000014 0.0000006 0.0000003
0.0000376 0.0000075 0.0000024 0.0000009 0.0000004
0.00000464 0.00000143 0.00000047 0.00000018 0.00000008
(n)
(n)
|G 4 − 1|
|G 5 − 1|
n
yn = 1
yn = 1/an
yn = 1
yn = 1/an
1 2 3
0.0000012 0.0000011 0.0000004
0.000000430 0.000000041 0.000000003
0.0000010
0.0000002
Example 2. Consider the sequence ( np )3
Sn = 0
2 1 + x2
dx
∞
which converges to S = 0 2/(1 + x2 ) dx = π . Take p = 10 and yn = 1/an = 1/(3(np )2 + 3np + 1). The corresponding numerical results are presented in the following tables.
(n)
(n)
|G 2 − π |
|G 3 − π |
n
yn = 1
yn = 1/an
yn = 1
yn = 1/an
1 2 3 4 5 6 7
0.00003497 0.00001403 0.00000673 0.00000370 0.00000223 0.00000145 0.00000099
0.00003388 0.00001156 0.00000522 0.00000278 0.00000165 0.00000106 0.00000072
0.0000089 0.0000043 0.0000024 0.0000014 0.0000009
0.0000017 0.0000010 0.0000006 0.0000004 0.0000002
(n)
(n)
|G 4 − π |
|G 5 − π |
n
yn = 1
yn = 1/an
yn = 1
yn = 1/an
1 2 3
0.00000281 0.00000160 0.00000098
0.000000175 0.000000002 0.000000018
0.0000011
0.0000006
5. Relations with generalizations of the ε -algorithm, and conclusions In this paper, we proposed a generalization of the G-transformation. The corresponding G-algorithm, based on the rsalgorithm, is given. It is shown that the generalized rs-algorithm related to the generalization of the G-transformation is equivalent to a generalized qd-algorithm. Numerical examples with the generalization of the G-transformation are given. By some computation, we find that the generalized rs-algorithm has a close relation to the first generalization of the ε -algorithm similar to the relation between the original rs-algorithm and the ε -algorithm [2]. For example, when v n = Run , then we have the relations (n+1) ε2k = yn+k+1
(n)
sk
(n+1)
sk
(n) ε2k
yn
(n)
+
r k +1
ε(n+1) , (n+1) 2k−2
rk
1228
C. Brezinski et al. / Applied Numerical Mathematics 60 (2010) 1221–1230
(n) ε2k +2
(n+1)
(n+1)
− ε2k
= − yn · yn+1 · . . . · yn+k+1
(n+1) (n) ε2k − ε2k = yn · yn+1 · . . . · yn+k+1
where the numbers (n) ε− 1 = 0,
r k +1 (n)
,
sk+1
(n)
r k +1
(n+1)
,
sk
(n) ε2k are computed by the first generalization of the ε -algorithm
ε0(n) = S n , n = 0, 1, . . . ,
εk(n+)1 = εk(n−+11) +
yn
εk(n+1) − εk(n)
k , n = 0, 1 , . . . .
,
Based on this observation, it is natural to inquire if there exits some algorithm similar to the rs-algorithm that has some relation to the second generalization of the ε -algorithm. From a generalization of the G-transformation, we found an algorithm similar to the rs-algorithm which gave a definite answer to this question. This generalization of the G-transformation is defined by
un+k v n+k .. . ¯ k −1 v n+k ¯G (n) = R k yn+k v n+k .. . ¯ k −1 R v n+k
.. .
R¯ k−1 v n+k+1 yn+k+1 v n+k+1
¯ k −1
R
.. . · · · R¯ k−1 v n+2k ··· yn+2k ··· v n+2k .. . k − 1 · · · R¯ v n+2k ··· ···
un+k+1 v n+k+1
.. .
v n+k+1
un+2k v n+2k
¯ ¯ where { yn }∞ −1 are auxiliary variables, un = S n yn , v n = xn yn , and the operator R is defined by R v n = ( v n−1 / yn−1 ) and R¯ k+1 v n = ( R¯ k v n−1 / yn−1 ). Set
v n+k v n+k−1 R¯ v R¯ v n+k n+k−1 .. .. . . ¯ k −1 R v n+k−1 R¯ k−1 v n+k (n) r¯k = yn+k yn+k−1 v v n+k n+k−1 .. .. . . ¯ k −2 R v n+k−1 R¯ k−2 v n+k yn+k yn+k−1 v v n+k n+k−1 .. .. . . ¯ k −1 R v n+k−1 R¯ k−1 v n+k (n) s¯k = v n+k+1 v n+k R¯ v R¯ v n+k+1 n+k .. .. . . ¯ k −1 R v n+k R¯ k−1 v n+k+1 (n)
Then, the quantities r¯k (n)
s¯ 0 = yn−1 , (n)
(n+1)
(n)
(n+1)
and s¯k
r¯1 = xn , (n)
yn+2k
−
··· ··· ··· ··· ··· ··· ···
n = 0, 1 , . . . ,
r¯k+1 1 − (n+1) , r¯k+1 1
··· ··· ···
¯R k−1 v n+2k−2 , yn+2k−2 v n+2k−2 .. . ¯R k−2 v n+2k−2 yn+2k−1 v n+2k−1 .. . ¯R k−1 v n+2k−1 . v n+2k−1 R¯ v n+2k−1 .. . R¯ k−1 v n+2k−1 v n+2k−2
R¯ v n+2k−2 .. .
defined above can be computed by the recursive algorithm
(n)
s¯k+1 = s¯k r¯k+1 = r¯k
(n)
··· ···
1
k , n = 0, 1 , . . . , (n)
s¯k
(10)
, yn+k−1 s¯ (n+1) k
k = 1 , 2 , . . . , n = 0, 1 , . . . .
(11)
C. Brezinski et al. / Applied Numerical Mathematics 60 (2010) 1221–1230
1229
(n) The quantities G¯ k+1 can be computed by the recursive algorithm
G¯ 0 = S n , (n)
n = 0, 1 , . . . ,
(n+1)
G¯ k+1 = G¯ k (n)
(n+2)
−
G¯ k
(n+2)
(n+1)
− G¯ k
(n+1)
r¯k+1 ,
(n+1)
r¯k+1 − r¯k+1
k = 1 , 2 , . . . , n = 0, 1 , . . .
After some computation, we find that the algorithm (10)–(11) has a close relation to the second generalization of the
¯ n , then we have the relations ε -algorithm. For example, when v n = Ru (n+1) (n+2) (n+1) rk+1 (n+2) sk ε2k (n) ε2k = yn+k (n+1) − (n+2) ε2k−2 , yn+2k+1
sk
rk
(n) (n+1) ε2k = − yn+k+1 · yn+k+2 · . . . · yn+2k+2 +2 − ε2k (n+1)
ε2k
(n+1)
r k +1
(n+1)
,
sk+1
(n+1)
(n)
− ε2k = yn+k · yn+k+1 · . . . · yn+2k+1
where the numbers (n) ε− 1 = 0,
r k +1
(n+1)
sk
(n) ε2k are computed by the second generalization of the ε -algorithm
ε0(n) = S n , n = 0, 1, . . . ,
εk(n+)1 = εk(n−+11) +
yn+k
εk(n+1) − εk(n)
,
k , n = 0, 1 , . . . .
The generalized G-algorithm proposed in the paper proved to be interesting, for example, for computing infinite integrals of the form
∞ S=
f (x) dx. a
Let (an ) be a strictly increasing sequence of numbers tending to infinity. Setting
an Sn =
f (x) dx,
xk = f (ak ),
a
we have
∞ S − Sn =
f (x) dx = an
a ∞ k+1 k=n ak
f (x) dx
∞
ak f (ak ) =
k=n
∞
ak xk .
(12)
k=n
On the other hand, the kernel of the generalized G-transformation is
S − Sn =
1 b 0 v n + b 1 R v n + · · · + b k −1 R k −1 v n yn
= xn b0 − b1 / yn + · · · + (−1)k−1 bk−1 / ykn−1 + c 1 xn+1 + · · · + ck−1 xn+k−1
(13)
where the coefficients c i are functions of yn . Assuming that an tends to infinity with n, and comparing the coefficient of xn of the first term in (12) and (13), we guess that yn should satisfy
b0 − b1 / yn + · · · + (−1)k−1 bk−1 / ykn−1 → ∞,
n → ∞.
In our numerical experiments, good results were obtained for the choice yn = 1/an , while, in the original G-algorithm, an is usually constant, and, thus, yn also. In particular, the original G-algorithm is recovered for yn = 1. Acknowledgements This work was partially supported by the National Natural Science Foundation of China (Grant no. 10771207), and the knowledge innovation program of LSEC and the Institute of Computational Math., AMSS, CAS. C. Brezinski thanks X.B. Hu and the Institute for inviting him for a stay during which part of this work was done. The authors thanks the reviewers for interesting comments which helped to improve the presentation.
1230
C. Brezinski et al. / Applied Numerical Mathematics 60 (2010) 1221–1230
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12]
C. Brezinski, A general extrapolation algorithm, Numer. Math. 35 (1980) 175–187. C. Brezinski, M. Redivo Zaglia, Extrapolation Methods. Theory and Practice, North–Holland, Amsterdam, 1991. J.P. Delahaye, B. Germain-Bonne, Résultats négatifs en accélération de la convergence, Numer. Math. 35 (1980) 443–457. H.L. Gray, T.A. Atchison, Nonlinear transformations related to the evaluation of improper integrals. I, SIAM J. Numer. Anal. 4 (1967) 363–371. H.L. Gray, T.A. Atchison, G.V. McWilliams, Higher order G-transformations, SIAM J. Numer. Anal. 8 (1971) 365–381. P. Henrici, Applied and Computational Complex Analysis, Vol. 1, John Wiley & Sons, New York, London, Sydney, Toronto, 1974. Y. Nakamura, in: Y. Nakamura (Ed.), Applied Integrable Systems, Shokabo, Tokyo, 2000, pp. 189–200 (in Japanese). W.C. Pye, T.A. Atchison, An algorithm for the computation of the higher order G-transformation, SIAM J. Numer. Anal. 10 (1973) 1–7. H. Rutishauser, Der Quotienten-Differenzen-Algorithmus, Z. Angew. Math. Phys. 5 (1954) 233–251. H. Rutishauser, Solution of eigenvalue problems with the LR-transformation, Nat. Bur. Standards Appl. Math. Ser. 49 (1958) 47–81. A. Salam, On a generalization of the ε -algorithm, J. Comput. Appl. Math. 34 (1991) 455–464. P. Wynn, On a device for computing the em ( S n ) transformation, MTAC 10 (1956) 91–96.
Contents lists available at ScienceDirect
Applied Numerical Mathematics www.elsevier.com/locate/apnum
A generalization of the G-transformation and the related algorithms Claude Brezinski a , Yi He b,c , Xing-Biao Hu b,∗ , Jian-Qing Sun b,c a
Laboratoire Paul Painlevé, UMR CNRS 8524, UFR de Mathématiques Pures et Appliquées, Université des Sciences et Technologies de Lille, 59655-Villeneuve d’Ascq cedex, France LSEC, Institute of Computational Mathematics and Scientific, Engineering Computing, AMSS, Chinese Academy of Sciences, P.O. Box 2719, Beijing 100190, PR China c Graduate School of the Chinese Academy of Sciences, Beijing, PR China b
a r t i c l e
i n f o
Article history: Received 30 November 2009 Received in revised form 4 March 2010 Accepted 10 March 2010 Available online 18 March 2010 Keywords: G-transformation G-algorithm rs-algorithm qd-algorithm
a b s t r a c t In this paper, a generalization of the G-transformation is proposed together with the corresponding auxiliary rs-algorithm for implementation. We show that the related generalization of the rs-algorithm is equivalent to a generalized qd-algorithm. Applications of the generalization of the G-transformation to the computation of infinite integrals are also given. © 2010 IMACS. Published by Elsevier B.V. All rights reserved.
1. Introduction When a sequence ( S n ) is slowly converging to its limit S, it can be transformed, without altering its terms, into a new sequence ( T n ) which, under some assumptions, converges faster to the same limit, that is such that
lim
n→∞
Tn − S Sn − S
= 0.
Sequence transformations are based on the idea of extrapolation. There exist many such sequence transformations since a general transformation able to accelerate the convergence of all sequences cannot exist as proved in [3]. Each sequence transformation can only accelerate the convergence of some classes of sequences. In many cases, the terms of the sequence ( T n ), which also often depend on a second index denoted by k, are expressed as ratios of determinants, thus needing the use of a recursive algorithm for the implementation of the transformation. On sequence transformations and the related algorithms, consult [2]. In [4], the authors proposed a nonlinear sequence transformation called the G-transformation, and applied it to the computation of improper integrals. In [5], a higher order G-transformation was introduced, which was later proved to be a particular case of the E-transformation [1], the most general sequence transformation known so far. These transformations are implemented by the G-algorithm and the E-algorithm, respectively. The former was proposed, together with the rs-algorithm as its auxiliary algorithm, for implementing this higher transformation [8], while the latter was presented with the help of another set of recursion relations. As far as the G-transformation is concerned, the auxiliary rs-algorithm is less costly than the auxiliary algorithm for the E-algorithm.
*
Corresponding author. E-mail addresses: [email protected] (C. Brezinski), [email protected] (Y. He), [email protected] (X.-B. Hu), [email protected] (J.-Q. Sun).
0168-9274/$30.00 © 2010 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2010.03.008
1222
C. Brezinski et al. / Applied Numerical Mathematics 60 (2010) 1221–1230
It has been shown that the rs-algorithm is equivalent to the qd-algorithm [9,2]. The ε -algorithm [12] is a well-known extrapolation method, and two generalizations of this algorithm were proposed and studied in [2,11]. They all are particular cases of the G-transformation. Similarly to the generalizations of the ε -algorithm, we will now study a generalization of the G-transformation. By using determinantal identities, we will obtain a generalized G-algorithm and the related rs-algorithm for implementing the transformation. The generalized rs-algorithm reduces to the original rs-algorithm for a particular choice of the auxiliary sequence involved. The generalized rs-algorithm is equivalent to a generalized qd-algorithm. Imposing boundary value conditions on the generalized qd-algorithm, we will obtain a generalized qd-equation with a Lax representation. There are many connections between extrapolation algorithms and discrete integrable systems which could be considered as discretization of partial differential equations. An important property of such systems is their integrability, that is the possibility of obtaining their solution into a close form such as, for example, ratios of determinants. This is, in particular, the case of extrapolation algorithms such as the ε -algorithm and its generalizations, and the G-transformation which have a determinantal representation. Thus, the results presented into this paper are motivated by this connection. This paper is organized as follows. In Section 2, we give a generalization of the G-transformation, and the related G-algorithm and rs-algorithm. In Section 3, we obtain a generalized qd-algorithm, and show that it is equivalent to the generalized rs-algorithm related to the generalization of the G-algorithm. In Section 4, applications of the generalization of the G-transformation to the computation of improper integrals are given. Conclusions and a discussion are presented in Section 5. 2. A generalization of the G -transformation In this section, we give a generalization of the G-transformation. Then, we derive the generalized G-algorithm to implement it on the basis of a generalized rs-algorithm. This generalized G-transformation has the form
un v n .. . k −1 R vn (n) Gk = yn v n .. . k −1 R
vn
u n +1 v n +1
.. .
R k −1 v n +1 y n +1 v n +1
.. .
R k −1 v n +1
.. . k −1 ··· R v n+k , ··· yn+k ··· v n+k .. . k −1 ··· R v n+k ··· ···
un+k v n+k
where (xn ) and ( yn ) are auxiliary sequences, un = S n yn , v n = xn yn , and the operator R is defined by R v n = ( v n / yn ) and R k+1 v n = ( R k v n / yn ). Obviously, the condition yn = 0 for all n has to be satisfied. When ∀n, yn = 1, the general Gtransformation reduces to the original one. It must be noticed that this transformation enters into the general framework of the E-algorithm [1] by dividing the first columns in the numerator and the denominator by yn , the second ones by yn+1 , and so on. Set
vn Rv n .. . k −1 R vn (n) rk = yn v n .. . k −2 R
(n)
sk =
vn
yn vn R
vn Rv n .. . k −1 R
.. .
R k −1 v n +1 y n +1 v n +1
.. .
R k −2 v n +1 y n +1 v n +1
.. .
k −1
v n +1 R v n +1
.. .
vn
vn
k −1
R v n +1 v n +1 R v n +1
.. .
R k −1 v n +1
··· ···
v n+k−1 R v n+k−1
k −1 ··· R v n+k−1 , ··· yn+k−1 ··· v n+k−1 .. . · · · R k−2 v n+k−1 ··· yn+k ··· v n+k .. . k −1 ··· R v n+k . ··· v n+k−1 ··· R v n+k−1 .. . k −1 ··· R v n+k−1
Then, a generalized rs-algorithm holds
.. .
C. Brezinski et al. / Applied Numerical Mathematics 60 (2010) 1221–1230
1223
Theorem 1. (n) (n) The quantities rk and sk defined above can be recursively computed by the algorithm (n)
(n)
s 0 = 1,
r1 = xn ,
(n)
(n+1)
(n) r k +1
(n+1)
sk+1 = sk
= rk
(n+1) r k +1 (n) r k +1
n = 0, 1, . . . ,
−1 ,
k = 0, 1 , . . . , n = 0, 1 , . . . ,
(n+1)
sk
1
−
yn+k+1 s(n) k
1 yn
,
k = 1 , 2 , . . . , n = 0, 1 , . . . .
(n)
Proof. From the expressions of the quantities rk
(n)
sk+1
(n+1)
sk
y n vn . .. k R vn = vn R vn . .. k
··· ···
y n +1 v n +1
R v n +1
.. . R k v n +1
R vn
(n)
v n+k+1
.. .
R k v n +1 v n +1
(2)
and sk , we have
yn+k+1
.. .
(1)
· · · R k v n+k+1 · · · v n+k · · · R v n+k .. . k · · · R v n+k
v n +1 R v n +1 .. . k −1
··· ···
v n +2 R v n +2
v n+k R v n+k
.. .
.. .
R v n +1 y n +1
R k −1 v n +2 · · · y n +2 ···
v n +1
v n +2
.. .
.. .
R k −1 v n +1
R k −1 v n +2
R k−1 v n+k yn+k+1
···
v n+k+1
.. . · · · R k−1 v n+k+1
Applying Sylvester’s determinantal identity, we get
y n +1 ··· yn+k+1 v n+1 v n +2 ··· yn v v n +1 · · · v n+k+1 R v n+1 R v n +2 ··· n . .. .. .. .. . . . . . . k R v n R k v n+1 · · · R k v n+k+1 R k−1 v n+1 R k−1 v n+2 · · · v n +2 · · · v n+k+1 yn y n +1 v n +1 Rv R v n+2 · · · R v n+k+1 v n v n +1 n +1 = .. .. .. .. .. . . . . . k R v n+1 R k v n+2 · · · R k v n+k+1 R k−1 v n R k−1 v n+1 v n +1 · · · v n+k yn+1 y n +2 vn Rv R v n+1 · · · R v n+k v n+1 v n +2 n − . .. .. .. .. .. . . . . k R v n R k v n+1 · · · R k v n+k R k−1 v n+1 R k−1 v n+2
v n+k R v n+k
.. .
R k−1 v n+k ··· yn+k ··· v n+k
.. .
· · · R k−1 v n+k ··· yn+k+1 ··· v n+k+1 . .. . k −1 ··· R v n+k+1
Then, the relation (1) is proved. Let us now turn to the proof of the second relation. We have
(n) r k +1 (n+1)
rk
=
vn Rv n . . . k
R vn yn vn
.. .
R k −1 v n
v n +1 R v n +1
··· ···
R k v n +1 y n +1 v n +1
· · · R k v n+k R k−2 v n+1 R k−2 v n+2 · · · ··· yn+k v n+1 v n +2 ··· v n+k R v n+1 R v n +2 .. .. .. . . . · · · R k−1 v n+k R k−1 v n+1 R k−1 v n+2
.. .
.. .
R k −1 v n +1
v n+k R v n+k
.. .
y n +1 v n +1
.. .
y n +2 v n +2
.. .
··· ···
yn+k v n+k
.. .
R k−2 v n+k ··· v n+k ··· R v n+k
.. .
· · · R k−1 v n+k
.
By definition of the operator R, it holds
y n +1 v n +1
y n +2 v n +2
R k −2 v n +1
R k −2 v n +2
.. .
.. .
··· ···
yn+k v n+k
.. .
· · · R k−2 v n+k
··· R v n+k−1 R v n +1 .. .. = yn+1 yn+2 · · · yn+k . . . k −1 k −1 R v n +1 · · · R v n+k−1
.
1224
C. Brezinski et al. / Applied Numerical Mathematics 60 (2010) 1221–1230
Then, Sylvester’s identity gives
vn Rv n . . . k
··· ···
v n +1 R v n +1
.. .
v n+k R v n+k
.. .
y n +1 v n +1
.. .
y n +2 v n +2
··· ···
.. .
yn+k v n+k
.. .
· · · R k v n+k R k−2 v n+1 R k−2 v n+2 · · · R k−2 v n+k v n +1 · · · v n+k vn ··· R v n+k−1 Rv R v n +1 R v n+1 · · · R v n+k n .. .. = yn+1 yn+2 · · · yn+k . . . . . .. .. .. k −1 k −1 k R v n +1 · · · R v n+k−1 k k R v n R v n+1 · · · R v n+k ⎛ v n +1 ··· v n+k−1 vn R v n+1 · · · R v n+k R v R v n +1 ··· R v n+k−1 ⎜ n .. .. = yn+1 yn+2 · · · yn+k ⎜ .. .. .. . . ⎝ k . . . k R v n+1 · · · R v n+k k−1 R v n R k−1 v n+1 · · · R k−1 v n+k−1 ⎞ v n +2 ··· v n+k v n +1 R v n · · · R v n+k−1 R v R v n +2 · · · R v n+k ⎟ . n + 1 . ⎟ − .. .. .. .. ⎠ .. k . . . k R v n · · · R v n+k−1 k−1 R v n+1 R k−1 v n+2 · · · R k−1 v n+k y n +2 ··· yn+k+1 v n v n +1 ··· v n+k−1 y n +1 v v n +2 ··· v n+k+1 R v n R v n +1 ··· R v n+k−1 n +1 1 = .. .. .. .. .. .. yn+k+1 . . . . . . k −1 k −1 k −1 k −1 k −1 k −1 R v n +1 R v n +2 · · · R v n+k+1 R vn R v n +1 · · · R v n+k−1 y n +1 ··· yn+k v n+1 v n +2 ··· v n+k yn v v n +1 ··· v n+k R v n+1 R v n +2 ··· R v n+k n 1 . − . . . . . . .. .. .. .. .. .. yn k −1 k −1 k −1 k −1 k −1 k −1 R vn R v n +1 · · · R v n+k R v n +1 R v n +2 · · · R v n+k
R vn
R k v n +1
From the above relations, we obtain (2), which ends the proof of the theorem.
2
(n)
The quantities G k+1 can be recursively computed by (n)
G 0 = Sn, (n)
(n)
G k = G k −1 −
(3) (n+1) (n) G k−1 − G k−1 (n) r . (n+1) (n) k
rk
− rk
(4)
The kernel of the transformation is given by the following theorem Theorem 2. (n) A necessary and sufficient condition that ∀n, G k = S is that, ∀n,
( S n − S ) y n = b 0 v n + b 1 R v n + · · · + b k −1 R k −1 v n . 3. A generalized qd-algorithm In this section, we present a generalized qd-algorithm related to the generalized rs-algorithm (1)–(2). Set (n)
qk = (n)
ek =
(n+1) (n+1) sk−1 , (n) (n) rk sk−1 (n) (n) rk+1 sk . (n+1) (n+1) rk sk−1
rk
Then, the generalized qd-algorithm is as follows
(5)
(6)
C. Brezinski et al. / Applied Numerical Mathematics 60 (2010) 1221–1230
1225
Theorem 3. (n) (n) The qk and ek defined above can be computed by the generalized qd-algorithm (n)
(n)
e 0 = 0, 1
(n)
ek +
q1 = xn+1 /xn ,
yn+k
(n) (n)
1
(n)
qk −
yn
(7)
(n+1) = e k −1
1
+
(n+1)
yn+k+1
1
−
qk
y n +1
(8)
,
(n−1) (n−1) qk+1 .
qk e k = e k
(9)
Proof. Substituting (5) and (6) into both sides of (9) immediately shows that this formula holds. (n) (n) In order to prove formula (8), we first rewrite qk and ek using (1)–(2) as (n)
qk = (n)
ek =
(n+1) (n+1) sk−1 (n) (n) rk sk−1 (n) rk+1 sk(n) (n+1) (n+1) rk sk−1
rk
(n+1)
=
rk
r k +1
(n+1)
r k −1
yn
(n+1)
(n+1)
r yn+k rk yn+k = yn+k k(n+1) + , yn r (n) r k −1 k
(n) (n) r r 1 − 1 = k(+n)1 − (nk+ +1 )
(n+1)
rk
(n)
rk
1
+
(n+1)
(n)
rk
(n)
=
(n)
rk
rk
rk
rk
and (n)
qk = (n)
ek =
(n+1) (n+1) sk−1 (n) (n) rk sk−1 (n) rk+1 sk(n) (n+1) (n+1) rk sk−1
rk
=
(n)
sk
(n+1)
+1
1
sk
sk−1
=
(n+1)
sk−1 (n)
sk−1
(n+1)
1
−
yn+k+1 s(n) k
(n+1)
(n)
sk
=
(n)
sk−1
+
,
(n)
sk−1
(n)
sk
(n+1)
yn
sk−1
sk−1
=
(n+1)
sk
1
yn+k+1 s(n+1) k −1
−
1
(n)
sk
. yn s(n+1) k −1
Then, for the left-hand side of (8), we get
LHS of (8) =
1
yn+k
(n+1)
rk
yn+k (n+1) + · (n) yn r k −1 rk
(n+1)
=
rk
(n+1)
r k −1
+
1
rk
(n+1) r k −1
+
(n+1)
sk
yn+k+1 s(n+1) k −1
(n+1)
=
(n+1)
rk
yn+k
1
(n+1)
sk
1
+
+
rk
(n+1)
sk
1
1
(n)
sk
− yn+k+1 s(n+1) yn s(n+1) k −1 k −1
(n+1)
rk
(n)
yn
yn+k+1 s(n+1) k −1
(n)
−1−
sk
−
1 yn
(n+1)
sk−1
by (1) .
For the right-hand side of (8), we have (n+1)
(n+1)
r r 1 RHS of (8) = k(n+1) − k(n+2) + y n + k +1 r k −1 r k −1 (n+1)
=
rk
(n+1)
r k −1
+
rk
(n+1) r k −1
+
(n+1)
1
sk
1
sk
+ yn+k+1 s(n+1) k −1
(n+1)
=
(n+1)
yn+k+1 s(n+1) k −1
(n+2)
(n+1)
sk
(n+1) sk−1
+
sk−1
−
1
(n+1) y n +1 sk−1 (n+2) (n+1) sk−1 rk 1 1 − − ( n + 1 ) ( n + 2 ) yn+k+1 s y n +1 r k −1 k −1
by (2)
= LHS of (8), which completes our proof.
2 (n)
(n)
If we impose the boundary condition e 0 ≡ em ≡ 0 on the generalized qd-algorithm (8) and (9), where m is a positive integer, then we obtain a Lax pair for this generalized qd-equation
⎧ (n) (n) ⎪ e 0 ≡ em ≡ 0 ⎪ ⎪ ⎪ ⎨ 1 (n) 1 1 1 (n) (n+1) (n+1) qk − = e k −1 + qk − , ek + y y y y ⎪ n n +1 n+k n+k+1 ⎪ ⎪ ⎪ ⎩ q(n) e (n) = e (n−1) q(n−1) , k = 1, . . . , m. k k k k +1
k = 1, . . . , m − 1,
1226
C. Brezinski et al. / Applied Numerical Mathematics 60 (2010) 1221–1230
In fact, setting
⎛
(n)
q1 y n +1
⎜ ⎜ ⎜ ⎜ ⎜ ⎜ Ln = ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝
⎞
(n) (n)
q1 e 1 y n +1
1 yn
−
(n)
q2 y n +2
1
(n)
+ e1 −
1 yn
.. .. . . .. .
..
. ..
..
. (n)
(n)
⎜ ⎜ ⎜ Rn = ⎜ ⎜ ⎜ ⎝
1
en1
..
.
0
⎞
(n)
q m −1 e m −1 yn+m−1
.
1
⎛
⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
qm yn+m
(n) + em −1 −
1 yn
,
m×m
⎟ ⎟ ⎟ ⎟ .. (n) ⎟ . e m −1 ⎟ ⎠
..
.
1
m×m
then, we can write down its Lax form as follows Theorem 4. The generalized qd-equation has the Lax representation
L n +1 R n = R n L n . Proof. The result can be easily verified by definition of L n , R n , and the qd-equation.
2
As shown in [6,7,10], the original qd-algorithm can be used to compute the eigenvalues of a given tridiagonal matrix. Similarly, from Theorem 4, we see that the generalized qd-algorithm can also be used to compute the eigenvalues of a matrix having the same form as L n . 4. Numerical experiments In this section, we apply the generalized G-transformation based on the generalized rs-algorithm to compute infinite integrals of the form
∞ S=
f (x) dx. a
a
We set S n = a n f (x) dx, where a = a0 < a1 < · · · < an < · · · satisfying an → ∞, n → ∞, and xn = f (an ). 0 ∞In our algorithm, { yn } is an auxiliary sequence which can influence the convergence rate. It is obvious that S − S n = an f (x) dx. From Theorem 2 and the definition of the integral, we guess that a proper choice of y n may be yn = O (1/an ). The numerous numerical examples performed confirm that this choice leads to a much higher precision. Then, in the following numerical experiments, this is the choice which is made. Example 1. Consider the sequence
n3 Sn =
2 x3
dx
1
+∞
which converges to S = 1 (2/x3 ) dx = 1. Set yn = 1/an = 1/(3n2 + 3n + 1). The corresponding numerical results are presented in the following tables.
C. Brezinski et al. / Applied Numerical Mathematics 60 (2010) 1221–1230
(n)
1227
(n)
|G 2 − 1|
|G 3 − 1|
n
yn = 1
yn = 1/an
yn = 1
yn = 1/an
1 2 3 4 5 6 7
0.000260 0.0000643 0.0000164 0.0000052 0.0000019 0.0000008 0.00000042
0.0005331 0.0000667 0.0000140 0.0000040 0.0000014 0.0000006 0.0000003
0.0000376 0.0000075 0.0000024 0.0000009 0.0000004
0.00000464 0.00000143 0.00000047 0.00000018 0.00000008
(n)
(n)
|G 4 − 1|
|G 5 − 1|
n
yn = 1
yn = 1/an
yn = 1
yn = 1/an
1 2 3
0.0000012 0.0000011 0.0000004
0.000000430 0.000000041 0.000000003
0.0000010
0.0000002
Example 2. Consider the sequence ( np )3
Sn = 0
2 1 + x2
dx
∞
which converges to S = 0 2/(1 + x2 ) dx = π . Take p = 10 and yn = 1/an = 1/(3(np )2 + 3np + 1). The corresponding numerical results are presented in the following tables.
(n)
(n)
|G 2 − π |
|G 3 − π |
n
yn = 1
yn = 1/an
yn = 1
yn = 1/an
1 2 3 4 5 6 7
0.00003497 0.00001403 0.00000673 0.00000370 0.00000223 0.00000145 0.00000099
0.00003388 0.00001156 0.00000522 0.00000278 0.00000165 0.00000106 0.00000072
0.0000089 0.0000043 0.0000024 0.0000014 0.0000009
0.0000017 0.0000010 0.0000006 0.0000004 0.0000002
(n)
(n)
|G 4 − π |
|G 5 − π |
n
yn = 1
yn = 1/an
yn = 1
yn = 1/an
1 2 3
0.00000281 0.00000160 0.00000098
0.000000175 0.000000002 0.000000018
0.0000011
0.0000006
5. Relations with generalizations of the ε -algorithm, and conclusions In this paper, we proposed a generalization of the G-transformation. The corresponding G-algorithm, based on the rsalgorithm, is given. It is shown that the generalized rs-algorithm related to the generalization of the G-transformation is equivalent to a generalized qd-algorithm. Numerical examples with the generalization of the G-transformation are given. By some computation, we find that the generalized rs-algorithm has a close relation to the first generalization of the ε -algorithm similar to the relation between the original rs-algorithm and the ε -algorithm [2]. For example, when v n = Run , then we have the relations (n+1) ε2k = yn+k+1
(n)
sk
(n+1)
sk
(n) ε2k
yn
(n)
+
r k +1
ε(n+1) , (n+1) 2k−2
rk
1228
C. Brezinski et al. / Applied Numerical Mathematics 60 (2010) 1221–1230
(n) ε2k +2
(n+1)
(n+1)
− ε2k
= − yn · yn+1 · . . . · yn+k+1
(n+1) (n) ε2k − ε2k = yn · yn+1 · . . . · yn+k+1
where the numbers (n) ε− 1 = 0,
r k +1 (n)
,
sk+1
(n)
r k +1
(n+1)
,
sk
(n) ε2k are computed by the first generalization of the ε -algorithm
ε0(n) = S n , n = 0, 1, . . . ,
εk(n+)1 = εk(n−+11) +
yn
εk(n+1) − εk(n)
k , n = 0, 1 , . . . .
,
Based on this observation, it is natural to inquire if there exits some algorithm similar to the rs-algorithm that has some relation to the second generalization of the ε -algorithm. From a generalization of the G-transformation, we found an algorithm similar to the rs-algorithm which gave a definite answer to this question. This generalization of the G-transformation is defined by
un+k v n+k .. . ¯ k −1 v n+k ¯G (n) = R k yn+k v n+k .. . ¯ k −1 R v n+k
.. .
R¯ k−1 v n+k+1 yn+k+1 v n+k+1
¯ k −1
R
.. . · · · R¯ k−1 v n+2k ··· yn+2k ··· v n+2k .. . k − 1 · · · R¯ v n+2k ··· ···
un+k+1 v n+k+1
.. .
v n+k+1
un+2k v n+2k
¯ ¯ where { yn }∞ −1 are auxiliary variables, un = S n yn , v n = xn yn , and the operator R is defined by R v n = ( v n−1 / yn−1 ) and R¯ k+1 v n = ( R¯ k v n−1 / yn−1 ). Set
v n+k v n+k−1 R¯ v R¯ v n+k n+k−1 .. .. . . ¯ k −1 R v n+k−1 R¯ k−1 v n+k (n) r¯k = yn+k yn+k−1 v v n+k n+k−1 .. .. . . ¯ k −2 R v n+k−1 R¯ k−2 v n+k yn+k yn+k−1 v v n+k n+k−1 .. .. . . ¯ k −1 R v n+k−1 R¯ k−1 v n+k (n) s¯k = v n+k+1 v n+k R¯ v R¯ v n+k+1 n+k .. .. . . ¯ k −1 R v n+k R¯ k−1 v n+k+1 (n)
Then, the quantities r¯k (n)
s¯ 0 = yn−1 , (n)
(n+1)
(n)
(n+1)
and s¯k
r¯1 = xn , (n)
yn+2k
−
··· ··· ··· ··· ··· ··· ···
n = 0, 1 , . . . ,
r¯k+1 1 − (n+1) , r¯k+1 1
··· ··· ···
¯R k−1 v n+2k−2 , yn+2k−2 v n+2k−2 .. . ¯R k−2 v n+2k−2 yn+2k−1 v n+2k−1 .. . ¯R k−1 v n+2k−1 . v n+2k−1 R¯ v n+2k−1 .. . R¯ k−1 v n+2k−1 v n+2k−2
R¯ v n+2k−2 .. .
defined above can be computed by the recursive algorithm
(n)
s¯k+1 = s¯k r¯k+1 = r¯k
(n)
··· ···
1
k , n = 0, 1 , . . . , (n)
s¯k
(10)
, yn+k−1 s¯ (n+1) k
k = 1 , 2 , . . . , n = 0, 1 , . . . .
(11)
C. Brezinski et al. / Applied Numerical Mathematics 60 (2010) 1221–1230
1229
(n) The quantities G¯ k+1 can be computed by the recursive algorithm
G¯ 0 = S n , (n)
n = 0, 1 , . . . ,
(n+1)
G¯ k+1 = G¯ k (n)
(n+2)
−
G¯ k
(n+2)
(n+1)
− G¯ k
(n+1)
r¯k+1 ,
(n+1)
r¯k+1 − r¯k+1
k = 1 , 2 , . . . , n = 0, 1 , . . .
After some computation, we find that the algorithm (10)–(11) has a close relation to the second generalization of the
¯ n , then we have the relations ε -algorithm. For example, when v n = Ru (n+1) (n+2) (n+1) rk+1 (n+2) sk ε2k (n) ε2k = yn+k (n+1) − (n+2) ε2k−2 , yn+2k+1
sk
rk
(n) (n+1) ε2k = − yn+k+1 · yn+k+2 · . . . · yn+2k+2 +2 − ε2k (n+1)
ε2k
(n+1)
r k +1
(n+1)
,
sk+1
(n+1)
(n)
− ε2k = yn+k · yn+k+1 · . . . · yn+2k+1
where the numbers (n) ε− 1 = 0,
r k +1
(n+1)
sk
(n) ε2k are computed by the second generalization of the ε -algorithm
ε0(n) = S n , n = 0, 1, . . . ,
εk(n+)1 = εk(n−+11) +
yn+k
εk(n+1) − εk(n)
,
k , n = 0, 1 , . . . .
The generalized G-algorithm proposed in the paper proved to be interesting, for example, for computing infinite integrals of the form
∞ S=
f (x) dx. a
Let (an ) be a strictly increasing sequence of numbers tending to infinity. Setting
an Sn =
f (x) dx,
xk = f (ak ),
a
we have
∞ S − Sn =
f (x) dx = an
a ∞ k+1 k=n ak
f (x) dx
∞
ak f (ak ) =
k=n
∞
ak xk .
(12)
k=n
On the other hand, the kernel of the generalized G-transformation is
S − Sn =
1 b 0 v n + b 1 R v n + · · · + b k −1 R k −1 v n yn
= xn b0 − b1 / yn + · · · + (−1)k−1 bk−1 / ykn−1 + c 1 xn+1 + · · · + ck−1 xn+k−1
(13)
where the coefficients c i are functions of yn . Assuming that an tends to infinity with n, and comparing the coefficient of xn of the first term in (12) and (13), we guess that yn should satisfy
b0 − b1 / yn + · · · + (−1)k−1 bk−1 / ykn−1 → ∞,
n → ∞.
In our numerical experiments, good results were obtained for the choice yn = 1/an , while, in the original G-algorithm, an is usually constant, and, thus, yn also. In particular, the original G-algorithm is recovered for yn = 1. Acknowledgements This work was partially supported by the National Natural Science Foundation of China (Grant no. 10771207), and the knowledge innovation program of LSEC and the Institute of Computational Math., AMSS, CAS. C. Brezinski thanks X.B. Hu and the Institute for inviting him for a stay during which part of this work was done. The authors thanks the reviewers for interesting comments which helped to improve the presentation.
1230
C. Brezinski et al. / Applied Numerical Mathematics 60 (2010) 1221–1230
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12]
C. Brezinski, A general extrapolation algorithm, Numer. Math. 35 (1980) 175–187. C. Brezinski, M. Redivo Zaglia, Extrapolation Methods. Theory and Practice, North–Holland, Amsterdam, 1991. J.P. Delahaye, B. Germain-Bonne, Résultats négatifs en accélération de la convergence, Numer. Math. 35 (1980) 443–457. H.L. Gray, T.A. Atchison, Nonlinear transformations related to the evaluation of improper integrals. I, SIAM J. Numer. Anal. 4 (1967) 363–371. H.L. Gray, T.A. Atchison, G.V. McWilliams, Higher order G-transformations, SIAM J. Numer. Anal. 8 (1971) 365–381. P. Henrici, Applied and Computational Complex Analysis, Vol. 1, John Wiley & Sons, New York, London, Sydney, Toronto, 1974. Y. Nakamura, in: Y. Nakamura (Ed.), Applied Integrable Systems, Shokabo, Tokyo, 2000, pp. 189–200 (in Japanese). W.C. Pye, T.A. Atchison, An algorithm for the computation of the higher order G-transformation, SIAM J. Numer. Anal. 10 (1973) 1–7. H. Rutishauser, Der Quotienten-Differenzen-Algorithmus, Z. Angew. Math. Phys. 5 (1954) 233–251. H. Rutishauser, Solution of eigenvalue problems with the LR-transformation, Nat. Bur. Standards Appl. Math. Ser. 49 (1958) 47–81. A. Salam, On a generalization of the ε -algorithm, J. Comput. Appl. Math. 34 (1991) 455–464. P. Wynn, On a device for computing the em ( S n ) transformation, MTAC 10 (1956) 91–96.