# How to get the most accurate volume of a special solid?

The solid is defined by implicit function below:

$$\frac{51}{100} (\cos x \cos y+\cos x \cos z+\cos y \cos z)+\cos x+\cos y+\cos z+1\le0$$ where $x,y,z\in[0,2\pi]$ I tried three different approaches with Mathematica to calculate the solid's volume:

1. Monte Carlo :

ClearAll["Global*"]; n = 1025; SeedRandom;

f[x_] := (Cos /@ x // Total) +  51/100*Dot[Cos /@ RotateRight @ x, Cos /@ x] + 1

Count[Flatten @ (f /@ RandomReal[{0, 2 Pi}, {n, 3}]), u_ /; u <= 0]/
n*(2.0 Pi)^3

2. immediate integral over implicit region:

NIntegrate[1, {x, y, z} ∈
ImplicitRegion[ 1 + Cos[x] + Cos[y] + Cos[z] +
51/100 (Cos[x] Cos[y] + Cos[x] Cos[z] + Cos[y] Cos[z]) <= 0 &&
0 <= x <= 2 Pi && 0 <= y <= 2 Pi && 0 <= z <= 2 Pi, {x, y, z}]]

3. Boole integrand over cube

NIntegrate[
Boole[1 + Cos[x] + Cos[y] + Cos[z] +
51/100 (Cos[x] Cos[y] + Cos[x] Cos[z] + Cos[y] Cos[z]) <= 0], {x,
0, 2 Pi}, {y, 0, 2 Pi}, {z, 0, 2 Pi}, WorkingPrecision -> 50]


Different approaches lead to different results.

Which one is more reliable and why since it seems there is no closed form result for this problem? Is it possible to evaluate the volume of this solid at arbitrary precision with Mathematica?

• Your code for 1. Monte Carlo doesn't return any value (it apparently just defines an f[x]). May 24, 2017 at 3:33
• Is your second one copied correctly? Trying to run it in Mathematica gives me an error. The {x,y,z}[Element]Implicit... part looks incorrectly entered. May 24, 2017 at 3:34
• @QuantumDot It does if you add in the missing newlines the OP left out with their copying. May 24, 2017 at 3:35
• Updated. thanks May 24, 2017 at 4:08
• Consider simplifying the problem by shifting and then treating just one octant: 8 NIntegrate[Boole[1 - Cos[x] - Cos[y] - Cos[z] + 51 (Cos[x] Cos[y] + Cos[x] Cos[z] + Cos[y] Cos[z])/100 <= 0], {x, 0, π}, {y, 0, π}, {z, 0, π}] May 24, 2017 at 8:08

Use symmetry, as J.M. suggested, split up the region....after some work....

8 (NIntegrate[1,
{x, 0, ArcCos[-(49/100)]},
{y, ArcCos[-((49 Cos[x])/(49 + 51 Cos[x]))], Pi},
{z,
Piecewise[{{ArcCos[(-100 - 100 Cos[x] - 100 Cos[y] -
51 Cos[x] Cos[y])/(
100 + 51 Cos[x] + 51 Cos[y])], -1 <= (-100 - 100 Cos[x] -
100 Cos[y] - 51 Cos[x] Cos[y])/(
100 + 51 Cos[x] + 51 Cos[y]) <= 1}}, 0],
Pi}
] +
NIntegrate[1,
{x, ArcCos[-(49/100)], Pi},
{y, 0, ArcCos[(-200 - 151 Cos[x])/(151 + 51 Cos[x])]},
{z,
ArcCos[(-100 - 100 Cos[x] - 100 Cos[y] - 51 Cos[x] Cos[y])/(
100 + 51 Cos[x] + 51 Cos[y])],
Pi}
] +
NIntegrate[1,
{x, ArcCos[-(49/100)], Pi},
{y, ArcCos[(-200 - 151 Cos[x])/(151 + 51 Cos[x])], Pi},
{z, 0, Pi}
])

(*  67.5935  *)

• Thank you. But can you also indicate the some work mentioned? :D May 25, 2017 at 0:21
• @LCFactorization Solve for the intersections in terms of Cos[z], then Cos[y] to determine the boundaries of the integration regions. I used Solve, but Reduce[inequality, {Cos[x], Cos[y], Cos[z]}] might work. It looks a little different, but there's probably more than one way. May 25, 2017 at 1:07
• Reduce[inequality && -1 < Cos[x] < 1 && -1 < Cos[y] < 1 && -1 < Cos[z] < 1, {Cos[x], Cos[y], Cos[z]}] give a better result for setting up the integrals. May 25, 2017 at 1:09
• I tried to increase the WorkingPrecision, and still got warning message like The global error of the strategy GlobalAdaptive has increased more than 2000 times. So this does not seem like a method assure numerical robustness. Am I right? May 25, 2017 at 1:51
• @LCFactorization Method -> {"GaussKronrodRule", "Points" -> 11}, WorkingPrecision -> 20 yields 67.5935383648647714596145728425065713295420., which agrees with the machine precision result to about 8 digits. After several minutes, Method -> {"GaussKronrodRule", "Points" -> 21}, WorkingPrecision -> 24 yields 67.5935383648622322295315138637653404890424.`, which agrees with the 20-digit result to 13 digits. The default method is much faster, but not as accurate. May 25, 2017 at 3:04