An integral with a fractional part in 3 dimensions

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.

-
@Nasser this is what I received from someone (I don't have Maple) - i.stack.imgur.com/NzS8u.png – A_math_ninja Jan 1 '14 at 14:01
It's obvious that the integral is significantly smaller than $1$ therfore your Maple result is wrong if you've calculated it correctly. – Artes Jan 1 '14 at 14:17
@Artes that's true. I was very surprised with that result since, as you said, the integral is expected to be far smaller. – A_math_ninja Jan 1 '14 at 14:24
I believe there is a closed form, but the calculation is too tedious... Basic idea would be integral over every continuity regions respectively, and calculate the infinite summation. – Silvia Jan 2 '14 at 11:34
@Silvia I know, but I wanted Mathematica confirms that. – A_math_ninja Jan 2 '14 at 13:47

2 Answers

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

-
A very good job! Thanks! (+1) – A_math_ninja Jan 2 '14 at 16:28

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.

-
Thanks for your useful post! (+1) Indeed, some things out there are pretty strange ... – A_math_ninja Jan 1 '14 at 16:59
Nasser, I made a small spelling correction one minute after your last edit, figuring it would count as part of your edit. It didn't. Sorry, I wouldn't have used up an edit if I had known it would be charged. – DavidC Jan 1 '14 at 17:12
We also know the answer can't be 1 since the region of integration is contained in the unit cube. – Chip Hurst Jan 1 '14 at 23:25