Volume Monte Carlo Method

My textbook handed us this code for a sphere with radius 1. I have to change this to estimate the mass of a sphere of radius 1 with mass density x^2*y^2*z^2. If anyone could help me better understand what the code they gave me means and how I manipulate it to do this, it would be greatly appreciated.

function vol = volMonteCarlo(region,x0,x1,y0,y1,z0,z1,n)

X = rand(n,3)

X(:,1)= x0+(x1-x0)*X(:,1);

X(:,2)= y0+(y1-y0)*X(:,2);

X(:,3)= z0+(z1-z0)*X(:,3);

vol = sum(region(X))*(x1-x0)*(y1-y0)*(z1-z0)/n;

end

sphere = @(X) heaviside(ones(size(X,1),1) ...

- X(:,1).*X(:,2).*X(:,3));

volMonteCarlo(sphere,-1,1,-1,1,-1,1,100000)

[Cassandra Meyer]

Oh so sorry, thats my fault. They were given as this, I forgot I changed it trying to figure it out.

-X(:,1).^2 - X(:,2).^2 - X(:,3).^2);

Answer 1

[Torsten]

I'm surprised that the sphere is characterized as

sphere = @(X) heaviside(ones(size(X,1),1) - X(:,1).*X(:,2).*X(:,3));

For me it seems that this is the complete cube [-1,1] x [-1,1] x [-1,1] from which the samples are chosen.

In my opinion,

sphere = @(X) (X(:,1).^2 + X(:,2).^2 + X(:,3).^2 <= 1)

so that you get your integral by

function_over_sphere = @(X) (X(:,1).^2 + X(:,2).^2 + X(:,3).^2 <= 1).*(X(:,1).^2.*X(:,2).^2.*X(:,3).^2);
volMonteCarlo(function_over_sphere,-1,1,-1,1,-1,1,100000)

To make one more test, you could try to calculate the volume of the sphere with radius 1 by integrating the constant function 1 over the sphere:

sphere = @(X)(X(:,1).^2 + X(:,2).^2 + X(:,3).^2 <= 1)*1

As you know, the result should approximately be V = 4/3*pi.

The example under

https://en.wikipedia.org/wiki/Monte_Carlo_integration

should demonstrate how Monte Carlo Integration works and help to understand the code from above.