Finding the volume of a sphere using the Monte Carlo algorithm

Question

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.

Answer 1

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.

Answer 2

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!

Answer 3

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} *)

Answer 4

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