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.

Szabolcs · Accepted Answer · 2013-11-09T18:04:44.913

25

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.

edited Nov 09 '13 at 18:04

answered Nov 09 '13 at 17:58

Szabolcs

234,956
30
623
1,263

1

+1 for providing an answer to the actual question. – Jon Purdy Nov 10 '13 at 03:23

ybeltukov · Answer 2 · 2013-11-10T13:19:08.883

18

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!

edited Nov 10 '13 at 13:19

answered Nov 09 '13 at 22:17

ybeltukov

43,673
5
108
212

To give an analogy with NIntegrate[], the second approach is equivalent to using "QuasiMonteCarlo", while the first one is akin to using "MonteCarlo". – J. M.'s missing motivation Mar 24 '16 at 05:48

score 17 · Answer 3 · answered Nov 09 '13 at 20:27

17

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

answered Nov 09 '13 at 20:27

Ray Koopman

3,306
14
13

+1. This question (and my answer there, as well as further links in there) seem related. – Leonid Shifrin Nov 09 '13 at 22:12
1

+1. A bit more compact: N@Mean@UnitStep.... – ybeltukov Nov 09 '13 at 22:25

score 0 · Answer 4 · edited Aug 04 '14 at 06:36

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

Finding the volume of a sphere using the Monte Carlo algorithm

4 Answers4

Linked

Related