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.

I used the following code to find the volume of the sphere $x^2+y^2+z^2 = 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.

share|improve this question

4 Answers 4

up vote 20 down vote accepted

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.

share|improve this answer
1  
+1 for providing an answer to the actual question. –  Jon Purdy Nov 10 '13 at 3:23

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} *)
share|improve this answer
    
+1. This question (and my answer there, as well as further links in there) seem related. –  Leonid Shifrin Nov 9 '13 at 22:12
1  
+1. A bit more compact: N@Mean@UnitStep.... –  ybeltukov Nov 9 '13 at 22:25

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!

share|improve this answer

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]]]
share|improve this answer

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.