Convergence acceleration via combined nonlinear-condensation transformations

J5& --f!B

Computer Physics Communications

ELSEMER

Physics Communications I16 ( 1999) 28-54

Computer

Convergence acceleration via combined nonlinear-condensation transformations Ulrich D. Jentschuraa~b~‘, Peter J. Mohra,‘, Gerhard Soff b,3, Ernst Joachim Weniger c,d*4 Insrirute

il National

h Institut ’ Center

Jijr

Theoretical

’ Irutituf

jiir

Studies Physikzlisrhe

of Standards

.ftir Theoretisrhe oj Physical und

and

Te&~lo~y,

Physik,

Gaithersburg,

TlJ Dresden.

Systems. Clark Atlanta Atlantu, GA 30314,

Theoretische

Chemie.

University USA

Universitiit

Received

30

MD

D-01062

20899-0001.

Dresden, James

Regensburg.

USA

Germany P. Bnrwiey

D-93040

Drive Regensburg,

at Fair

Street,

S. W,

Germany

July 1998

Abstract A method of numerically evaluating slowly convergent monotone series is described. First, we apply a condensation transformation due to Van Wijngaarden to the original series. This transforms the original monotone series into an alternating series. In the second step, the convergence of the transformed series is accelerated with the help of suitable nonlinear sequence transformations that are known lo be particularlypowerfulfor alternatingseries.Sometheoreticalaspectsof our approach are discussed. The efficiency, numerical stability, and wide applicability of the combined nonlinear-condensation transformation is illustrated by a number of examples. We discuss the evaluation of special functions close to or on the boundary of the circle of convergence, even in the vicinity of singularities. We also consider a series of products of spherical Bessel functions, which serves as a model for partial wave expansions occurring in quantum electrodynamic bound state calculations. @ 1999 Elsevier Science B.V. PAC.% 02.70.-c: 12.20.D~; Kqwords: Computational

atomic

and molecular

31.15.-p; 02.60:x techniques; Quantum electrodynamics physics; Numerical approximation and

(specific

1. Introduction Divergent

and slowly

convergent

series occur abun-

dantly in the mathematical and physical sciences.Accordingly, there is an extensive literature on numerical techniques which convert a divergent or slowly convergent seriesinto a new serieswith hopefully better ’ ’ ’ ’

E-mail: E-mail: E-mail: E-mail:

uljQnist.gov [email protected] [email protected] [email protected]

OOIO-4655/99/$ PII

SOOlO-465.5(98)001

- see front

matter I l-8

@ 1999 Elsevier

calculations);

Calculations

and mathematical

techniques

in

analysis

