Your simulation's convergence to an incorrect result is because your comparison accepts any negative y, thus inflating the result. You are in essence computing the area of the union of the semi-unit circle that lies above the x-axis and the rectangle with top-left corner at {-1,0} and lower-right corner at {1,-1}, which is 2 + π/2 = 3.5708
.
This is the minimal repair to your code.
McArea[Num_] :=
Module[{hit, miss, index, x, y},
hit = 0; miss = 0;
For[index = 1, index <= Num, index = index + 1,
x = Random[Real, {-1, 1}];
y = Random[Real, {-1, 1}];
If[Abs[y] <= Sqrt[1 - x^2], hit = hit + 1, miss = miss + 1];];
Return[(hit/Num) 4];]
Timing[N @ McArea[100000]]
{0.590944, 3.13596}
Now let's explore some improvements. The first still uses a For-loop, but with a slight better comparison expression and some improved Mathematica practice.
mcArea2[num_] :=
Module[{hit, index, x, y},
For[hit = 0; iindex = 1, index <= num, ++index,
{x, y} = RandomReal[{-1, 1}, 2];
If[x^2 + y^2 <= 1, ++hit]];
4. hit/num]
Timing[mcArea2[100000]]
{0.571963, 3.14488}
Note the use of lowercase initial letters on all user defined identifiers. Using uppercase initial letters is bad practice because it can get you into conflicts with Mathematica's predefined, built-in variables. Also, note the use ++
and the removal of the unnecessary Return
.
Do
is better than For
in this case since we don't need any of fine control over iteration that For
provides.
mcArea3[num_] :=
Module[{hit, x, y},
hit = 0;
Do[{x, y} = RandomReal[{-1, 1}, 2]; If[x^2 + y^2 <= 1, ++hit], {num}];
4. hit/num]
Timing[mcArea3[100000]]
{0.476420, 3.14056}
Edit
People have started to submit functional programming answers where the connection to the original question is no longer obvious, so I will add mine.
mcArea4[num_] :=
4. Plus @@ Table[Boole[Norm[RandomReal[{-1, 1}, 2]] <= 1], {num}]/num
Timing[mcArea4[100000]]
{0.032864, 3.14664}
If[Abs[y] <= ...
or two Ifs, one for positive y and one for negative y. $\endgroup$ – Peltio Nov 15 '13 at 11:21RandomReal
. Also take a look at loops in Mathematica what may help you in improving the performance of the code. p.s. why are you countingmiss
? $\endgroup$ – Kuba♦ Nov 15 '13 at 11:27RegionPlot[y <= Sqrt[1 - x^2], {x, -1, 1}, {y, -1, 1}]
$\endgroup$ – Simon Woods Nov 15 '13 at 19:12