# Evaluation of a triple sum does not finish in reasonable time

I'm trying to compute the following triple sum, but no result is produced within a reasonable amount of time. What to do?

Sum[1/( i j k (i + j + k + 1)), {i, 1, Infinity}, {j, 1, Infinity}, {k, 1, Infinity}]


I use Mathematica 8.0.

-
This works... ParallelSum[Sum[Sum[1/(i j k (i + j + k + 1)), {i, 1, Infinity}] /. j -> j1, {k, 1, Infinity}], {j1, 1, 50}]  but note that it stops at j=50 –  chris May 8 '13 at 9:39
FYI its value is 5.35132... does not seem to converge quickly though... –  chris May 8 '13 at 10:30
note that the value I gave you is for jmax=50 not infinity! –  chris May 8 '13 at 10:38
@chris You can get a closed form for one summation, it might help convergence. –  b.gatessucks May 8 '13 at 12:39
@chris and are the two of you siblings ? –  b.gatessucks May 8 '13 at 12:40

@whuber gave me another idea : first change the summand to an exponential and then do the sums using the symmetry.

Integrate[Exp[-t i j k (i + j + k + 1)], {t, 0, Infinity}, Assumptions -> {i > 0, j > 0, k > 0}]
(* 1/(i j k (1 + i + j + k)) *)


so we need to do the sums of the exponential and then integrate at the end. Now we can make a change of variables t i j k = y which will simplify the exponent and bring a factor 1/(i j k) from the Jacobian. This makes it easier to see that

$$\sum_{i,j,k=1}^{\infty} \frac{1}{i j k} \exp(-y\ (i+j+k+1)) = \exp(-y) \left( \sum_{i=1}^{\infty} \frac{1}{i}\exp(-y\ i)\right)^3$$

Simplify[Sum[Exp[-y i]/i, {i, 1, Infinity}, Assumptions -> {y > 0}],
Assumptions -> {y > 0}]
(* -Log[1 - E^-y] *)


and we are left with the integration :

Integrate[Exp[-y] (-Log[1 - E^-y])^3, {y, 0, Infinity}]
(* 6 *)


It's also straightforward to generalize to the sum over n integers and the result is :

Integrate[Exp[-y] (-Log[1 - E^-y])^n, {y, 0, Infinity},
Assumptions -> {n \[Element] Integers, n >= 1}]
(* Gamma[1 + n] *)

-
Aside from justifying switching the order of integration and summation and a typo (missing $i$ in second summation), it looks good. Numerically I got an upper bound of about 6.08 on the sum. (+1) –  Michael E2 May 8 '13 at 19:26
@MichaelE2 Fixed, thanks a lot. –  b.gatessucks May 8 '13 at 19:28
Very nice idea! @Michael E2 has clearly identified my mistake, which now squares our two results. (I promised only to find a lower bound for the sum, so a numerical value of the integral around $3.5$ is just fine). I apologize for making you work so hard :-) and will make up for it by transferring all the points from my answer to yours, which deserves it. (There's a day's wait before I can do that, though...) –  whuber May 8 '13 at 19:37
@whuber You had the key insight and deserve most of the credit. –  b.gatessucks May 8 '13 at 19:44
btw, congrats on your 10k! –  The Toad May 11 '13 at 22:20

Here is another variant along the same lines as @whuber proposed. Consider function $$f(x,y,z,t)=\sum _{i=1}^{\infty } \sum _{j=1}^{\infty } \sum _{k=1}^{\infty } x^{i-1} y^{j-1} z^{k-1} t^{i+j+k}.$$ Mma gives

f[x_, y_, z_, t_] =
Sum[x^(i - 1) y^(j - 1) z^(k - 1) t^(i + j + k), {i,
1, \[Infinity]}, {j, 1, \[Infinity]}, {k, 1, \[Infinity]}]


$$-\frac{t^3}{(t x-1) (t y-1) (t z-1)}.$$ Then the required sum is equal to $$\int _0^1\int _0^1\int _0^1\int _0^1f(x,y,z,t)\,dzdydxdt,$$

Integrate[f[x, y, z, t], {t, 0, 1}, {x, 0, 1}, {y, 0, 1}, {z, 0, 1}]


which gives 6 as in the answer of @b.gatessucks.

-

If I take the challenge to be to get Mathematica to do the sum in the form given (as distinct from applying some mathematics insightfully), then here is a way to get the answer. The only trick I use is to break the triple sum into a double and single sum.

s2 = Sum[1/(i j (i + j + k + 1)), {i, 1, Infinity}, {j, 1, Infinity}]; // Timing
(s2 = FullSimplify[s2]) // Timing
Sum[s2/k, {k, 1, Infinity}] // Timing

{93.696694, Null}
{0.234598, (\[Pi]^2 + 6 HarmonicNumber[1 + k]^2 - 6 PolyGamma[1, 2 + k])/(6 + 6 k)}
{12.027743, 6}


The intermediate step of simplifying s2 seems to be necessary.

Interestingly Mathematica caches auxiliary expressions in computing the sum. The second time around is faster. (Almost all the time is spent in FullSimplify)

(s2 = Sum[1/(i j (i + j + k + 1)), {i, 1, Infinity}, {j, 1, Infinity}];
s2 = FullSimplify[s2];
Sum[s2/k, {k, 1, Infinity}]) // Timing

{0.228539, 6}

-

Not an answer, but it might help. You can get one summation in closed form :

sum1 = Sum[1/(i j k (i + j + k + 1)), {i, 1, Infinity}]
(* (EulerGamma + PolyGamma[0, 2 + j + k])/(j k (1 + j + k)) *)


You can also do the easy bit :

sum2 = Sum[Sum[EulerGamma/(j k (1 + j + k)), {j, 1, Infinity}], {k, 1, Infinity}]
(* 2 EulerGamma *)


so you're left with a double sum

sum = 2 EulerGamma +
Sum[PolyGamma[0, 2 + j + k]/(j k (1 + j + k)), {j, 1, Infinity}, {k, 1, Infinity}]


You can make some further simplification by using the symmetry of the summand or proceed numerically.

-