An integral with a fractional part in 3 dimensions

Question

The evaluation with Maple suggests the triple integral is around $1$, but Mathematica tells
it's $0.0958758$. When using the code

N[Integrate[FractionalPart[x/y] FractionalPart[y/z] FractionalPart[z/x], {x, 0, 1}, {y, 0, 1}, {z, 0, 1}]]

it returns:

"NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. >> NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 2000 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 0.09587582633950331and 0.004696948831748114 for the integral and error estimates. >>"

I'd like to know if Mathematica's approximation is correct, and if there is a possible closed form.

Answer 1

The analytic answer is $$ %\sum_{n=1}^\infty \frac{1}{2n(n+1)^2}\sum_{n=1}^\infty \frac{1+3n}{6n^2(n+1)^3}+\sum_{n=1}^\infty \frac{3n^2-1}{6n^2(n+1)^3}=\\ 1+\pi ^2\frac{2 \zeta (3)-9}{72} \approx 0.0958502 $$

Therefore, Mathematica is correct.

Proof

The 3D integral

NIntegrate[FractionalPart[x/y] FractionalPart[y/z] FractionalPart[z/x], 
   {x, 0, 1}, {y, 0, 1}, {z, 0, 1}]

0.0958741

depends only on ratios between x, y, and z. Therefore, it can be reduced to the 2D integral

NIntegrate[FractionalPart[x/y] y FractionalPart[1/x], {x, 0, 1}, {y, 0, 1}]

0.0958383

It has the following structure

ArrayPlot@Table[FractionalPart[x/y] y FractionalPart[1/x], 
    {x, 0.001, 1, 0.001}, {y, 0.001, 1, 0.001}]

We can split the integral to two parts

NIntegrate[FractionalPart[x/y] y FractionalPart[1/x], {x, 0, 1}, {y, 0, x}]
NIntegrate[FractionalPart[x/y] y FractionalPart[1/x], {x, 0, 1}, {y, x, 1}]

0.0176263

0.0782212

After substitution $y/x = ξ$ the first one becomes

NIntegrate[x FractionalPart[1/x], {x, 0, 1}] NIntegrate[x^2 FractionalPart[1/x], {x, 0, 1}]

0.0176233

The second one is equivalent to

NIntegrate[x (1 - x) FractionalPart[1/x], {x, 0, 1}]

0.0781823

These integrals can be written as a sum of

Integrate[x (1/x - n), {x, 1/(n + 1), 1/n}] // Simplify
Integrate[x^2 (1/x - n), {x, 1/(n + 1), 1/n}] // Simplify
Integrate[x (1 - x) (1/x - n), {x, 1/(n + 1), 1/n}] // Simplify

1/(2 n (1 + n)^2)
(1 + 3 n)/(6 n^2 (1 + n)^3)
(-1 + 3 n^2)/(6 n^2 (1 + n)^3)

Sum[1/(2 n (1 + n)^2), {n, 1, ∞}]
Sum[(1 + 3 n)/(6 n^2 (1 + n)^3), {n, 1, ∞}]
Sum[(3 n^2 - 1)/(6 n^2 (1 + n)^3), {n, 1, ∞}]

1/2 (2 - π^2/6)
1/6 (3 - 2 Zeta[3])
1/12 (6 - π^2 + 4 Zeta[3])

Finally, we obtain

1/2 (2 - π^2/6) 1/6 (3 - 2 Zeta[3]) + 1/12 (6 - π^2 + 4 Zeta[3]) // Simplify
N[%]

1 + 1/72 π^2 (-9 + 2 Zeta[3])
0.0958502

Answer 2

The Maple image you showed in your comment above was done in Maple 2D math (document mode). little hard to replicate it using classical interface (worksheet mode) which I use since one does not see the actual Maple commands that way.

The Maple commands you ave also might be using another package that was loaded before and not shown, such as Maple's Student. It is best that you post complete self contained Maple code to be able to compare.

But I have Maple 17, and not able to get the same result. Maple is not able to do this integral either.

restart:
f:=x->x-floor(x);
Int(f(x/y)*f(y/z)*f(z/x),[x=0..1,y=0..1,z=0..1]);
value(%);

Error, (in floor) too many levels of recursion

There is one thing I noticed though, is that Maple has a rule to integrate floor function, and Mathematica does not do it: Maple's definition is same as M:

floor - greatest integer less than or equal to a number

hence

restart:
int(floor(x),x);

    (*  floor(x)*x  *)

and in Mathematica, it remains unevaluated

Integrate[Floor[x], x]

But do you want to see something bizarre? I gave Integrate a rule integrate floor(x/y) to x*floor(x/y), which is also how it is defined in Maple:

 int(floor(x/y),x);
 (*  floor(x/y)*x *)

Unprotect Integrate and add the rule:

 (Unprotect[Integrate]; Integrate[Floor[x_/y_],x_] := x Floor[x/y];Protect[Integrate];)
 Integrate[x-Floor[x/y] , {x, 0, 1}]

and this is the result

Where did the 99999999999999999999/(100000000000000000000 y) come from? This should be 1

 N[%]

Very strange result. Here is Maple's result:

 int(x-floor(x/y),x=0..1);
 (* -floor(1/y)+1/2 *)

So same final result, but M exact answer has this 99999999999999999999/(100000000000000000000 y) term which only becomes 1.0 with using N on it. This might be related to the WorkingPrecision messages you are seeing.

The bottom line is: Maple does not give 1 as you showed in your screen image.