Accelerating convergence for some double sums

Question

I am interested in the following general double sums, for integers $a\geq 1$ and $b\geq 2$,

$$Z(a,b) = \sum_{k,\ell \geq 0} \frac{2k+3}{\binom{k+2}{2}^a} \frac{2\ell+3}{(\binom{k+2}{2}+\binom{\ell+2}{2})^b},$$

which are converging very slowly. For these sums, there is also an alternative expression as an iterated integral in dimension $a+b$, similar to multiple zeta values.

I would like more particularly to find $Z(2,2)$ and $Z(1,3)$ with precision as large as possible. Using the naive summation, I could only obtain 4 decimal digits, namely $Z(2,2) \simeq 4.7058$ and $Z(1,3) \simeq 1.6470$.

It is known that $2 Z(2,2) + 4 Z(1,3) = 16$.

What would be a smart way to accelerate the convergence, in general and in the special case using maybe the previous formula ?

Answer 1

This is easy to do with PARI/GP. Here is my code

p(n) = binomial(n+2,2);
Y(k,b) = sumnum(l=0, (2*l+3)/(p(k)+p(l))^b,sumtable);
Z(a,b) = sumnum(k=0, (2*k+3)/p(k)^a*Y(k,b));
default(realprecision,57);
sumtable = sumnuminit();
print(2*Z(2,2)+4*Z(1,3))
/* 16.0000000000000000000000000000000000000000000000000000000 */

It takes 1.361 CPU seconds for 57 decimal places. The documentation has some information about the methods being used and there are several summation functions other than sumnum. For example, I originally used sumpos but sumnum is faster. Thanks to Henri Cohen for his comment to use sumnuminit to speed up the calculation of Y(k,b). Thanks to Jorge Zuniga for his comment that replacing sumnum with summonien is much faster. And especially thanks to the developers of PARI/GP who provided these fast numerical summation functions.

Some details are in the arXiv paper Gaussian Summation: An Exponentially Converging Summation Scheme by Hartmut Monien. (published in 2010 in Mathematics of Computation).

Answer 2

With Mathematica I can first sum the series over $\ell$ to get a closed-form expression in terms of polygamma functions, _{$$Z(2,2)= \sum_{k,\ell \geq 0} \frac{2k+3}{\binom{k+2}{2}^2} \frac{2\ell+3}{(\binom{k+2}{2}+\binom{\ell+2}{2})^2}$$ $$\qquad=\sum_{k\geq 0}\frac{16 (2 k+3) }{(k+1)^2 (k+2)^2 \sqrt{-4 k (k+3)-7}}\left(\psi ^{(1)}\left[\frac{1}{2} \left(3-\sqrt{-4 k (k+3)-7}\right)\right)-\psi ^{(1)}\left(\frac{1}{2} \left(\sqrt{-4 k (k+3)-7}+3\right)\right)\right],$$} and then obtain a high-precision result using the Wynn epsilon method.. To fifteen decimal places I find

$$Z(2,2)=4.705905644127174\cdots$$

NSum[Z22,{k,0,Infinity},NSumTerms->2000,WorkingPrecision->250,Method->"WynnEpsilon"]

As a test, I also computed $Z(1,3)$ to fifteen decimal places, which gave

$$Z(1,3)=1.647047177936412\cdots.$$

And thus I find

$$2Z(2,2)+4Z(1,3)=15.999999999999996\cdots$$ correct to fifteen decimal places.

Answer 3

Here is a little Maple. Note that using "sum" on the inside causes it to find an algebraic expression for the sum over $\ell$ and using "Sum" on the outside tells it to not try to sum that algebraically over $k$.

I'll request 120 digits then round to 100 for confidence but I didn't need to. Maple increases the working precision automatically, and this is confirmed by the fact that exactly the same answers are obtained if 100 digits are requested directly.

The documentation says that Levin's u-transform is used, with these references:

Fessler, T.; Ford, W.F.; and Smith, D.A. "HURRY: An acceleration algorithm for scalar sequences and series." ACM Trans. Math. Software, Vol. 9, (1983): 346-354.

Levin, D. "Development of non-linear transformations for improving convergence of sequences". Internat. J. Comput. Math, Vol. B3, (1973): 371-388.

evalf[120](Sum((2*k+3)/((k+2)*(k+1)/2)^2
    * sum( (2*l+3)/((k+2)*(k+1)/2 + (l+2)*(l+1)/2)^2,l=0..infinity),
       k=0..infinity)):
Z22 := evalf[100](%);
Z22 := 4.7059056441271748413982167408282763463686149628226287913\
         07611234496885020388225363260883986399848630 + 0. I

evalf[120](Sum((2*k+3)/((k+2)*(k+1)/2)^1
    * sum( (2*l+3)/((k+2)*(k+1)/2 + (l+2)*(l+1)/2)^3,l=0..infinity),
       k=0..infinity)):
Z13 := evalf[100](%);
Z13 := 1.6470471779364125793008916295858618268156925185886856043\
         46194382751557489805887318369558006800075685 + 0. I

evalf[100](2*Z22 + 4*Z13);
16.0000000000000000000000000000000000000000000000000000000000000\
   0000000000000000000000000000000000000 + 0. I