numerical properties. An overview of the existing sequencetransformationsaswell asmany referencescan be found in books by Wimp [ I] and Brezinski and Rcdivo Zaglia [ 21. The historical development of these techniquesup to 1945 is describedin a monograph by Brezinski [ 3 1, and the more recent developmentsare discussedin an article by Brezinski [4]. A very important classof sequences{sn}z+ is characterized by the asymptotic condition (1.1)

Science

B.V. All rights reserved

U.D. Jentschuru

et d/Computer

Physics

which closely resembles the well known ratio test for infinite series. Here, s = s, is the limit of this sequence as n ---f 00. A convergent sequence satisfying condition ( 1.1) with ]pI < 1 is called linearly convergent, and it is called logarithmically convergent if p= 1. If the elements of the sequence { sn}rb in Eq. ( 1.1) are the partial sums sn = C”,=, Uk of an infinite series, and if p satisfies either p = 1 or 0 < p < 1 (these are the only cases which will be considered in this article), then there exists an integer N > 0 such that all terms ak with k > N have the same sign. Hence, series of this kind can be split up into a finite sum containing the leading terms with the irregular signs, and a monotone series whose terms all have the same sign. The subject of this article is the efficient and reliable evaluation of series that exhibit these properties. The partial sums of a nonterminating Gaussian hypergeometric series, zFi(a,b;c;z)

=c

cG(a),(b),, zm m=O

(c)m

Z’

(1.2)

where (a), = r(a + m)/r(a) = a(a + 1) ..( (a + m - 1) is a Pochhammersymbol, or its generalization

,,+&(% =

c m

,?I*

. . . . ap+l;P1 ,...I (al), (Pl)nt

Pp;z)

. . . (ql+,)n2

zm

*. . (P,,,,

G-2

(1.3) ’

which converge for ]z 1 < 1 and diverge for Iz 1 > 1, are typical examplesof linearly convergent sequences with p = z. The partial sumsof the Dirichlet series

i(z&m+l)-’ nr=O

(1.4)

for the Riemann zeta function converge logarithmically if Re( z.) > 1. This follows from the following asymptotic estimate (n --+ 00) of the truncation error (p. 21 of Ref. [l]):

Communications

116 (1999)

28-54

29

The Dirichlet series (1.4) is notorious for its extremely slow convergence if Re( z ) is only slightly larger than one. In this case, the series can only be used for the computation of C(z) if it is combined with suitable convergence acceleration methods like the Euler-Maclaurin summation (see for instance Chapter 8 of Ref. [ 51, p. 379 of Ref. [ 61, or Chapter 6 of Ref. [7]). The acceleration of linear convergence is comparatively simple, both theoretically and practically, as long as p in Eq. (1.1) is not too close to one. With the help of Germain-Bonne’sformal theory of convergence acceleration [ 81 and its extension (Section 12 of Ref. [9]), it can be decided whether a sequence transformation is capable of accelerating linear convergence or not. Moreover, many sequencetransformations are known which are capable of accelerating linear convergence effectively. Examples are the A2 process,which is usually attributed to Aitken [ 101 although it is in fact much older (p. 90 of Ref. [ 3) ), Wynn’s epsilon algorithm [ 11 1, which produces Pad& approximants if the input data are the partial sums of a formal power series,or Levin’s sequencetransformation [ 121 and generalizations (Sections 7-9 of Ref. [ 9]), which require as input data not only the elementsof the sequenceto be transformed but also explicit estimatesfor the correspondingtruncation errors. The acceleration of logarithmic convergence is much more difficult than the acceleration of linear convergence. Delahaye and Germain-Bonne [ 131 showed that no sequence transformation can exist which is able to accelerate the convergence of all logarithmically convergent sequences.Consequently, in the caseof logarithmic convergence the successof a convergenceacceleration processcannot be guaranteed unlessadditional information is available. Also, an analogue of Germain-Bonne’s formal theory of linear convergence acceleration [ 81 cannot exist. In spite of these complications, many sequence transformations are known which work reasonably well at least for suitably restricted subsets of the class of logarithmically convergent sequences,Examples are Richardson extrapolation [ 141, Wynn’s rho algorithm [ 151 and its iteration (Section 6 of Ref. [ 91) as well as Osada’smodification of the rho algorithm [ 161, Brezinski’s theta algorithm [ 17] and its iteration (Section 10 of Ref. [ 9]), Levin’s u and

30

U.D. Jent.rchuru

et al. /Computer

Physics

I! transformations I 121 and related transformations (Sections 7-9 of Ref. [ 91)) and the modification of the Ail’ process by Bjerstad, Dahlquist, and Grosse [ 181. However, there is a considerable amount of theoretical and empirical evidence that sequence transformations are in general less effective in the case of logarithmic convergence than in the case of linear convergence. Numerical stability is a very important issue. A sequence transformation can only accelerate convergence if it succeeds in extracting some additional information about the index-dependence of the truncation errors from a finite set s,, s,-1, . . . , s,+k of input data. Normally, this is done by forming arithmetic expressions involving higher weighted differences. However, forming higher weighted differences is a potentially unstable process which can easily lead to a serious loss of significant digits or even to completely nonsensical results. If the input data are the partial sums of a strictly alternating series, the formation of higher weighted differences is normally a remarkably stable process. Hence, a serious loss of significant digits is not to be expected if the partial sums of a strictly alternating series are used as input data in a convergence acceleration or summation process. If, however, the input data are the partial sums of a monotone series, numerical instabilities due to cancellation are quite likely, in particular if convergence is very slow. Thus, if the sequence to be transformed either converges linearly with a value of p in Eq. ( 1.1) that is only slightly smaller than one, or if it converges logarithmically (p = I), numerical instabilities are a serious problem and at least some loss of significant digits is to be expected. Generally, the sequence transformations mentioned above are not able to determine the limit of a logarithmically convergent series, whose terms ultimately all have the same sign, with an accuracy close to machine accuracy. This restricts the practical usefulness of sequence transformations severely, e.g., in FORTRAN calculations with a fixed precision. In this paper, we show that these stability problems can often be overcome by transforming slowly convergent monotone series not by straightforward application of a single sequence transformation but by a combination of two di@rent transformations. In the.first step, a monotone series is transformed into a strictly alternating series with the help of a

Communicatiom

I16 (I 994) 28-54

condensation transformation. This transformation was first mentioned on p. 126 of Ref. [ 191 and only later published by Van Wijngaarden [20]. Later, the Van Wijngaarden transformation was studied by Daniel [ 211, and recently it was rederived by Pelzl and King [ 221, who used it for the high-precision calculation of atomic three-electron interaction integrals of explicitly correlated wave functions. In the second step, the convergence of the resulting alternating series is accelerated by suitable nonlineur sequence transformations [ 9,121 which are known to be very powerful in the case of alternating series. Since the transformation of alternating series is a remarkably stable process, the limits of even extremely slowly convergent monotone series can be determined with an accuracy close to machine accuracy. Conceptually, but not technically our approach resembles that of Brezinski, Delahaye, and Germain-Bonne [ 23 ] who proposed to extract a linearly convergent subsequence from a logarithmically convergent input sequence by a selection process. In this article, we will call our approach, which consists of the Van Wijngaarden condensation transformation and the subsequent nonlinear sequence transformation, the combined nonlinear-condensation trans,formation (CNCT) . The CNCT is not restricted to logarithmically convergent series. It can also be used in the case of a linearly convergent monotone series with a value of p in Eq. ( 1. I ) close to one. Typically, this corresponds to a power series whose coefficients ultimately all have the same sign and whose argument is positive and close to the boundary of the circle of convergence. We will present some examples which show that the CNCT works even if the argument of the power series is very close to a singularity. Our approach requires the evaluation of terms of the original monotone series with high indices. Consequently, our two-step approach is computationally more demanding than the application of a single sequence transformation, and it cannot be applied if only a few terms of a slowly convergent series are available. In spite of these restrictions, we believe that the CNCT is very useful at leas1 for certain problems since it is able to produce highly accurate results that can only be accomplished otherwise with a considerably greater numerical effort. Special functions are defined and, in many cases,

U.D.

Jenischura

et al. /Computer

Physics

also evaluated via series expansions. The evaluation of special functions is an old problem of numerical mathematics with a very extensive literature (compare for example the books by Luke [ 24,251, Van Der Laan and Temme [26], and Zhang and Jin [ 271). Nevertheless, there is still a considerable amount of research going on and many new algorithms for the computation of special functions have been developed recently (compare for example the papers by Lozier and Olver [28] and Lozier [ 291 and the long lists of references therein). We believe that sequence transformations in general and the CNCT in particular are useful tools for the evaluation of special functions. Of course, this is also true for problems in theoretical and computational physics. This paper is organized as follows: In Section 2, we describe the Van Wijngaarden transformation. In Section 3, we discuss the nonlinear sequence transformations which we use for the acceleration of the resulting alternating series. In Section 4, we apply the CNCT to the Riemann zeta function. In Section 5, we consider the evaluation of the Lerch transcendent and related functions with arguments on or close to the boundary of the circle of convergence. In Section 6, we examine the evaluation of the generalized hyper. . geometric series [>+I F, (p 1 2) with unit argument or with an argument z which is only slightly smaller than one. In Section 7, we discuss the evaluation of an infinite series involving Bessel and Hankel functions. In Appendix A, we discuss the efficiency of sequence transformations in the case of monotone series or strictly alternating series in more detail. Finally, in Appendix B we discuss exactness properties of the sequence transformations which we use in the second step for the acceleration of the convergence of the resulting alternating series. The example of Section 7 serves as a model problem for slowly convergent series, which occur in quantum electrodynamic bound state calculations and which were treated successfully with the methods discussed in this paper [30]. Therefore, we expect the CNCT to be a general computational tool for the evaluation of slowly convergent sums over intermediate angular momenta which arise from the decomposition of relativistic propagators in QED bound state calculations into partial waves.

Communications

I I6 (I 999)

28-54

31

All calculations were done in Mathematica5 with a relative accuracy of 16 decimal digits [ 311. In this way, we simulate the usual DOUBLE PRECISION accuracy in FORTRAN.

2. The Van Wijngaarden

transformation

Let us assume that the partial sums ,I rTn = lE

a(k)

(2.1)

k=O

of an infinite series converge either linearly or logarithmically to some limit c = go0 as n --+ 30. We also assume that all terms a(k) have the same sign, i.e.. the series CE? a(k) is a monotone series. Following Van Wijngaarden 1201, we transform the original series into an alternating series CFo( - 1),jA.i according to Fa(k)

= g(-l)‘Aj,

k=O

(2.2)

j=O cc

(2.3) b:” = 2ka(2k(j

+ 1) - 1).

(2.4)

Obviously, the terms A/ defined in Eq. (2.3) all have the same sign if the terms a(k) of the original series all have the same sign. In the sequel, the quantities Aj will be referred to as the condensed series, and the series CF,,( - 1)‘Aj will be referred to as the transformed alternating series, or alternatively as the Van Wijngaarden transformed series. We call the Van Wijngaarden transformation a condensation transformation because its close connection to Cauchy’s condensation theorem (p. 28 of Ref. [ 321 or p. 121 of Ref. [ 331). Given a monotone series I:? a(k) with terms that satisfy la( k+ 1) 1 < la(k) 1, Cauchy’s condensation theorem states that CEO a(k) s Certain commercialequipment,instruments, or materials are identified in this paper to foster understanding. Such identification does not imply recommendation or endorsement by the National Institute of Standards and Technology, nor does it imply that the materials or equipment identified are necessarily the best available for the purpose.

32

U. D. Jentsrhura

Table I Demonstration (-!)‘A, All

that Van Wijngaarden’s Leading a(0)

-Al

et al. /Computer

transformation,

terms of the original

Physics

Eqs. (2.2)-(2.4),

I16

is a rearrangement

(1999)

28-54

of the original

+40(3)

--Ai do)

a( k)

+

-2 a(3)

-

-n(3)

+ -

+a(2)

+a( I)

series c:,

series

+20(l) -a( 1)

A2

;;,iiW,

Communications

+a(?)

converges if and only if the first condensedseriesA0 defined in Eq. (2.3) converges. In Eq. (2.3)) the indices of the terms of the original seriesare chosen in such a way that sampling at very high indices takes place (according to Eq. (2.4), the indices of the terms of the original seriesgrow exponentially). In this way, we obtain information about the behaviour of the terms of the original seriesat high indices. Moreover, if the terms of the original seriesbehave asymptotically (n --+ co) either like a(n) N n-l-” with E > 0 or like a(n) N n@rnwith 0 < r < 1 and /3 real, then the terms of the original seriesbecome negligibly small after a few evaluations. Specifically, the series(2.3) for the terms A,j convergeslinearly in thesecases. When summingover k in Eq. (2.3), we found that it is normally sufficient to terminate this sumwhen the last term is smaller than the desired accuracy. Typically, 20 to 30 terms are neededfor a relative accuracy of lo-‘” in the final result for the condensedsum. It should be possibleto acceleratethe convergence of the series (2.3) for A/. Since, however, this is a monotoneseries,it is not clear whether and how many digits would be lost in the convergence acceleration process.Consequently,we prefer to perform a safeand straightforward evaluation of the condensedsumsA,i and add up the terms of the seriesuntil convergence is reached, although this is most likely not the most efficient approach. The transformation from a monotone series to a strictly alternating series according to Eqs. (2.2)(2.4) is essentially a reordering of the terms a(k) of the original series. This is seen as follows. We first define the partial sums

fd4) +n(4)

faC.7)

Sn = 2(-l j=o

+

jJAj

(2.5)

of the Van Wijngaarden transformed original series.It can be shown easily that S, with n > 0 reproducesthe partial sumg,,, Eq. (2.1), which containsthe first n+ 1 termsof the original series.To illustrate this procedure, we present Table 1. Formal proofs of the correctness of this rearrangementcan be found in Ref. [ 2 1] or in the appendix of Ref. ] 221. Daniel [ 211 was able to formulate some mild conditions which guaranteethat the limits u and S of the partial sums(2.1) and ( 2.5), respectively, simultaneously exist and are equal, lim (+, = fl= lim S, = S. ,I-00 n-+cx2

(2.6)

For example, in the Corollary on p. 92 of Ref. [ 2 I] it was shown that if a strictly decreasing sequence {Mk}& of positive bounds exists which satisfy la( k)l < Mk for all k 2 0 and if cEO it’fk < co holds, then the original monotone seriesand the Van Wijngaarden transformed series on the right-hand side of Eq. (2.2) both converge to the same limit (i.c., g = S holds). This useful criterion is fulfilled by all seriesconsideredin this paper.

3. Nonlinear sequencetransformations The series(2.3) for the terms of the Van Wijngaarden transformed seriescan be rewritten as follows: A, = a(j) + 2a(2j + 1) + 4a(4j + 3) +. ‘. . (3.1) Since the terms u(k) of the original series have by assumptionthe samesign, we immediately observe

U.D.

Jentschum

et d/Computer

Physics

(3.2) Consequently, the Van Wijngaarden transformation. Eqs. (2.2)-( 2.4)) does not lead to an alternating series whose terms decay more rapidly in magnitude than the terms of the original monotone series. Thus, an acceleration of convergence can only be achieved ii the partial sums (2.5) of the Van Wijngaarden transformed series are used as input data in a convergence acceleration process. Since the Van Wijngaarden transformed series on the right-hand side of Eq. (2.2) is alternating if the terms of the original series all have the same sign, it is recommended to choose a suitable convergence acceleration method which is particularly powerful in the case of alternating series. A judicious choice of the convergence accelerator is of utmost importance since it ultimately decides whether our approach is numerically useful or not. Daniel [ 211 used the Euler transformation (Eq. (3.6.27) on p. 16 of Ref. [34]),

O”Fakes. (-l)k (I Jkuk = c x=0 k=O

Communication.~

I16

(1999)

28-54

new sequence {s~}~~ with better numerical properties [ 1,2,9]. The basic assumption of a sequence transformation is that the elements of the sequence {.sn}Eo to be transformed can be partitioned into a limit s and a remainder r, according to s, = s + r,

= (-l&l)“’ ix=0

k u,. 0m

(3.4)

The Euler transformation, which was published in its original version in 1755 on p. 281 of Ref. [ 351, is a series transformation which is specially designed for alternating series. It is treated in many books on numerical mathematics. Nevertheless, the Euler transformation is not a particularly efficient accelerator for the Van Wijngaarden transformed series considered in this article. In Table 2, we show some explicit results obtained by the Euler transformation. We also applied the Euler transformation to all other Van Wijngaarden transformed series presented in this paper, and we consistently observed that it is clearly less powerful than the transformations we used. Much better results can be expected from the more modern nonlinear sequence transformations which transform a sequence { s,}zo, whose elements may be the partial sums of an infinite series, into a

(3.5)

for all n > 0. Only in the case of some more or less trivial model cases will a sequence transformation be able to determine the limit s of the sequence {s,,},~~ exactly after a jnite number of steps. Hence, the elements of the transformed sequence {s~}~~ can be partitioned into the same limit s and a transformed remainder I:, according to s: = s + r;

(3.6)

for all n 2 0. In general, the transformed remainders r; are nonzero for all finite values of n. However, convergence is obviously increased if the transformed remainders {rh}Ee vanish more rapidly than the original remainders { r,}zo,

Here, A is the (forward) difference operator defined by Af(n) = f(n + 1) - f(n), and A$

33

(3.7) The best known nonlinear sequence transformation is probably Wynn’s epsilon algorithm [ 111 which produces Padt approximants if the input data are the partial sums of a formal power series. The epsilon algorithm is also known to work well in the case of alternating series. Consequently, it is an obvious idea to use the epsilon algorithm for the acceleration of the alternating series which we obtain via the Van Wijngaarden transformation. However, much better results can be obtained by applying sequence transformations which use explicit remainder estimates as input data in addition to the partial sums of the series to be transformed. Consequently, in this article we use exclusively the following two sequence transformations by Levin [ 121 and Weniger [ 91, respectively,

+ ,Wk-’sn/w,} L$qp,sn, w,)= Ak{(n Ak{(n+ P>k-‘/w,)

U.D. Jentschura

(p+n+j>k-l

(fl+n+k)k-1

et al. /Computer

Physics

Sn+j

ii---n+j (3.9)

(P+n+j)k-1

(P+n+k)k-1

1

Wn+j

Here, {s,,}z? is the sequence to be transformed, and {w,~}~~ is a sequence of truncation error estimates. The shift parameter /I has to be positive in order to permit n = 0 in Eqs. (3.8) and (3.9). The most obvious choice is p = 1. This has been the choice in virtually all previous applications of these sequence transformations, and we will use p = 1 unless explicitly stated. However, we show in Appendix B that in some cases it may be very advantageous to choose other values for 0. The numerator and denominator sums in Eqs. (3.8) and (3.9) can also be computed with the help of the three-term recursions (Eq. (7.2-8) of Ref. [ 91) L:‘;,(p)

= L:“+“(/3)

_ (P + II + k)(B + n + kjk-’ (/3+n+kfljk and (Eq. (8.3-7)

Lh”hP)

(3.10)

116 (I 999) 28-54

truncation errors of the input sequence {sn}zO. In most cases of practical interest this is not possible, since only the numerical values of a finite number of sequence elements is available. Hence, in practice one has to determine the remainder estimates from these numerical values. On the basis of purely heuristic and asymptotic arguments, Levin [ 121 proposed some simple remainder estimates which - when used in Eq. (3.8) - lead to Levin’s u, t, and u transformations. These remainder estimates can also be used in Eq. (3.9) (Section 8.4 of Ref. [9]). However, in the case of a strictly alternating series, the best simple truncation error estimate is the first term neglected in the partial sum (p. 259 of Ref. [33]). This is also the best simple estimate for the truncation error of a strictly alternating nonterminating hypergeometric series ~Fo( cy,,& -2) with (Y,/?, z > 0, which diverges strongly for every z It’ 0 (Theorem 5.12-5 of Ref. [ 361) . Consequently, in the caseof convergent or divergent alternating seriesit is a natural idea to use the first term neglected in the partial sum as the remainder estimate,as proposedby Smith and Ford [ 371. In this article, we always accelerate the convergence of the partial sumsS, of the Van Wijngaarden transformed series(2.3) which is strictly alternating if the original seriesis monotone. Thus, we useexclusively the sequencetransformations (3.8) and (3.9) in combination with the remainderestimate proposed by Smith and Ford [ 37I, w, = AS, = (- l)“+lA,+l

of Ref. [91)

.

(3.12)

This yields the following sequencetransformations (Eqs. (7.3-9) and (8.4-4) of Ref. [9]):

sg, (p) = $+“(p) (P+n+k>(P+n+k-(p+n+2k)(P+n+2k-

Communications

1) sjtt)(p)e 1)

(3.11)

The initial values Lg’( /?) = $“’ (p) = s,/o, and Lr) (p) = Sh”’ (p) = 1/o” produce the numerator and denominator sums, respectively, of lr)(P, s,, w,) and Si”)(P, s,,, w,). The performance of the sequence transformations li”‘(p, s,, w,,) and Si”‘(p, s,, w,) depends crucially on the remainder estimates {~,,}z~. Under exceptionally favourable circumstances it may be possible to construct explicit approximations to the

cd;“’ (p, s,) = P)k (,& sm A%> = C:“‘(P,S,, S:“‘(P,S,)

(-lY+‘A,+,)

,

(3.13)

.

(3.14)

=S~“‘(PS,,AS,)

= Scn)(p k r”,S (-l)“+‘A,+,)

Unless explicitly stated, we use p = 1. In the applications described in the following, it will be clear from the context which monotone series was transformed according to Van Wijngaarden such as to produce the input data for the nonlinear sequence transformations

U. D. Jentschuru

et al. /Computer

Physics

which are the partial sums S, of the alternating series (2.3). Alternative remainder estimates for the sequence transformations (3.8) and (3.9) were discussed in Sections 7 and 8 of Ref. [9] or in Refs. [38,39]. From a purely theoretical point of view, the sequence transformations d:“’ and Sp’ as well as their parent transformations Cy) and $“’ have the disadvantage that no general convergence proof is known. Only for some special model problems could rigorous convergence proofs be obtained (Refs. [ 40-421 or Sections 13 and 14 of Ref. [ 91) . However, there is overwhelming empirical evidence that d:“’ and SF) work very well in the case of convergent or divergent alternating series for instance as they occur in special function theory [ 2,9,12,37,43-481 or in quantum mechanical perturbation theory [ 44,45,49-551. Pelzl and King [22] only used the Levin transformation df’ (p, S,) f or t h e acceleration of the convergence of Van Wijngaarden transformed alternating series. However, we shall show that the closely related transformation Sifl) (/3, S,) frequently gives better results. Application of the sequence transformations (3.13) and (3.14) to the partial sums S,, of the infinite series (2.3) produces doubly indexed sets of transforms, which depend on the starting index n of the input data and the transformation order k, {S,1,S,,+-~r...rSn+k,Sn+ki-I}

+Tk(“).

(3.15)

Here, Tj”’ stands for either dk”) (p, S,) or 8:) (0, S,). The transforms can be displayed in a rectangular scheme called the table of the transformation. In principle, there is an unlimited variety of different paths, on which we can move through the table in a convergence acceleration or summation process. In this article, we always proceed on such a path that for a given set of input data, we use the transforms with the highest possible transformation order as approximations to the limit of the series to be transformed. Thus, given the input data {Ss,SI,. . . ,S,,,,S,,,+t}, we always use the sequence {Td’), 7’,(O), . . . , Ti,@} of transforms as approximations to the limit (compare for instance Eqs. (7-5.4) and (7-5.5) of Ref. [9] ), where T(O) ,,I should provide the best approximation since it has the highest transformation order.

Communications

I16

(I 999)

2X-54

35

In Appendix A, we explain why the transformations Lr’ (p, s,, w,) and Sin) (p, s,, w,) are much more effective in the case of an alternating series than in the case of a monotone series. 4. Tbe Riemann zeta function The Riemann zeta function is discussed in most books on special functions. Particularly detailed treatments can be. found in monographs by Edwards [ 71, Titchmarsh [ 561, IviC [ 571, and Patterson [ 581. Many applications of the zeta function in theoretical physics are described in books by Elizalde, Odintsov, Romeo, Bytsenko, and Zerbini [59] and by Elizalde [ 601. In this section, we discuss how the CNCT can be used for the evaluation of the Riemann zeta function, starting from the logarithmically convergent Dirichlet series (1.4). We do not claim that our approach is necessarily more powerful than other, more specialized techniques for the evaluation of the Riemann zeta function. However, because of its simplicity, the Dirichlet series ( 1.4) is particularly well suited for an illustration of our approach. In the case of the Riemann zeta function, the terms by’ of the series (2.3) for Aj are given by

,,j/” =

(21-‘)k (4.1)

(j’

Thus, the series (2.3) is a power series in 2td1 which converges linearly. Moreover, it is essentially the geometric series in the variable 2’-” that can be expressed in closed form according to A, =

1 1-

1 (j + 1)~

21-7

(4.2)

Inserting this into the infinite series on the right-hand side of Eq. (2.2) yields the following transformed alternating series for the Riemann zeta function:

[(.7)x

l

l-21-2

E”‘. J=o

(j+l)z

(4.3)

An alternative proof of the validity of this series expansion can be found in Section 2.2 of Titchmarsh’s book [56]. The terms of the alternating series (4.3) decay as slowly in magnitude as the terms of the logarithmi-

U.D. Jenrschuru

36

et al. /Computer

Phvics

tally convergent Dirichlet series ( I .4). Nevertheless, the series (4.3) offers some substantial computational advantages because it is alternating. For example, it converges for Re( z ) > 0, whereas the Dirichlet series converges only for Re( z ) > 1. For Re( z ) 2 0 both the Dirichlet series ( I .4) and the alternating series (4.3) diverge. However, the alternating series can provide an analytic continuation of the Riemann zeta function for Re( z) < 0 if it is combined with a suitable summation process. For this purpose, we can use the same nonlinear transformations (3.13) and (3.14) as in the case of a convergent Van Wijngaarden transformed alternating series. This observation extends the applicability of’ the CNCT, which was originally designed for convergent monotone series only. Of course, this extension to divergent series is only possible if the monotone series (2.3) for A,j can be summed. Normally, the summation of a monotone divergent series is very difficult. However, in the case of the Riemann zeta function this is trivial since the series (2.3) is - as remarked above - essentially the geometric series in the variable 2’-’ which can be expressed in closed form according to Eq. (4.2), no matter whether /2’-” / < 1 or 12’-‘/ > 1 holds. Our first numerical example is the zeta function with argument z = 1.Ol. The Dirichlet series ( 1.4) converges for this argument, but its convergence is so slow that it is computationally useless. It follows from the truncation error estimate (1 S) that on the order of 10600terms would be needed to achieve the modest relative accuracy of 10m6 if the Dirichlet series were summed by adding up the terms. A much more efficient evaluation of [(z ) near Re( z ) = I is based on the Euler-Maclaurin formula (Eys. (2.01) and (2.02) on p. 285 of Ref. [5])

I16 (I 999) 28-54

Here, [[xl] is the integral part of x, i.e., the largestinteger m satisfying m 5 X. Bk ( X) is a Bernoulli polynomial defined by the generatingfunction (Eq. ( 1.06) on p. 281 of Ref. IS]) texp(xt) exp(t)

- I

=g &n(x) f m=o

,

It/ < 297.) (4.5)

(4.6)

B, = B,,(O)

is a Bernoulli number (p. 28 1 of Ref. [ 51). In the next step,we rewrite the Dirichlet series( 1.4) asfollows:

N-l l(z)=C

’ +g ’ ,,,* Cm+ I)2 n,=NCm+ 1)’ ’

(4.7)

and apply the Euler-Maclaurin formula (4.4) to the infinite serieson the right-hand side. With N = 100, we easily obtain the value of the zeta function accurate to more than 15 significant decimal digits. The value of J( I .Ol ) to I5 decimal placesis lo-“l(

I .Ol ) = 0.100577 943 338497,

(4.8)

We presentthis result (as well as all other subsequent numerical results) normalized to a number in the interval (0, 1) by a multiplicative power of 10. In Table 2, we apply the Euler transformation (3.3) and the nonlinear sequencetransformations (3.13) and (3.14) to the partial sumsof the alternating series (4.3) with ;: = 1.01. We list as a function of n the partial sums

Sn= , -;+

2k=* (k+ (-ljklj2

(4.9)

of the alternating series(4.3), the partial sums

$&d = r” fW dx+ 4[f(N) + f(M)] n=N

Comrnunimtions

N

+t: ,I, &

[f’2J-“(M)

- f@-‘)(N)]

+ Rq’ (4.4a)

M

1

___ /B& Rq = - (2q) ! N

- Uxll)f’2q)(~)

E,, = 1 -:I-;

i$g;:$

0

(4.10)

of the Euler transformed series (3.3), and the nonlinear sequence transformations d,(,‘)(I, So) and

sp(l,s,).

dx. (4.4b)

The results in Table 2 show that the Euler transformation is in the case of the alternating series (4.3)

U.D.

Table 2 Evaluation of the Riemann is given as 10-j [( 1.01)

Jentschuru

zeta function

et al. /Computer

of argument

Physics

Communications

z = I .Ol by accelerating

116

(1999)

the convergence

28-54

37

of the alternating

series (4.3);

II

S,,

FJn

4”’

0 I 2 3 4 5 6 7 8 9 IO II I2 13 14 IS

0.144770081711084 0.072 885 040 855 542 0.120614482322980 0.084 920 235 019 068 0.113411984373.518 0.089712 109307 161 0.109994997614328 0.092 27 1 153 050 434 0.108007 136313467 0.093 859 665 080 579 0.106 708 750 240 925 0.094 940 666 205 327 0.105794821569601 0.095 723 429487 784 0.105 116912361076 0.096316203847728

0.144770081 711084 0.090 606 301069 428 0.096 697 48 1 2.52 857 0.098 985 546 018 036 0.099 90 I 837 494 047 0.100283957662399 0.100447835715031 0.100519S72S724S4 0.1005Sl4706S3282 0.100565 83059981 I 0.100572360 153981 0.100575353801 276 0.100576735870616 0.100 s77 377 707 430 0.100 577 677 299 954 0.100577817763434

0.144770081711084 0.101 569133 143252 0.100456642533059 0. IO0 587 783 459 042 0.100577428 866203 0.100577954415585 0.100577944 116204 0.100 577 943 249 050 0.100577943 342049 0.100 577 943 338 553 0.100577943 338482 0.100577943338498 0.100 577 943 338 497 0.100577943 338497 O.lOOS77943 338497 0.100 577 943 338 497

0.144770081711084 0.101569 133 143252 0.1004S6642S33OS9 O.lOOS79332613649 0.100578083572921 0.100 577 949 566 834 0.100.577943567122 0.100577943346 IS0 0.100577 943 338734 0.100 577 943 338 so3 0.100 577 943 338 497 0. IO0 577 943 338 497 0.100577943338497 0.100577943 338497 0.100577943338497 0.100577943338497

exact

0. too 577 943 338 497

O.lOOS77943338497

O.lOOS77943338497

0.100577943338497

much less effective than the nonlinear sequence transformations dh”) (I, Sa) and 8:“) (1, SO). When we applied Wynn’s epsilon algorithm [ 111 to the alternating series (4.3), we also obtained clearly inferior results. This is in agreement with the observations of Pelzl and King [ 22 1. Thus, in the following we will only consider the nonlinear transformations (3.13) and (3.14) which use explicit remainder estimates. In Table 2, we present apparently redundant data since our transformation results in the last two columns do not change for n 2 12. However, we include these results here as well as below in order to demonstrate that our transformation remains stable even if we increase the transformation order. If the argument z is zero or a negative integer, the Riemann zeta function satisfies (p. 807 of Ref. [ 341) i(O)

= -;, ‘.

3( 1 - 2m) = -2

(4.1 laj

J( --2m) = 0, (

m= 1,2,...

.

(4.1 lb)

Here, Bl,,, is a Bernoulli number defined in Eq. (4.6). Our second numerical example is the zeta function with argument e = -1 which, because B2 = l/6, satisfies

(1 >so)

[(--l)=-l/12=-0.0833333....

sy

the result

(I I St,)

(4.12)

As remarked above, the alternating series (4.3) diverges for z = -1. However, our results in Table 3 show that this series can be summed effectively by the nonlinear sequence transformations (3.13) and (3.14). The results in Table 3 indicate that the sequence transformation S$“) (1, Se) is exact for the alternating series (4.3) with z = - 1. The exactness of the sequence transformations cy) (0, s,, w,,), Eq. (3.8)) and Sk”){/?, s,, w,), Eq. (3.9), for the infinite series (4.3) and for related series is discussed in Appendix B. There is considerable research going on in connection with the Riemann-Siegel conjecture ] 61-64 1. In this context, it is necessary to evaluate the Riemann zeta function on the so-called critical line z = l/2 + ti (-CC < t < 00). This can also be accomplished by applying the sequence transformations (3.13) and (3.14) to the partial sums of the alternating series (4.3). In Table 4, where we treat the case z = l/2+ 13.7i, the Levin transformation outperforms the Weniger transformation. For limitations of space, we do not present the partial sums of the transformed alternating series in Table 4.

38

U.D. Jentschuru

Table 3 Evaluation lOL(-1)

of the Riemann

I,

Sn

0 I 2 3 4 5 6 7 8 9 IO II 12 13 I4 IS

-3.333 3.333 -h.666 6.666 - tn.000 10.000 - 13.333 13.333 - 16.666 16.666 -20.000 2(1.000 -23.333 23.333 -26.666 26.666

exact

-0.833

zeta function

333 333 666 666 000 000 333 333 666 666 000 000 333 333 666 666

333 333 666 666 000 000 333 333 666 666 000 000 333 333 666 666

333 333 666 666 000 000 333 333 666 666 000 000 333 333 666 666

et ~1. /Computer

of argument

; = -I

333 333 667 667 000 000 333 333 667 667 000 000 333 330 670 670

333 333 333 333

5. The Lerch transcendent

Ph,ysics Communicutions

by summing

the result is given

-3.333 333 333 333 -0.666 666 666 666 -0.860215053763441 -0.833 333 333 333 -0.833 333 333 333 -0.833 333 333 333 -0.833 333 333 333 -0.833 333 333 333 -0.833 333 333 333 -0.833 333 333 333 -0.833 333 333 333 -0.833 333 333 333 -0.833 333 333 333 -0.833 333 333 333 -0.833 333 333 333 -0.833 333 333 333

-0.833

-0.833

(5.2)

or the Jonquiere function (p. 33 of Ref. [ 65 ] )

flil =Z@(i,S, F(Z,S)$ ____ z I). &)(n+I)”

series (4.3);

-3.333 333 333 333 333 -0.666 666 666 666 667 -0.860215 053763441 -0.830449826989619 -0.833557890954819 -0.833319627418409 -0.833334020666741 -0.833 333 304 093 535 -0.833333334413 139 -0.833 333 333 298 109 -0.833 333 333 334 362 -0.833 333 333 333 306 -0.833 333 333 333 334 -0.833 333 333 333 333 -0.833 333 333 333 333 -0.833 333 333 333 333

and related functions

1)

alternating

fii”’ ( I * so)

contains many special functions as special cases, for example the Riemann zeta function =@(l,s,

28-54

d:“’ ( I1 so)

The Lerch transcendent (p. 32 of Ref. [ 651)

l(s)

the divergent

I I6 (1999)

(5.3)

333 333 333 333

(5.4)

333 667 333 333 333 333 333 333 333 333 333 333 333 333 333

333 333 333 333

physics. They play a role in Bose-Einstein condensation [ 67,681, and they are particularly important in quantum field theory. For example, they occur in integrands of one-dimensional numerical integrations in quantum electrodynamic bound state calculations [ 69-731. The series expansion (5.1) for @( Z, S, LY) converges very slowly if the argument z is only slightly smaller than one. However, the CNCT can be used for an efficient and reliable evaluation of the Lerch transcendent and its special cases. In our first example (Table 5 ) , we evaluate 10W2Lii(0.99999)

= -IO-*1n(0.00001)

=0.115 129254649702.

In the physics literature, the Jonquiere function is usually called generalized logarithm or polylogarithm, and the following notation is used:

as

The series (5.4) lz/ < 1, but =

(5.5)

for Li,s( z) converges linearly for

zk+I

Li,(z)=C-

=-ln(1

-z)

(5.6)

f(=Ok + ’ This is also the notation and terminology used in the book by Lewin [ 661. Lerch transcendents and their special cases, the polylogarithms. are ubiquitous in theoretical

has a singularity at z = 1, Thus, in Table 5 we evaluate Lii (z) in the immediate vicinity of a singularity. In contrast, the series expansions

U.D.

Table 4 Evaluation

of [(l/2

+ 13.7i)

Jentschum

et al./Computer

Physics

Communications

116

(1999)

28-54

39

with the CNCT

sy (1,so)

7 8 9 IO II 12 13 14 IS I6 17 18 19 20 21 12 23 25 2s

0.414 107543949 134 0.575 871239097 112 0.523424912 174020 0.481953715 196 159 0.012442899246184 0.123074021609358 0.105 377 569 236 I75 0.107635288 132001 0.107 429 326 957 578 0.107 438 933 679469 0.107439640888613 0.107439434 156 190 0.107439457 152512 0.107439455840 151 0.107 439 455 825 365 0.107439455836355 0. IO7 439 455 835 269 0.107439455835311 0.1074394SS83S313 0.107439455835313 0.1074394SS835313 0.10743945.583S313 0.107439455835313 0.1074394SS83S313 0.107439455835313 0.107439455835313

exact

0.10743945S83S313

0 2 3 4 5 6

+ +

-

0.017316297 1257901 0.042690435 565 7.58 i 0. IS2 835 043959 961 i 0.288 400 086 031 923 i 0.237 603 260 694 125 i 0.316357718264423i 0.313246538788829; 0.312843 180304712i 0.312999708812577; 0.312974 1843086751 0.312976813188762i 0.312976661398796i 0.312976659181027i 0.3129766607172SSi 0.312976660547552i

-

0.312976660556014i 0.312976660556233i 0.312976660556 0.312976660556 0.312976660556 0.312976660556 0.3 I2 976 660 5.56 0.3 I2 976 660 556 0.312976660556 0.3 I2 976 660 556 0.3 I2 976 660 5.56

163i I63 i 163 i 163i 163 i 163 i

-

0.312976660556

163;

-

-

157i 163i 163i

(5.7)

0.414 107 543 949 134 0.57.5871239097112 0.567 958 887 269 553 0.474129917411618 -0.180 827 868 994 142 0. I26 392 529 409 594 0.107 386 234 787 298 0.107 124 24 1668 490 0.107489 121 873498 0.107436 183812222 0.107439 393557 558 0.107439488450447 0.107439 453 065 793 0.107 439 455 877 677 0.107 439 455 848 368 0.107439455833989 0. IO7 439 455 835 348 0.107439455835 317 0.107439455835312 0.10743945.5835313 0.107439455835313 0.107439455 835 313 0.107439455 835 313 0.107439455835313 0.107439455835313 0.107439455835313 0.107439455835313

+

-

0.017316297 12579Oi 0.042 690435 565 758 i 0.12991338648622Oi 0.168775423138291i 0.367 542 940 737 OS 1 i 0.290 235 127 404 228 i 0.317087400005856i 0.312622662817549i 0312978336617038i 0.312980 180873835 i 0.3 12 976 229 866 877 i 0.3 12 976 678 675 422 i 0.312976661 841 904i 0.312976660319609i 0.312976660568977i 0.31297666OS5644Oi 0.312976660556 072i 0.312976660556 169i 0.312976660556 163i 0.312976660556 163i 0.312976660556 163i O.312976660SS6163i 0.312976660556 l63i 0.312976660556163i r).312976660556 163i 0.312976660556 163i

-

0.312976660556

+ -

163i

and 7 somewhat better results than the Levin transform dA”’ (1, SO). In our last example of this section, we evaluate

and

104@(0.99999,2,10000)

Zk+l

Li,(z)

= fi: k=O (k+

1)”

converge logarithmically we consider 10-l L&(0.99999)

(5.8) for z = 1. In Tables 6 and 7,

= 0.164480893699293..

. (5.9)

and lo-‘Li3(0.99999)

= 0.798 585 139 222 548

(5.11)

with the help of the CNCT. The prefactor of 10000 is introduced so that the final result is of order one.

6. The generalized The

Gaussian

hypergeometric

series

hypergeometric

function

2Fl (a, b; c; z ), defined by the series expansion ( 1.2)

= 0.120204045438733, (5.10)

respectively. As in most previous examples, the Weniger transform 8:‘) (1, SO) produces in Tables 6

for Iz / < 1, is one of the most used and best understood special functions. It is discussed in virtually all books on special functions. A particularly detailed treatment can be found in the book by Slater [ 741.

Table 5 Evaluation

of IO-*Li1(0.99999)

= -IO-'

ln(0.00001)

with

the CNCT

Ei .n

4”

I

0.162768973713089 0.086384436856544

0.162768973713089 0.116225388785336

0.162768973713089 0.116225388785336

2 3 4 5

0.135357615351010 0.099665296922988 0.12757.5317431702 0.1047S53448S1635

0.114995657006664 0.115 140148 148939 0.115 128665188679 0.115 129271216942

0.114995657006664 0.115131002772470 0.I15129400970919 0.115129261517373

6 I 8

0.123997631730577 0.107401422517315 0.121964815378398

0.115 0.115 0.115

129254851 129254746908 129254612355

O.llS1292S4830226 0.115 129254725403 0.115129254664668

9 10 II

0.109009755 0.120662087466136 0.110085384510686

0.115 0.115 0.115

129254647759 129254650661 129254650507

0.115 129254647723 0.115129254648068 0.1151292S4649602

12 IS I4

0.119759673589890 0.110852765866205 0.119099529780039

0.115 129254649358 O.llS129254649714 0.115129254649727

0.115 129254649772 0.11512925464971I 0.115 129254649700

IS I6 17 IX

0.111426375175159 0.118596725346641 0.111 870534473655 0.118201666639655

0.115129254649696 O.llSl292S4649702 O.llSl29254649703 0.115 129254649702

0.115 129254649702 0.11Sl29254649702 0.115 129254Gl9702 0.115 129254649702

19 20

0.112224086SlS227 0.117883SOS7S271S

0.115129254649702 O.llS 129254649702

0.115129254649702 0.11512Y254649702

exact

0.115129254649702

O.ll5

0.115

IZ

0

Table 6 Evaluation

of

IO-'Liz(O.99999)

12.5042

with

( 1 I SC,)

129254649702

082

( I1 so)

129254649702

the CNCT

0.199982280324442 0.149990640162221 0.172207484145972 0.1597114140661I1 0.167708334Sl4884

0.199982280324442

0.162155301413564 0.166234803732817 0.163111643694762

0.1644808971283S3 0.l64480894717325 0.164480893606728

0.165S791628SSS83 0.163580602633 0.165232198806363

0.164480893702656 0.164480893699380 0.164480893699272

0.163844487813342 0.165026841360043 0.164007426937253 0.164895394970447

s?’

0.165371886328955 0.164381453915497 0.164488426505073 0.16448OS38599000

I96

0.164480893699295 0.164480893699292 0.164480893699293

0.199982280324442 0.165371886328955 0.164381453915497 0.164482760527739 0.164481025806042 O.lh4480896552227 0.164480893640937 0.164480893688676 0.164480893698442 0.164480893699234 0.1644X0893699288 0.164480893699292 0.164480893699293

0.164115002453607

0.164480893699293 0.164480893699293

O.l6448OX93699293 0.164480893699293 0.164480893699293

0.164480893699293

0.164480893699293

0.164480

893699293

U.D. Jentschuru

Table 7 Evaluation

of IO-’

Li3(0.99999)

et al. /Computer

Physic.7 Communications

I16 (1999)

28-54

with the CNCT

St?

‘f~"'(I,S(,)

6;"'

1 2 3 4 5 6 7 8 9 10 II 12 13 14 1s

0.13333133341s539 0.116665 166707769 0.121603216117468 0.1 I9 520 007 764 208 0.120 586 594 446 594 0. I 19 969 366 038 597 0.120358052 152572 0.120097666726411 0. I20 280 540 99 I 745 0. I20 147 227 650 952 0.120247 386435540 0.120 170239824481 0.1202309168138S7 0.120 182336 147846 0.120221 833436597 0.120 189289 I61 290

0.133331333415539 0.120474532 168000 0.120176326936846 0.120206042152677 0.120 203 955 328 380 0.120204046033711 0. I20 204 045 72.5 809 0.120204045413 120 0. I20 204 04s 439 707 0. I20 204 045 438 749 0. I20 204 045 438 729 0.120204045438733 0. I20 204 045 438 733 0. I20 204 04.5 438 733 0. I20 204 045 438 733 0. I20 204 045 438 733

0.133331333415539 0.120474532 168000 0. I20 176 326 936 846 0.120 204 748 497 388 0.120204079128106 0. I20 204 045 387 208 0. I20 204 045 378 284 0.120 204 045 434 802 0.120 204 045 438 553 0.120204045438726 0. I20 204 045 438 733 0.120204OJS438733 0.120 204 04.5 438 733 0. I20 204 045 438 733 0. I20 204 045 438 733 0.120204045438733

exact

0.120204045438733

0. I20 204 045 438 733

0.120 204 045 438 733

8 9 IO II 12 13 14 IS 16 17 18 I9 20

1.152086970 131424 O.S7609348506S712 0.960 055 803 546 36 I 0.672 109OSO5lS 104 0.902 446 45.5 72 1367 0.7lOSlS275487446 0.875 0 13 443 138 958 0.731090 03s 137 739 0.859010851 725411 0.743 892 107 I47 896 0.848536431 310 I59 0.752620 788 733 222 0.841 150625351432 0.758 951478 S83 304 0.835 664 028 102 224 0.763 752 250 680 046 0.831428053244539 0.767517561 OS3 133 0.828 0.59 092 035 324 0.770 549 625 376 190 0.825315796529872

1.152086970131424 0.806478876912452 0.797618 192 129 198 0.798663645011412 0.798 58 I 028 864 897 0.798 585 227 987 408 0.798 585 145 208 936 0.798 585 138 553 527 0.798 585 139 276 063 0.798585 139218618 0.798585 139219617 0.798585 139223310 0.798 585 I39 222 720 0.798 585 139 222 444 0.798 585 I39 222 555 0.798585 139222555 0.798 585 139 222 546 0.798 585 139 222 548 0.798 585 I39 222 548 0.798 585 I39 222 548 0.798 585 I39 222 548

I.152086970 131424 0.806478876912452 0.797618 192 129 198 0.798596 144946064 0.798 586 253 716 867 0.798585 188634 170 0.798 585 140 857 888 0.798 585 139 249 075 0.798585 139237667 0.798 585 I39 229 500 0.798 585 139 222 908 0.798 585 I39 222 I75 0.798585 139222491 0.798 585 I39 222 559 0.798 585 I39 222 550 0.798 585 I39 222 548 0.798 585 I39 222 548 0.798 585 I39 222 548 0.798 585 I39 222 548 0.798 585 I39 222 548 0.798 585 139 222 548

exact

0.798 585 I39 222 548

0.798 585 139222 548

0.798 58.5 139 222 548

!I 0

Table 8 Evaluation

0

4 5 6

of lo4 a(O.99999.2,

10000)

( I, S(j)

with the CNCT

41

U.D. Jentschura

42

et al. /Computer

Ph,ysics Cmmunications

The fact that the mathematical properties of the hypergeometric function 2Fl are so well understood, greatly facilitates its computation. The series ( 1.2) does not suffice for the computation of a nonterminating 2Fl since it converges only if ]Z 1 < 1. By contrast, the hypergeometric function 2Fj is a multivalued function defined in the whole complex plane with branch points at z = I and 00. However, a hypergeometric function 2FI (a, 6; c; z ) can be transformed into the sum of two other ~FI’s withargumentsw=l-z,w=l/z,w=l/(l-z), or w = I - 1/z, respectively. Thus, the argument w of the two resulting hypergeometric series can normally be chosen in such a way that the two new series in vv either converge, if the original series in z diverges, or they converge more rapidly if the original series converges too slowly to be numerically useful. For example, if 11- z 1 < 1 and if c - a - b is not a positive or negative integer, then we can use the analytic continuation formula (Eq. ( 15.3.6) of Ref. [ 341) I-(c)l-(c - a - b) I-(c - a)l-(c - 6) x2F,(a,b;u+b-c+ 1;l -z) l’(c)r(u+ b - c) (1 _ Z)c-a-h + r(a)r(b) xZFl(c-qc-b;c-a--b+];1

116 (1999)

28-54

computation of a generalized hypergeometric function p+t Fp with p > 2 via its series expansion ( 1.3) is still a more or less open problem. For example, in Theorem 2 of Ref. [ 801 explicit expressions for the analytic continuation formulas of a ,,+I F,, with p = 2,3, . . . around z = 1 were constructed. However, these expansions are by some orders of magnitude more complicated than the analogous formula (6.1) for a IF,. We show here that the CNCT can be very useful if the argument z of a generalized hypergeometric series P+.tFp is only slightly smaller than one, even if z = 1 is a singularity. The function 3F2( 1,3/2,5; 9/8,47/8; z ) has a singularity at z = I because the sum of the real parts of the numerator parameters is greater than the sum of the real parts of the denominator parameters (see p. 45 of Ref. [ 741). Thus, in Table 9 we evaluate the function IO-’ 3F2 ( 1,3/2,5;

9/8,47/8:

0.99999)

= 0.238 434 298 763 330

2F, (a, b;c; z) =

-iI.

(6.1)

Obviously, the two hypergeometric series in 1 - z will converge rapidly in the vicinity of z = 1. With the help of this or other analytic continuation formulas it is possible to compute a hypergeometric function 2F1 (a, b; c; z) with arbitrary real argument z E ( -00, +x) effectively via the resulting hypergeometric series (see p. 127 of Ref. [ 7.51 or Table I of Ref. [76]). The situation is much less favourable in the case of the generalized hypergeometric series

O” (a)n,(b)nr(c)nrm! z”’ 3 (6.2) .?F&b, c;d, e;z) = 5 (d),,,(e)n, and it gets progressively worse in the case of higher generalized hypergeometric series ,,+I Fp with p 2 3. The analytic continuation formulas of the type of Eq. (6.1) are not always known, and those that are known become more and more complicated with increasing p [77-801. Thus, the efficient and reliable

(6.3)

in the immediate vicinity of the singularity. With the help of the CNCT, it is nevertheless possible to evaluate this function without any noticeable loss of significant digits. The function 3F2( 1,3,7; 5/2, 14; Z) does not have a singularity at z = 1, i.e., the hypergeometric series (6.2) converges for z = 1, albeit very slowly. In Table 10 we compute the function IO-’ 3F2( 1,3,7; 5/2,14;0.99999) = 0.267 102 823 984 762

(6.4)

with the help of the CNCT. There is another major difference between a Gaussian hypergeometric series ~FI and a generalized hy. pergeometrtc series ,,.+I Fp with p 1 2. The series ( 1.2) for a hypergeometric function 2Fr (a, b; c; z ) with unit argument z = 1 converges if Re(c) > Re(u + b), and its value is given by the Gauss summation theorem (compare for instance Section 1.7 of Ref. [74]), 2F,(u,b;c;

1) =

r(C)r(C

-

U -

r(c - u)r(c -

6) 6)

(6.5)

The series (1.3) for ,)+I F,,(cYI, . . . , a~~+~;p,, . . ., P,,; z) with p L 2 converges at unit argument z = I

U.D. Jentschum Table 9 Evaluation ,l

of low4

~F2(1.3/2,5;9/8,47/8;0.99999)

et al. /Computer

with

Physics

Communicutions

I16 (I 999) 28-54

the CNCT tl~"'(l,so)

S" 195881

0.343961

s:"'

0.343961195

1 2 3

0.172030597S97940 0.286614046845766 0.200705485068072

0.240789533227428 0.238145631122015 0.238457646856530

0.240789533227428 0.238145631122015 0.238437505168542

4 5 6

0.269409131416317 0.212 175660263795 0.261216282820282

0.238433043073649 0.238434334888202 0.238434307597861

0.238434595265258 0.238434314202617 0.238434305531575

7 8 9

0.218319995918563 0.256437590806900 0.222142795720768

0.238434299034546 0.238434297734183 0.238434298685404

0.238434301380046 0.238434298863521 0.238434298539687

10 1I I?

0.253309971367287 0.224749013245827 0.251 104853004796

0.238434298885696 0.238434298760156 0.238434298751345

0.238434298711196 0.238434298769407 0.238434298766552

13 14 1s 16

0.226638970976291 0.249467024022419 0.228071952328669 0.248202731012486

0.238434298765877 0.238434298763853 0.238434298762980 0.238434298763382

0.238434298763311 0.238434298763230 0.23843429876333i 0.238434298763332

17 18 19

0.229 195680986673 0.247197365060 0.230100443575937

0.238434298763343 0.238434298763322 0.238434298763332

0.238434298763330 0.238434298763330 0.238434298763330

20

0.246378830994286

0.238434298763330

0.238434298763330

exact

0.238434298763330

0.238434298763330

0.238434298763330

if Re(/3r + . . + pr,) > Re( cyr + . + (Y,,.+.I) (compare p. 45 of Ref. [ 741) In Theorem 1 of Ref. [SO], an analogue of the Gauss summation theorem for a generalized hypergeometric series ,,+I F,, with unit argument was formulated. However, this expression is not a simple ratio of gamma functions as in Eq. (6.5), but an infinite series with complicated terms. Otherwise, there is a variety of different simple summation theorems for generalized hypergeometric series r,+.rFI, with unit argument which are ratios of gammafunction. However, thesesummationtheorems dependupon the value of p, and are usually only valid for certain values and/or combinations of the parameters cyr, ~2, . . ., CY,, ,.I and pr, &, . . ., pr,. A list of thesesimplesummationtheoremscan be found in Appendix III of Slater’s book [ 741. The reconstruction of known simple summation theorems for nonterminating generalized hypergeometric seriesand the construction of new oneswith the help of computer algebra systemsis discussedin books by PetkovSek, Wilf, and Zeilberger [ 8 1] and Koepf [ 82 1. With the help of Watson’s summation theorem (p. 245 Ref. [ 74]), we obtain

195 195 881

( 1. so)

0

I90

43

0.343961195195881

IO-‘&(1,3,7;5/2,

14;l) = ;o;7;;;;

=0.267108047538428.

(6.6)

We show in Table 11 that this result can also be obtained with the help of the CNCT, starting from the hypergeometric series(6.2). 7. Products of Besseland Hankel functions In this section, we want to apply the CNCT to series whose terms involve more complicated entities and which are more typical of the problems encountered in theoretical physics. As an example we consider here a series whose terms are products of spherical Bessel and Hankel functions (p. 435 of Ref. [ 34]),

j,(z) = [74(2z)l"2J1+1,2(z) p(z)

= ~?r/(2Z)l"Qf&~Z).

I

(7.1) (7.2)

Here, the index I is a nonnegative integer, and J~_.,,Q and H\!i12 are Besseland Hankel functions, respec-

44

U.D. Jentschura

Table 10 Evaluation

of IO- ’ 3 F2 ( 1,3,7;

5/2, 14; 0.99999)

et al. /Computer

Physics

Communications

I16 (I 999) 28-54

with the CNCT

n

S*

dp

( I, S”)

0 I 2 3 4 5 6 7 a 9 10 II I2 13 14 15 16 17 18 19 20

0.354205299 194014 0.227 102 649 597 007 0.288 360 357 978 268 0.254 808 733 179 764 0.274 632 798 536 574 0.262289292919201 0.270285 887 617646 0.264938303793252 0.268 610033 794438 0.266 03 1550 394 689 0.267878 112865 182 0.266 532 664 737 038 0.267 528 211474 989 0.266 78 I 286 870 185 0.267 348 762 560 I30 0.266912 657 747 879 0.267 25 1339 082 039 0.266985 765 115 182 0.267 195 878 973 722 0.267028262510 1.52 0.267 163 009 854 165

0.354205 299 194014 0.268 438 401594 236 0.266 943 489 834494 0.267 12 I 290 068 426 0.267 100 105730057 0.267 103 299 447 369 0.267 102742 191670 0.267 102 836 668 825 0.267 102 82i 243 079 0.267 102824 198445 0.267 102 823 960 840 0.267 102 823 987 269 0.267 102 823 984 5 IO 0.267 102 823 984 786 0.267 102 823 984 759 0.267 102 823 984 762 0.267 102 823 984 762 0.267 102 823 984 761 0.267 102 823 984 762 0.267 IO2 823 984 762 0.267 102 823 984 762

0.354 0.268 0.266 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267

exact

0.267 102 823 984 762

0.267 IO2 823 984 762

0.267 102 823 984 762

tively (pp. 358 and 360 of Ref. [34] ). In physical applications, the index 1 usually finds a natural interpretation as an angular momentum quantum number. Here, we investigate the following model problem:

@’

(I, so) 205 299 194 014 438 401594 236 943 489 834494 112224310036 101 775442210 102 948 107 378 102815 189 177 102 824 285 564 102 823 985 155 102 823 985 352 102 823 984 75 1 IO2 823 984 758 IO2 823 984 76 1 IO2 823 984 762 IO2 823 984 762 102 823 984 762 IO2 823 984 76 1 102 823 984 762 102 823 984 762 102 823 984 762 IO? 823 984 762

exp( -YW) = (xy)-‘!2 W m X

(2if

c

1)~/(cos~)~1+1/2(YX)~/+1/2(Y.V)

3

I=0

(7.4)

exp(-y[l - rl) = - O”(21+ ])j,(iry)/p(iy) , y;Il -rl

c I=0

(7.3)

where y is positive and 0 < r < 1. The spherical Besseland Hankel functions jl (iry) and hi” (iy) can alsobe expressedin termsof modijed sphericalBessel functions. However, we prefer to retain the given notation becausethere are some inconsistenciesin the literature regarding the prefactors to be assignedto the modi$ed spherical Bessel functions, whereasthe sphericalBesselandHankel functions aredefined consistently in most textbooks. With the help of known properties of Besselfunctions and Legendrepolynomials it can be showneasily that the seriesexpansion (7.3) is a specialcaseof the well known addition theorem of the Yukawa potential (p. 107 of Ref. [65]),

where w = (x2 + y2 - 2xy cosc$)“‘, 0 < x < y and y > 0. This addition theorem is also the Green function of the three-dimensionalmodified Helmholtz equation. First, we want to analyze the behaviour of the terms of the series on the right-hand side of Eq. (7.3) if the index I becomeslarge. The leading orders of the asymptotic expansions of the spherical Bessel and Hankel functions j,( iry) and h\” (iy) as 1 -+ oo with r and y fixed can be obtained easily from the series expansion for j/ (Eq. ( 10.1.2) of Ref. [ 341 or Eq. (E.17) of Ref. [ 831) and from the explicit expressionfor hj” (Eqs. (10.1.3) and (10.1.16) of Ref. [ 341 or Eq. (E. 18) of Ref. [ 831)) yielding

jr(h)

(iv) N (21+ I)!!

[1+0(1/l>]

3

(7.5)

U.D. Jentschura Table 11 Evaluation

of

10-l

3F2( 1.3,7;

S/2,14;

et al. /Compurer

Physics

Communications

S”

4”

0 1 2 3 4 5 6 7 8 9 I# II 12 13 14 15 16 17 18 19 20

0.354212896979703 0.227 106448 489 852 0.288366524620961 0.254 8 13 300 376 035 0.274638493612452 0.262294 169 832612 0.270291 370710880 0.264943330016988 0.268615409392370 0.266036 655 402 122 0.267 883 429 838 93 1 0.266 537 8 13 952 027 0.267 533 494 7 13 522 0.266 786 462 107 999 0.267354025525918 0.266917848923282 0.267256589412242 0.266 990 966 387 224 0.267201 121 176708 0.267 033 470 369 207 0.267168246683931

0.354 0.268 0.266 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267

exact

0.267 108 047 538 428

0.267 108 047 538 428

N -i

(21- I)!!e--J (iy)‘+’

[1 +0(1/l)]

.

( 1 , so)

(7.6)

Thus, we obtain for the leading order of the asymptotic expansion (1 -+ 00) of the product of the spherical Bessel and the spherical Hankel function, r/e-!

j,(iry)h,(‘)(iy)

N -

(21-l- l)Y

X3-54

45

1 f with the CNCT

II

h;“(iy)

116 (1999)

[’ +0(1/01

(7.7)

This asymptotic result showsthat the infinite serieson the right-hand side of Eq. (7.3) converges only for /F( < l,andthatthefunctionexp(-y[l-r])/(yllr] ) has a singularity at r = 1. Thus, it is interesting to usethe CNCT for an evaluation of the serieson the right-hand side of Eq. (7.3) in the immediate vicinity of the singularity. This series is similar to the sum over angular momenta encountered in QED bound statecalculations (cf. Ref. [30] and Eq. (4.2) in Ref. [ 831). In principle, the spherical Besseland Hankel functions in the products j/(iry) h!” (iy) can be evaluated recursively. If, however, the CNCT is used, this is not useful. The products belonging to all angular

,‘)

212 896 979 703 443 680 394 043 948 705 538 902 126514679686 105 329 111381 108523033 192 107 965 739 446 108 060 223 478 108 045 796 597 108047752 131 108047S14SO3 108 047 540 935 108 047 538 176 108 047 538 452 108 047 538 425 108 047 538 428 108 047 538 428 108 047 538 428 108 047 538 428 108 047 538 428 108 047 538 428

0.354 0.268 0.266 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267 0.267

( 1 . so) 212 896 979 703 443 680 394 043 948 705 538 902 117448402 341 106998932606 108 171668 587 108 038 742 29 I 108 047 839 255 108 047 538 82 1 108047539018 108047538417 108 047 538 424 108 047 538 427 108 047 538 428 108 047 538 428 108 047 538 428 108047538428 108 047 538 428 108 047 538 428 108047538428 108 047 S38 428

0.267 108 047 538 428

momenta 1are not needed. Instead, only the products for some speci$c angular momenta have to be evaluated. Accordingly, we prefer to evaluate the spherical Bessel function jl (iry) from its series expansion (Eq. (10.1.2) of Ref. [34] or Eq. (E.17) of Ref. [ 831) and the sphericalHankel function h!” (iy) via its explicit expression(Eq. ( 10.1.16) of Ref. [ 341 or Eq. (E.18) of Ref. 1831). It should be taken into account that the recursive evaluation of the terms of angular momentum decompositionsis normally possibleonly in the case of simple model problems. We shortly outline the applications of the CNCT to a more involved computational problem occurring in QED bound state calculations [ 301. The seriesto be evaluated in the QED calculation is given by Eq. (4.2) in Ref. [ 831. This seriescan be rewritten as

where J( r, y, t, y) has a simple mathematical structure. The term r,(r,y,t,~) (K = 1,2,.. .) can be rewritten as (see Eq. (5.7) in Ref. [ 841)

46

U.D.

Jentschura

et al. /Cotnprter

Physics

Communications

116

(1999)

28-54

a = e*/(4 z-). The series(7.8) occurs in the evaluation of an energy shift (the Lamb shift) of atomic levels due to the self energy of the electron. The self energy is an effect describedby second-orderperturbation theory within the framework of QED. The asymptotic behaviour of the terms T, in Eq. (7.8))

TKW~[ilO(l;)], (7.9) where the functions fi are the radial Dirac-Coulomb wave functions (see E$. (A.8) in Ref. [84]), l the functions Gi K are related to the radial components of the relativistic Green function of the bound electron (see Eq. (5.5) in Ref. [ 84]), and l the functions A, and A! are related to the angular momentum decomposition of the Green function of the virtual photon and are defined in Eq. (5.8) of Ref. [84]. The quantity a is a scaling variable for the subsequent radial integration over y (see Eq. (4.1) in Ref. [ 831) The variable r denotes the ratio of the two radial coordinates (0 < Y < I), and t is related to the (complex) energy of the virtual photon (0 < t < 1). The dependence on the coupling y of the electron to the central field, which appears on the left-hand side of Eq. (7.9), is implicitly contained in the functions Gi,, on the right-hand side of Eq. (7.9) (see Eq. (A. 16) in Ref. [ 841) . The propagator GB is defiqed in Eq. (5.1) in Ref. [ 841; the radial components G& of this propagator are explicitly given in Eqs. (5.4), (A.16) and (A.20) of Ref. [ 841 in terms of Bessel functions and Whittaker functions with non-integer parameters. Recurrence formulas relating the Gi,, for different values of K are not known. Therefore, the terms in the series (7.8) cannot be evaluated by recursion, but each one of these terms can be computed using techniques described in Refs. [ 30,831. Note that the angular momentum decomposition of Eq. (7.8) does not correspondto an expansion in terms of the QED perturbation theory. The perturbative parameter in QED is the elementary charge P or the fine structure constanth l

‘We use natural units as it is customary calculations ( fi = c = q = 1).

for QED

bound

state

Kim,

(7.10)

is similar to that of the model series(7.7) and to that of the Lil -seriesdiscussedin Section 5 (see Eq. (5.5) and the results presentedin Table 5). Slow convergenceof the seriesin Eq. (7.8) is observedin the limit r+ 1. Returning to our model problem given in Eq. (7.3), we consider here the example of r = 0.9999 and y = 0.7. This yields low5 5(21+ l)j,(iO.9999 X 0.7)hj”(i0.7) l=O = -0.142847 143207 13.5.

(7.1 I)

Evaluation of the seriesover the spherical Besseland Hankel functions via the CNCT is presented in Table 12.If this seriesis evaluatedby adding up its terms, then about 450000 products of Bessel and Hankel functions have to be evaluated for a relative accuracy of IO-l6 in the final result. This is to be compared with the approximately 300 evaluations of products of Besseland Hankel functions which are neededto compute the 25 Van Wijngaarden terms presentedin Table 12 with the necessaryaccuracy. According to our experience, the use of the CNCT in QED bound state calculations reducesthe amount of computer time for slowly convergent seriesof the type of Eq. (7.8) by three orders of magnitude. 8. Summary and conclusions Our approach for the acceleration of the convergenceof monotoneseriesconsistsof two steps.In the first step, the Van Wijngaarden condensationtransformation, Eqs. (2.2)-(2.4), is applied to the original monotone series.In the secondstep, the convergence of the resulting alternating series (2.3) is accelerated with the help of suitable nonlinear sequencetransformations.

Table 12 Evaluation

of 10es cE,(21+

U.D. Jentschura

et al. /Computer

1) jt(i0.9999

x 0.7) Iz~“(iO.7)

Physics

Communications

116 (I999)

28-54

with the CNCT

n

S”

4t0) ( 1, so)

@’ ( 1, so)

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 1.5 16 17 18 19 20 21 22 23 24 2s

-0.206 084 520 894 668 -0.103 046 104 279 554 -0.17 1733 352 076 042 -0.120 220 473 724 774 -0.161427 969 242 195 -0.127091 187739732 -0.156520814934640 -0.130 771367 274 807 -0.153 658 28 1433 404 -0.133061584 756825 -0.151784408324881 -0.134 623 097 855 864 -0.150 463 209 08 1729 -0.135 755 492 037 282 -0.149481 826233342 -0.136614208810295 -0.148 724 109 578 581 -0.137287 765 215 359 -0.148 121427644456 -0.137830196213572 -0.147630649333917 -0.138 276 357 305 065 -0.147223291 777 104 -0.138649758 251940 -0.146 879 773 958 865 -0.138966841 391919

-0.206 084 520 894 668 -0.144 259 660 09 I 669 -0.142 674 25 1 704 499 -0.142 861066 165 942 -0.142 846 281288 224 -0.142847 117647583 -0.142 847 153435 206 -0.142 847 152 002 732 -0.142 847 142 941 523 -0.142 847 142 048 135 -0.142 847 143 286 397 -0.142 847 143 324 758 -0.142 847 143 I8 I 923 -0.142 847 143 200611 -0.142 847 143 210 852 -0.142 847 143 206780 -0.142 847 143 206921 -0.142 847 143 207 226 -0.142 847 I43 207 123 -0.142 847 143 207 132 -0.142 847 143 207 137 -0.142 847 143 207 134 -0.142847 143 207 135 -0.142 847 143 207 135 -0.142 847 143 207 135 -0.142 847 143 207 135

-0.206 084 520 894 668 -0.144 259 660 091669 -0.142 674 25 1704 499 -0.142 849 004 178 547 -0.142847217919030 -0.142 847 093 177 734 -0.142 847 I35 794 48 I -0.142 847 148 838 958 -0.142 847 14.5 983 994 -0.142 847 143 380 152 -0.142 847 143 026 940 -0.142 847 143 169466 -0.142847 143211999 -0.142847 143208928 -0.142 847 143 207 036 -0.142 847 143 207 092 -0.142 847 143 207 139 -0.142 847 143 207 135 -0.142 847 143 207 135 -0. I42 847 143 207 135 -0.142 847 143 207 135 -0.142 847 143 207 135 -0.142 847 143 207 135 -0.142 847 143 207 135 -0.142847 143207 135 -0.142 847 143 207 135

exact

-0.142847

-0.142

-0. I42 847 143 207 135

143207

135

Daniel [ 211 showed that the transformed alternating series (2.3) converges under relatively mild conditions to the same limit as the original series. The condensation transformation cannot solve the problems due to slow convergence on its own since the transformed alternating series does not converge more rapidly than the original series (see Eq. (3.2) ) . However, the convergence of alternating series can be accelerated much more effectively than the convergence of the original monotone series. Moreover, the transformation of alternating series by sequence transformations is numerically a relatively stable process, whereas the direct transformation of monotone series inevitably leads to the loss of accuracy due to round-off. Many sequence transformations are known which are in principle capable of accelerating the convergence of alternating series [ 1,2,9]. However, their ef-

41

847 143 207 135

ficiency varies considerably. In all numerical examples studied in the context of this paper, we found that certain sequence transformations [ 9,121 (which also use explicit estimates for the truncation errors as input data) are significantly more powerful than other, better known sequence transformations as for instance Wynn’s epsilon algorithm [ 111 or the Euler transformation [ 351 (which only use the partial sums of the infinite series to be transformed as input data). Consequently, we use in this article only the sequence transformations G:“’ (P, s,, w,) and Si”’ (,f3, s,, w,) defined in Eqs. (3.8) and (3.9) in combination with the truncation error estimate (3.12). Those properties of these sequence transformations, which are needed to apply them successfully in a convergence acceleration process, are described in Section 3. In Appendix A. we explain why these transformations are much more effective in the case of alternating series than in the

48

U.D.

Jentschura

et al. /Computer

Physics

case of monotone series, and in Appendix B we discuss exactness properties of these transformations. We consistently observed that the combined transformation is a remarkably stable numerical process and entails virtually no loss of numerical significance at intermediate stages. The evaluation of the condensed series (2.4) is a computationally simple task, provided that the evaluation of terms of sufficiently high indices in the original series is feasible. The nonlinear sequence transformations used in the second step are also remarkably stable. The final results, which we present in the tables, indicate that a relative accuracy of lo-l4 can be achieved with floating point arithmetic of 16 decimal digits. Moreover, our method remains stable in higher transformation orders. In Sections 4-6 we demonstrate that our approach is very useful for the computation of special functions which are defined by logarithmically convergent series. Our approach also works very well if the argument of a power series with monotone coefficients is very close to the boundary of the circle of convergence (even in the vicinity of singularities) I In Section 4, we treat the Dirichlet series ( 1.4) for the Riemann zeta function which because of its simplicity is excellently suited to illustrate the mechanism of our combined transformation. The application of the Van Wijngaarden transformation to the logarithmically convergent series ( 1.4) yields the known alternating series (4.3) which does not converge more rapidly than the original monotone series ( 1.4). However, the alternating series can be transformed effectively by the sequence transformations (3.13) or (3.14) in the case of slow convergence or even summed in the case of divergence. Consequently, the Riemann zeta function can be evaluated effectively and reliably even if its argument z is only slightly larger than one (Table 2)) if both the monotone series ( 1.4) and the alternating series (4.3) diverge (Table 3)) or if the argument z is complex (Table 4). In Section 5 we treat the Lerch transcendent 9( Z, s, a) which contains many other special functions as special cases, for example the Riemann zeta function or the polylogarithms. In Tables 5-8 we show that the combined transformation makes it possible to evaluate @J(z, s, au) effectively and reliable via the power series (5.1) even if z is very close to the boundary of the circle of convergence. In Section 6 we discuss the evaluation of a non-

Communications

116

(1999)

28-54

terminating generalized hypergeometric series ,,+i Fp via its defining power series (1.3) which converges for ]Z 1 < 1. In the case of a Gaussian hypergeometric series ~FI, explicit analytic continuation formulas are known which transform a hypergeometric series with argument z into the sum of two hypergeometric series 2Fi with argument 1 - 2. Thus, if the convergence of a hypergeometric series 2Fi is slow because its argument z is only slightly smaller then one, then the two transformed hypergeometric series 2FI with argument 1 - z will converge rapidly. Moreover, if a Gaussian hypergeometric series zF1 with unit argument z = 1 converges, then its value is given by the Gauss summation theorem (6.5). In the case of a generalized hypergeometric series ,,+l Fp with p 2 2, the situation is much more difficult. Explicit analytic continuation formulas are either unknown or they become increasingly complicated with increasing p. Similarly, simple analogues for the Gauss summation theorem (6.5) exist only for certain values of p and for certain combinations of the parameters of the p+i F,,. In Tables 9-l 1 we show that a generalized hypergeometric series ,]+, Flj (p 2 2) with an argument z that is only slightly smaller than one or with unit argument (Z = 1) can be computed effectively and reliably with the help of the combined transformation. Partial wave decompositions of Green’s functions, which occur for example in quantum electrodynamic bound state calculations, entail more complex mathematical entities than series expansions for special functions. Nevertheless, we show in Section 7 that the combined transformation can be applied successfully to these partial wave decompositions. The series over products of spherical Bessel and Hankel functions considered in this paper serves as a model problem for the angular momentum decomposition of more complex Green’s functions, as for example the relativistic Green’s function of the bound electron [X4]. In the model problem studied here as well as in a recent evaluation of self energy corrections in bound systems [ 301, we observed a reduction in computer time by three orders of magnitude. Thus, the combined transformations makes extensive and highly accurate calculations feasible in situations that could otherwise be too time-consuming.

V.D. Jentsrhura

et al. /Computer

Physics

Acknowledgements We thank Prof. F.W.J. Olver for illuminating discussions. U.D.J. thanks J. Baker, J. Devaney and J. Sims for stimulating discussions, and acknowledges support from the Deutsche Forschungsgemeinschaft (contract no. SO333-l/2) and the DAAD. P.J.M. acknowledges continued support by the Alexander von Humboldt Foundation. G.S. thanks the Deutsche Forschungsgemeinschaft and the Gesellschaft fur Schwerionenforschung for support. E.J.W. thanks Prof. N. Temme for providing information on the Van Wijngaarden transformation and Prof. A.Z. Msezane and Prof. CR Handy for their invitation to the Center for Theoretical Studies of Physical Systems at Clark Atlanta University. E.J.W. acknowledges support from the Fonds der Chemischen Industrie.

Communications

116 (1999)

49

28-54

Loosely speaking, this means that we have to choose the remainder estimates in such a way that the numerator of the ratio on the right-hand side of Eq. (A.2) becomes as small as possible and the denominator becomes as large as possible. If we can find remainder estimates such that S,-s=co”[C+u(n-‘)],

n+co,

(A.3)

then the weighted difference operator Akwk( n) will make the numerator small, no matter whether the input data {sn}go are the partial sums of an alternating or of a monotone series. The situation is quite different in the case of the denominator. Application of the weighted difference operator Akwk( n) to 1/o,, yields

Appendix A. On the efficiency of the transformation of alternating and monotone series In the case of the Levin transformation ckn’ (p, s,, w,), Eq. (3.8), and the closely related Weniger transformation Si”’ (p, s,,, w,), Eq. (3.9), it is relatively easy to understand why alternating series can be transformed much more effectively than monotone series. These two transformations are both special cases of the following transformation which is characterized by the positive weights wk (n):

?;(“yWk(n);Sn,oJ,)= Ak{w(n>sn/wn) Ak{w(n)/‘w,) ’

(A.1)

If wk (n) = (n +p)‘I, we obtain Levin’s transformation, and if wk (n) = (n + @)k- t, we obtain Weniger’s transformation. Since the difference operator A’ is linear, Eq. (A. 1) can also be rewritten as follows:

= s + Ak{Wk(n) [&I - sl/%} Ak{W/%>

’

(A.2)

Obviously, the sequence transformation Ii”’ converges to the (generalized) limit s of the input sequence {~,,}:e if the remainder estimates (wn}F? can be chosen in such a way that the ratio on the right-hand side becomes negligibly small.

If the remainder estimatesall have the samesign, the denominatorwill alsobecomerelatively smallbecause of cancellation due to differencing (hopefully not as small as the numerator becausethen the transformation processwould not converge). This cancellation processis also the major source of numerical instabilities which are magnified since they occur in the denominator. If the remainderestimatesare strictly alternating in sign, o, = (-l)“lw,l, then we obtain (wk(n) is by assumptionpositive)

I (A.5)

Thus, all termsof the denominator sum have the same sign, and there is no cancellation as in the case of monotoneremainder estimatesand also no source of numerical instabilities.Consequently,the denominator becomesrelatively large which improvesconvergence. As an example, we apply the transformation (A. 1) to the monotone mode1sequence

co sn=s+ (n+l)”

+

(n +cl,a+l

t...

,

as well as to the alternating model sequence

(A.6)

SO

U.D.

Jentschura

et al. /Computer

Physics

t,,=t+ (-1)” (n$a {

+(n+c;)“+l +...

and derive asymptotic ( IZ + 00) transformation error estimates. For that purpose, we assume cy > 0 - which implies that the two model sequences converge to their limits s and t, respectively - and wk(n) = U(nk-‘) as n -+ cc. In the case of the monotone input sequence {s,,):~, we choose w,, = (n + 1)-O, and in the case of the alternating input sequence {t,,}gO, we choose QI, = ( - 1)“( n-t 1) --n. Then, we obtain for the numerators

*k

wktn) [s, - sl (n + 1)-a

{ = A’

>

wktn) lfn - tl (-l)n(Yz

= u(n-k-‘).

+ I)--”

tA.8)

The estimates for the denominators differ considerably. In the case of the monotone input sequence (A.6), we obtain

(A.91 which in combination with Eq. (A.8) yields the following asymptotic transformation error estimate as II + 00: &(n+k(n);s,,

(n + I)-*) s, - s

- s

= U(n-k).

(A.10)

In the case of the alternating input sequence (A.7), we exploit the fact that the first term with j = 0 in the sum on the right-hand side of Eq. (AS) is smaller in magnitude than the whole sum, yielding 1 fAk{Witn)/[(-l)“tn

+ I>-“]}I 1

’

Iwk(n)/[(-l)“(n+

(A.1 1) I)-“11

’

Of course, we could also try to construct more sophisticated estimates for the denominators. However, the relatively crude estimate (A. 11) suffices for our purposes since it implies in combination with Eq. (A.8),

?;.‘“‘(Wk(~t);f,. (-l)“(n+ t, - t

I,-“)

-t

= u(n-2k). (A.12)

Communications

I16

(1999)

28-54

The two transformation error estimates (A.lO) and (A.12), which both hold as II --) 00, show that there is a substantial difference between the transformation of monotone and alternating series. There is a considerable amount of evidence that this conclusion is actually generally true, i.e., also in the case of Pad6 approximants and other sequence transformations. Appendix

B. Exactness

results

In this appendix, we analyze exactness properties of the Levin transformation f$“) (p, s,, o,), Eq. (3.8), and of the Weniger transformation Sk”’ (p, s,, a,,), Eq. (3.9), and their variants d:“’ (p, S,), Eq. (3.13), and Sy’(p, S,), Eq. (3.14), which both use the remainder estimate (3.12) proposed by Smith and Ford [37]. For that purpose we introduce in Eq. (A.2) the remainder r,, = s,, - s and obtain

Obviously, the general sequence transformation I>“‘, Eq. (A.1 ), which contains CL’) (p, s,, o,) and 3;“’ (p, s,, on) as special cases, is exact for a given input sequence {s,,}:~ if the difference operator Ak annihilates wk(rz)r,/w,, but not wk(n)/o,. Thus, if we want to prove the exactness of 5:“) for some sequence {s,,}zO, we need to know explicit expressions for the remainders {r,,}z$, of the input sequence as functions of n. In the case of the alternating series (4.3) for the Riemann zeta function with zero or negative integral argument, this can be accomplished easily. Replacing z by -I in Eq. (4.3) with I = 0, 1,2, _. ., we obtain

It-11=j&F

(-l).‘(j

+ 1)’

j=O

For 1 = 0, this series is the negative of the geometric series l/( 1 + X) = c,co( -n) j at x = 1. According to Theorem 12-9 of Ref. [ 91, the sequence transformations 13:“’ (/3, s,, w,,) and Sin’ ( p, s,, w,) are exact for the geometric series if the first term neglected in the partial sum is used as remainder estimate (which corresponds to the remainder estimate (3.12) ) . Thus,

U.D. Jentschuru

et ul. /Computer

Physics

the exactness of these sequence transformations for the infinite series (B.2) has to be analyzed only for 12 1. For the determination of its truncation error with arbitrary integral 1 2 1, we rewrite the infinite series (B.2) as follows:

5

(-l)‘(j .;=o

+ 1)’

= e(-l)j(j /=0

Crmtmunicutions

116 11999) 28-54

51

If we set in Eq. (B.8) n = - 1, then the partial sum is an empty sum and vanishes,yielding the following explicit expressionfor the Riemannzeta function with zero or negative integer argument,

5(-l)=&

[(~&+T] ljil. (B.9)

Equivalent explicit expressionscan be obtained via the substitution y = In(x) Obviously, we have

+ l)‘+

9

(-l),‘(j

+ 1)’

i=n+l

IIcII

(B.3)

(-l),‘(j

In this way, we obtain

+ 1)’

.i=O

00

c (-l)“(n+v+2)‘. v=o

+ (-l)n+’

(B.4)

The infinite series on the right-hand side of Eq. (B.4) obviously diverges. However, it can be summed easily since it can be represented as a derivative of the geometric series I/( 1 + X) = C,“f( -x)j at x = 1,

111. “=O

9 (-l)“(n+v+2)’ v=o

(-l)“(n

+ V + 2)‘Xn+v+’

(B.5)

(B. 12)

Combination of Eqs. (4.llb), (B.9), and (B.11) yields the following expressions for the Bernoulli numberswith even indices: (B.13)

(B.6) =fq

[;+q2’-1&J}~J=o.

CB.14)

(J3.7) By inserting this into Eq. (B.2)) we obtain the following representationof l( -1) as a partial sum plus an explicit expressionfor the truncation error,

5(-l> = j-J+S(-l)n+’

~(-l)‘(i 1 .i=o

[(is)’

Ei]

+ 1)’

l.r;,}.

(B*8)

It follows from Eq. (B. 1) that the general sequence transformation q”), Eq. (A.l), is exact for someinput sequence{sn}z+ if the expression wk( n)r,/w, with r, = s, - s is a polynomial of degreek - 1 in n. Thus, Eq. (B.7) implies that in the caseof the Weniger transformation 6:“) (1, Se), Eq. (3.14)) which usesthe remainder estimate (3.12), we have to show that the ratio (B.15)

52

U.D. Jentschura

et al. /Computer

Physics

is for sufficiently large values of k a polynomial of degree k - 1 in n. For 1 = 1, 2, 3, and 4, Eq. (B.7) yields the following explicit expressions:

cm cv=O(-l)V(nfv+2)i=v,

c c

Conununications

E!,Fo(lf

116 (I 999) 28-54

1; -1)

=e

(-l)%

+ 111

j=O

(B.16)

33 (-l)“(n+~+2)*= v=o

(n+1)2(n+2),

(B.17)

+ (-l)“+‘g

(B.22)

(-l)“(n+v+2)1

=k

(-l)‘(j + 1)1+ (-l)‘lf’(n .i=o x *FI(l,n+z+2;n+2;-1)

co

(-l)“(n+zJ+2j3

f2)/ (B.23)

Vdl

=

=k

(2n+3)(2n2+6n+3) 8

(B.18)

(-l>.‘(j

+ l)/

j=o

cc t3

(-l)“(n+v+2)4

+

C-l>"+'

+

2),&j( l/2)‘+‘. (B.24)

l)(n+2)(n2+3n+ 2

1)

(B.19)

Since(n+l)k-l=(n+l)(n+2)...(n+k-l), theratio(B.l5)isforI=1,2andk>3apolynomial of degree k - 1 in II, whereas for 1 = 3, 4 it is a rational expression in n which is not annihilated by the operator Ak. Thus, the exactness of the transformation Sj~’ (1, So) for l( - 1) as in Table 3, or for I( -2) is more or less accidental. If, however, we choose in the Levin transformation dh”) (/3, So), Eq. (3.13), with p = 2 instead of the usual /3 = 1, we find that the corresponding ratio (n+2>k-‘C~~(--l)v(.+.+2)’ (-l)“+‘(n + 2)’

(B.20)

is a polynomial of degree k - 1 in n for all I= 0, 1,2, . and k 2 I + 1. This follows at once from the fact that the differential operator in Eq. (B.7) produces a polynomial of degree 1 in n. In the same way, it can be shown that the Weniger transformation SA”) (p,So), Eq. (3.14), with p= 2 is forall/=O, 1,2,... and for k > 1 + 1 exact for the hypergeometric series l!1Fo(l+

+j

j=O

v=o = (n+

C(-l)j(n

1; -1)

= F(-l)‘(j

+ 1), = &,

j=O

(B.21) This can also be rewritten

as follows:

In this case, the corresponding (n+2)k-,C~~(-l)“(n+V+2), c-1),+1 (n + 2)/

ratio (B.25)

is a polynomial of degree k - 1 in n for all 1 = 0, 1, 2 . . and for k 2 1 + 1. This follows at once from tde fact that the second sum on the right-hand side of Eq. (B.24) is a polynomial of degree 1 in n. However, the value p = 2 in the Eqs. (3.13) and (3.14) cannot significantly improve convergence for the zeta function. We also investigated the performance of the transformations die) (2, SO) and c!$,‘) (2, So). Except for the cases of z = - 1 or z = -2 in which the Levin transformation becomes exact with p = 2, the transformations dA”’ (2, SO) and 6:‘) (2, So) yield results which differ only marginally from the results obtained by d,(,‘) (1, SO) and 6:‘) (1, SO), which are presented in Tables 2-4. References I 11 J. Wimp, Sequence Transformations and Their Applications (Academic Press, New York, 198 I). 121 C. Brezinski, M. Redivo Zaglia, Extrapolation Methods (North-Holland, Amsterdam, 1991). [ 31 C. Brezinski, History of Continued Fractions and Pad6 Approximants (Springer, Berlin, 199 1). 141 C. Brezinski, Appl. Numer. Math. 20 ( 1996) 299. 151 F.W.J. Olver. Asymptotics and Special Functions (Academic Press, New York, 1974).

U.D. Jentschura

et ~1. /Computer

Physics

[ 6 ] C.M. Bender, S.A. Orszag, Advanced Mathematical Methods for Scientists and Engineers (McGraw-Hill, New York. 1978). ] 7 ] H.M. Edwards, Riemann’s Zeta Function (Academic Press, New York, 1974). ] 8 ] B. Germain-Bonne, Rev. Francaise Automat. Informat. Rech. Operat. 7 (R-l) (1973) 84. ]9] E.J. Weniger. Comput. Phys. Rep. IO (1989) 189. [ IO] A.C. Aitken. Proc. Roy. Sot. Edinburgh 46 ( 1926) 289. [ I I ] I? Wynn, Math. Tables Aids Comput. IO ( 1956) 91. ] I2 ] D. Levitt, lnt. J. Comput. Math. B 3 (1973) 371. ] 131 J.P. Delahaye, B. Germain-Bonne, SIAM J. Nume?. Anal. 19 (1982) 840. ] 141 L.F. Richardson, Phil. Trans. Roy. Sot. London A 226 (1927) 229. ] IS ] P. Wynn, Proc. Camb. Phil. Sot. 52 ( 1956) 663. [ 161 N. Osada, SIAM J. Numer. Anal. 27 (1990) 178. [ 171 C. Brezinski, C.R. Acad. Sci. Paris 273 (1971) 727. ] I8 ] P. Bjorstad, G. Dahlquist, E. Grosse, BIT 21 (1981) 56. ] 191 C.W. Clenshaw, E.T. Goodwin, D.W. Martin, G.F. Miller. F.J.W. Olver, J.H. Wilkinson, eds., National Physical Laboratory, Modem Computing Methods, 2nd ed.. Notes on Applied Science, Vol. 16 (H.M. Stationary Office, London. 1961). [ 201 A. van Wijngaarden, Cursus: Wetenschappelijk Rekenen H. Process Analyse (Stichting Mathematisch Centrum, Amsterdam, 1965) pp. 51-60. 121 ] J.W. Daniel, Math. Comput. 23 ( 1969) 91. ] 22 ] P.J. Pelzl, F.E. King, Phys. Rev. E 57 ( 1998) 7268. 1231 C. Brezinski, J.P. De&aye, B. Germain-Bonne. SIAM J. Numer. Anal. 20 (1983) 1099. ] 24 ] Y.L. Luke. The Special Functions and Their Approximations I & II (Academic Press, New York, 1969). ] 25 ] Y.L. Luke, Algorithms for the Computation of Mathematical Functions (Academic Press, New York, 1977). ]26] C.G. van der Laan, N.M. Temme, Calculation of Special Functions: The Gamma Function, the Exponential Integrals and Error-Like Functions (Centrum voor Wiskunde en lnformatica, Amsterdam, 1980). ] 27 ] S. Zhang, J. Jin. Computation of Special Functions (Wiley. New York, 1996). [ 28 ] D.W. Lozier, F.W.J. Olver, in: Mathematics of Computation 1943-1993: A Half-Century of Computational Mathematics, W. Ciautschi, ed., Proc. Symp. Appl. Math. 48 (1994) 79. ]2Y] D.W. Lazier, J. Comput. Appl. Math. 66 (1996) 345. ] 30 ] U. Jentschura, P.J. Mohr, G. Soff, manuscript in preparation. ] 3 I ] S. Wolfram, Mathematics - A System for Doing Mathematics by Computer (Addison-Wesley, Reading, MA, 1988). ] 32 ] T.J. I’a. Bromwich. An Introduction to the Theory of Infinite Series, 3rd ed. (Chelseu, New York, 1991). ] 33 ] K. Knopp. Theorie und Anwendung der unendlichen Reihen (Springer, Berlin, 1964). ] 341 M. Abramowitz, I.A. Stegun. Handbook of Mathematical Functions, 10th printing (National Bureau of Standards, Washington, D.C., 1972). 135) L. Euler, Institutiones calculi differentialis cum eius usu in analysi finitorum ac doctrina serium. Part 11.1. De

Communications

[ 361 ] 37 ] ] 38 ]

] 39 ] 140) [41 ] (421 ]43] ]&I] (4.5 ] ]46] ]47] [ 481 ] 49 ] ] 50 ] ] Sl ] ] 521 [53] I.541 ]SS] ]S6] ]S7] [ 581 ]SY 1

] 601 161 ] [ 621 1631 ] 641 16.5 ]

] 661

