Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
    
@Nasser this is what I received from someone (I don't have Maple) - i.stack.imgur.com/NzS8u.png –  Chris's sis Jan 1 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 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. –  Chris's sis Jan 1 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 at 11:34
    
@Silvia I know, but I wanted Mathematica confirms that. –  Chris's sis Jan 2 at 13:47

2 Answers 2

up vote 6 down vote accepted

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}]

enter image description here

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
share|improve this answer
1  
A very good job! Thanks! (+1) –  Chris's sis Jan 2 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

Mathematica graphics

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]

Mathematica graphics

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

Mathematica graphics

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

 N[%]

Mathematica graphics

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.

share|improve this answer
    
Thanks for your useful post! (+1) Indeed, some things out there are pretty strange ... –  Chris's sis Jan 1 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. –  David Carraher Jan 1 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 at 23:25

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.