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

14 views (last 30 days)
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)

Accepted Answer

John D'Errico
John D'Errico on 27 Jun 2019
Edited: John D'Errico on 27 Jun 2019
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 Comments
Martin Avagyan
Martin Avagyan on 27 Jun 2019
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 on 27 Jun 2019
Edited: John D'Errico on 27 Jun 2019
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.

Sign in to comment.

More Answers (1)

Martin Avagyan
Martin Avagyan on 27 Jun 2019
That was very useful, thank you for your time! :)

Categories

Find more on Discrete Math in Help Center and File Exchange

Products


Release

R2017b

Community Treasure Hunt

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

Start Hunting!