I16 (1999)

28-54

53

transformatione serium (Academia Imperialis Scientiarum Petropolitana, 175.5). This book was reprinted as Vol. X of Leonardi Euleri Opera Omnia. Seria Prima (Teubner, Leipzig, Berlin, 1913). B.C. Carlson, Special Functions of Applied Mathematics (Academic Press, New York, 1977). D.A. Smith, W.F. Ford, SIAM J. Numer. Anal. 16 ( 1979) 223. E.J. Weniger, in: Nonlinear Numerical Methods and Rational Approximation II, A. Cuyt. ed. (Kluwer. Dordrecht, 1994) p. 269. H.H.H. Homeier. E.J. Weniger, Comput. Phys. Commun. 92 (1995) I. A. Sidi, Math. Comput. 33 (1979) 315. A. Sidi, Math. Comput. 35 ( 1980) 833. A. Sidi, SIAM J. Math. Anal. 17 (1986) 1222. D.A. Smith, W.F. Ford, Math. Comput. 38 (1982) 481. E.J. Weniger, J. Comput. Appl. Math. 32 (1990) 291. E.J. Weniger, J. &ek, Comput. Phys. Commun. 59 (1990) 471. E.J. Weniger, Comput. Phys. IO (1996) 496. D. Roy. R. Bhattacharya, S. Bhowmick. Comput. Phys. Commun. 93 ( 1996) 159. R. Bhattacharya, D. Roy. S. Bhowmick, Comput. Phys. Commun. IO1 (1997) 213. E.J. Weniger. J. &ek. F. Vinette, Phys. Lett. A 156 ( 199 1) 169. E.J. Weniger, Numer. Algor. 3 ( 1992) 477. E.J. Weniger, J. Ciiek, F. Vinette, J. Math. Phys. 34 ( 1993) 571. E.J. Weniger, bit. J. Quantum Chem. 57 (1996) 265; Erratum, Int. J. Quantum Chem. 58 (1996) 319. E.J. Weniger. Ann. Phys. (NY) 246 (1996) 133. E.J. Weniger, Phys. Rev. Lett. 77 (1996) 28.59. E.J. Weniger. Phys. Rev. A 56 (1997) 5165. E.C. Titchmarsh, The Theory of the Riemann Zeta-Function, 2nd ed. (Clarendon Press, Oxford, 1986). A. Ivic, The Riemann Zeta-Function (Wiley, New York, 1985). S.J. Patterson, An Introduction to the Theory of the Riemann Zeta-Function (Cambridge Univ. Press, Cambridge, 1988). E. Elizalde. S.D. Odintsov, A. Romeo, A.A. Bytsenko, S. Zerbini, Zeta Regularization Techniques with Applications (World Scientific, Singapore, 1994). E. Elizalde. Ten Physical Applications of Spectral Zeta Functions (Springer, Berlin, 1995). M.V. Berry, J.P. Keating, Proc. Roy. Sot. London A 437 (1992) 151. R.B. Paris, Proc. Roy. Sot. London A 446 (1994) 565. M.V. Berry. Proc. Roy. Sot. London A 450 ( 1995) 439. R.K. Bhaduri. A. Khare, S.M. Reimann, E.L. Tomusiak, Ann. Phys. 2.54 (1997) 25. W. Magnus, F. Oberhettinger, R.P. Soni, Formulas and Theorems for the Special Functions of Mathematical Physics (Springer, New York, 1966). L. Lewin, Polylogarithms and Associated Functions (NorthHolland, New York, 198 1 ).

