Monte Carlo integration (hit or miss) to find the area of a circle of radius R

5 次查看（过去 30 天）
Left Terry
Left Terry 2025-1-31
评论： Walter Roberson 2025-1-31
Ran in:
Hi everyone. Trying to solve an old exam topic regarding Monte Carlo integration, I wrote the following code for which I based it on a code from a professor in the C language. For R = 1 it worked so I assume I did it correctly. Then I tried R > 1 and the results were wrong compared to the formula πR^2 so I tried a few modifications by trial and error method. I fixed the problem and it seems to works correctly but there is a modification that I don't understand. If you could help me understand this it would be great. Thank you in advance.
clear, clc, format short
%---------- Computation for R = 1 ----------
R = 1;
N = 1e+7;
insideCircle = 0;
for i = 1:N
x = rand(1);
y = rand(1);
if (x^2 + y^2) < R
insideCircle = insideCircle + 1;
end
end
Area_MC = 4*insideCircle/N;
Area = pi*R^2;
error = 100*abs(Area - Area_MC)/Area;
fprintf('\n\tExact area\t\t%f\n\tMonte Carlo area\t%f\n\tPercentage error\t%.5f%%\n',Area,Area_MC,error);
Exact area 3.141593 Monte Carlo area 3.142050 Percentage error 0.01456%
%---------- Computation for R > 1 (e.g. R = 4 but you can try any real value) ----------
R = 4; %input('Enter radius R: R = ','s'); R = str2double(R);
% while R <= 0 || isreal(R) == 0 || isnan(R) == 1 || isnumeric(R) == 0
% disp('Invalid input')
% R = input('Enter radius R: R = ','s');
% R = str2double(R);
% end
N = 1e+7;
insideCircle = 0;
for i = 1:N
x = rand(1)*R;
y = rand(1)*R;
if (x^2 + y^2) < R
insideCircle = insideCircle + 1;
end
end
Area_MC = 4*R^3*insideCircle/N; % Modification needed that I don't understand (I had to multiply with R^3)
Area = pi*R^2;
error = 100*abs(Area - Area_MC)/Area;
fprintf('\n\tExact area\t\t%f\n\tMonte Carlo area\t%f\n\tPercentage error\t%.5f%%\n',Area,Area_MC,error);
Exact area 50.265482 Monte Carlo area 50.285594 Percentage error 0.04001%

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2025-1-31
编辑：Torsten 2025-1-31
You compare the area of the square with side length R and corner points (0,0), (R,0), (R,R) and (0,R) with the area of the quarter circle of radius R. This square has area R^2, and x and y are generated to cover this square uniformly.
N = 1e+7;
insideCircle = 0;
for i = 1:N
x = R*rand();
y = R*rand();
if (x^2 + y^2) < R^2
insideCircle = insideCircle + 1;
end
end
Area_MC = 4*R^2*insideCircle/N; % Modification needed that I don't understand (I had to multiply with R^3)
Area = pi*R^2;
error = 100*abs(Area - Area_MC)/Area;
fprintf('\n\tExact area\t\t%f\n\tMonte Carlo area\t%f\n\tPercentage error\t%.5f%%\n',Area,Area_MC,error);
  7 个评论
显示 5更早的评论
Torsten
Torsten 2025-1-31
编辑：Torsten 2025-1-31
I didn't realize at first that you are working in a quarter circle instead of a full circle.
For this approach, the only thing that was wrong in your original code is that you used "if (x^2 + y^2) < R" instead of "if (x^2 + y^2) < R^2".
Of course, "Area_MC = 4*R^3*insideCircle/N;" has to be reset to "Area_MC = 4*R^2*insideCircle/N;"
Walter Roberson
Walter Roberson 2025-1-31
By the way I don't see any difference between rand() and rand(1) in the plots
X = rand returns a random scalar drawn from the uniform distribution in the interval (0,1).
X = rand(n) returns an n-by-n matrix of uniformly distributed random numbers.
... so when n = 1, rand(1) returns a 1 x 1 matrix of uniformly distributed random numbers, which is the same as returning a random scalar. rand() and rand(1) mean the same thing.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Half-Normal Distribution 的更多信息

标签

产品


版本

R2016a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by