Finding the volume of a sphere using the Monte Carlo algorithm

I used the following code to find the volume of the sphere $x^2+y^2+z^2 \leq 1$ in the first octant:

McVolume[Num_] :=
Module[{hit, miss, index, x, y, z}, hit = 0; miss = 0;
For[index = 1, index <= Num, index = index + 1,
x = Random[Real, {0, 1}];
y = Random[Real, {0, 1}];
z = Random[Real, {0, 1}];
If[z <= Sqrt[1 - x^2 - y^2], hit = hit + 1, miss = miss + 1];];
Return[(hit/Num) 1];]

McVolume[100]


It turns out to be 11/25. Even thought the answer is correct but using this code i am getting an error:

LessEqual::nord: Invalid comparison with 0. +0.48202 I attempted. >>

I don't understand, why is this error popping up?

Any help would be appreciated.

The programming style you are using is not very fitting for Mathematica. Here's a better way (shorter, much faster):

n = 1000000; (* number of points to use *)
octantVolume = N[ Total@UnitStep[1 - Norm /@ RandomReal[1, {n, 3}]]/n ]


The reason why you get the error you mention is that for some x, y, the expression 1 - x^2 - y^2 is negative, thus its square root is not a real number and can't be use in inequalities. Use If[x^2+y^2+z^2 <= 1^2, ...] instead.

• +1 for providing an answer to the actual question. Commented Nov 10, 2013 at 3:23

Method of random number generation is also significant:

Default:

n = 10^6;
AbsoluteTiming[N@Mean@UnitStep[1. - Total[RandomReal[1, {3, n}]^2]] - π/6]


{0.197896, 0.000649224}

Niederreiter low-discrepancy sequence (see "methods" here):

SeedRandom[Method -> {"MKL", Method -> {"Niederreiter", "Dimension" -> 3}}];
AbsoluteTiming[N@Mean@UnitStep[1. - Total[RandomReal[1, {n, 3}]^2, {2}]] - π/6]


{0.139869, 0.0000132244}

Precision is >10 times better. It is equivalent to increasing the random sequence by 100-1000 times!

• To give an analogy with NIntegrate[], the second approach is equivalent to using "QuasiMonteCarlo", while the first one is akin to using "MonteCarlo". Commented Mar 24, 2016 at 5:48

When n is large it's much faster to operate on a 3 x n array than to process each of the n 3-vectors separately. This is one of the standard "tricks" to speed things up.

n = 10^6;  (*  Isn't that easier to read than 1000000 ?  *)

AbsoluteTiming @ N[ Total@UnitStep[ 1. - Norm/@RandomReal[1,{n,3}] ]/n ]
(* {4.555842, 0.524302} *)

AbsoluteTiming @ N[ Total@UnitStep[ 1. - Total[ RandomReal[1,{3,n}]^2 ] ]/n ]
(* {0.908908, 0.523554} *)

• +1. This question (and my answer there, as well as further links in there) seem related. Commented Nov 9, 2013 at 22:12
• +1. A bit more compact: N@Mean@UnitStep.... Commented Nov 9, 2013 at 22:25

I tried several ways to save time and make it more effective:

McVol1[num_] :=
Module[{hit, miss, index, x, y, z}, hit = 0; miss = 0;
For[index = 1, index <= num, ++index, {x, y, z} = RandomReal[{-1, 1}, 3];
If[x^2 + y^2 + z^2 <= 1, ++hit]];
hit/num]
Print["time and value...... :", Timing[N@McVol1[100000]]]


Result is

time and value...... :{0.452403, 0.52286}


The second is more effective

 McVol2[num_] :=
Plus @@ Table[Boole[Norm[RandomReal[{-1, 1}, 3]] <= 1], {num}]/num
Print["time and value...... :", Timing[N@McVol2[100000]]]