54

U.D. Jentschuru

et al. /Computer

Physics

167 I A. Griffin, W.-C. Wu, S. Stringari, Phys. Rev. Lett. 78 ( 1997) 1838. 1681 A.J. MacLeod, Comput. Phys. 1 I (1997) 385. [ 69 I K. Pachucki, Phys. Rev. A 46 ( 1992) 648. [ 70 1 K. Pachucki, Ann. Phys. (NY) 226 ( 1993) I. [ 711 U. Jentschura, K. Pachucki, Phys. Rev. A 54 (1996) 1853. 1721 U.D. Jentschura, G. Soff, P.J. Mohr, Phys. Rev. A 56 ( 1997) 1739. (731 P. Indelicate, P.J. Mohr. Phys. Rev. A 58 (1998) 165. ( 74 I L.J. Slater, Generalized Hypergeometric Functions (Cambridge Univ Press, Cambridge, 1960). 1751 N.M. Temme, Special Functions - An Introduction to the Chxsical Functions of Mathematical Physics (Wiley. New

Communications

116 (1999)

28-54

York, 1996). R.C. Forrey, J. Comput. Phys. 137 (1997) 79. W. Biihring, SIAM J. Math. Anal. 18 (1987) 884. W. Biihring, SIAM J. Math. Anal. 18 (1987) 1227. W. Biihring, SIAM J. Math. Anal. 19 (1988) 1249. W. Biihring, Proc. Amer. Math. Sot. 114 (1992) 14.5. M. PetkovSek, H.S. WiIf. D. Zeilberger, A = J3 (A.K. Peters, Wellesley, MA, 1996). [ 821 W. Koepf, Hypergeometric Summation (Vieweg. Braunschweig, 1998). 1831 P.J. Mohr, Ann. Phys. (NY) 88 (1974) 52. [841 P.J. Mohr, Ann. Phys. (NY) 88 (1974) 26.

1761 [771 [781 [ 791 [SO1 [81I