How to Calculate 100 times the volume of a 3D sphere using Monte Carlo simulation

14 ビュー (過去 30 日間)
Hi everybody, I have to calculate the volume of a 3D sphere using 50000 points, 100 times, and this is the code I have used to do it, but my teacher says the value of the vector with the results of all repetitions is wrong. I really can't find the issue, if somebody could help me I'd be very thankful. PD: There's text in spanish so I'm sorry for that, it's not necessary to understand the problem.
%% Borrado de datos
clear,clc;
%% Pedida del radio al usuario
%radio=input('El radio de la esfera a evaluar su volumen es:' ); % en el grader no se puede
r=2;
while isempty(r)==1 || isnumeric(r)~=1 || imag(r)~=0 || r-round(r)~=0 || r<=0
clc % Borrado de linea de comandos
disp('error, radio no valido')
r=input('Introduzca de nuevo el radio de la esfera a evaluar su volumen:' ); % en el grader no se puede
end
%% FUNCION A INTEGRAR PARA SACAR EL VOLUMEN, SOLO GRAFICA
x=linspace(0,r,50000);
y=4*pi*x.^2;
plot(x,y),grid()
%% FUNCION VOLUMEN REAL
v=@(r)(4/3)*pi*r.^3;
volumen=v(r);
%% Funcion area a integar
v1=@(x)4*pi*x.^2;
%% limites integracion
x1=v1(0);
x2=r;
y1=v1(x1);
y2=v1(x2);
%% CALCULO DE MONTECARLOS PARA HALLAR EL VOLUMEN
n=50000;
inside=0;
montecarlo=zeros(1,n);
for s=1:100
for k=1:n
xale=rand*(x2-x1)+x1;
yale=rand*(y2-y1)+y1;
if yale <= 4*pi*xale.^2 % FUNCION AREA SOBRE LA QUE SE INTEGRA
inside=inside+1;
end
montecarlo(k)=(dentro/k)*y2*x2;
end
montecarlo1(s)=montecarlo(k);
end
%% Calculo de la media de los 100 montecarlos
volumeM=length(montecarlo1); % valor del tamano vector
mediaM=montecarlo1(s)/volumeM; % valor 100 del vector entre 100
error=100*abs((mediaM-volumen)/volumen);
%% Resultados
disp('El volumen por la formula es: '),disp(volumen)
disp('El volumen por el metodo de Montecarlo 100 veces es:' ),disp(mediaM)
disp('El error del metodo es:'),disp(error)

採用された回答

John D'Errico
John D'Errico 2019 年 6 月 27 日
編集済み: John D'Errico 2019 年 6 月 27 日
First, you should already know the volume of a sphere. Is it producing reasonable numbers? If not, then why did you even bother to hand it in?
As for trying to read what you did, I won't even try, because what you did seems to make little sense. I can't even try to guess what you were thinking to write that code. For example, you are asked to comute the volume of a sphere. So why do you have formulas to compute the area?
But, you did do some work here. Sigh. And I can't fix what I can't even follow sufficiently to figure out what you are trying to do.
% The radius of the sphere. Assume the center lies at the origin.
r = 2;
% the volume of a cube of side length 2*r is 8*r^3.
cubevol = 8*r^3;
np = 50000;
nsim = 100;
predvol = zeros(1,nsim);
for isim = 1:nsim
pts = (rand(np,3)*2 - 1)*r;
% what fraction of the points from the cube of side length 2*r,
% happen to fall inside the sphere?
R = sqrt(sum(pts.^2,2));
insidecount = sum(R <= r);
% the estimated volume of the sphere is just the fraction of
% cube points that happen to lie inside the sphere, times the
% volume of the cube.
spherevol(isim) = cubevol*insidecount/np;
end
Now, does this work? Compare the set of predicted sphere volumnes, versus the actual volume.
Min 33.26
1.0% 33.26
5.0% 33.3
10.0% 33.32
25.0% 33.4
50.0% 33.52
75.0% 33.61
90.0% 33.67
95.0% 33.74
99.0% 33.83
Max 33.86
Mean 33.51
Std 0.1385
Rms 33.51
Range 0.6042
>> 4/3*pi*r^3
ans =
33.5103216382911
And that is the most important thing you should have done to test your reults, BEFORE you handed it in.
  2 件のコメント
Martin Avagyan
Martin Avagyan 2019 年 6 月 27 日
Well, I was trying to integrate the area so I could get back the volume of the sphere and then do the simulation, but now I see it was not a good decision. Thank you very much
John D'Errico
John D'Errico 2019 年 6 月 27 日
編集済み: John D'Errico 2019 年 6 月 27 日
I was wondering if you were trying to integrate the area. Hmm. Perhaps your mistake was in not doing the integral properly.
Lets see how that might have worked, first as an analytical integral.
syms r
>> int(4*pi*r^2,[0,2])
ans =
(32*pi)/3
>> vpa(ans)
ans =
33.510321638291127876934862754981
Ok, so that makes sense. Now, how might it work as a Monte Carlo? In the solution I wrote in my answer, I used the fraction of points that lie inside the sphere to compute the volume. But we can do the integral that you attempted as:
rmax = 2;
A = @(R) 4*pi*R.^2;
np = 5000000;
rsamp = rmax*rand(np,1);
mean(A(rsamp))*(rmax - 0)
ans =
33.5025574971549
So this Monte Carlo is a subtly different variety of Monte Carlo, although they can be shown to be related.
In this case, I simply evaluated the function a lot of times, then computed the mean value of that function, multiplying by the interval width.
Either way is valid.

サインインしてコメントする。

その他の回答 (1 件)

Martin Avagyan
Martin Avagyan 2019 年 6 月 27 日
That was very useful, thank you for your time! :)

カテゴリ

Help Center および File ExchangeDiscrete Math についてさらに検索

製品


リリース

R2017b

Community Treasure Hunt

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

Start Hunting!

Translated by