MATLAB: Monte carlo: intersection volume of cylinders

intersectionmonte carloplotvolume

I have to calculate the intersection volume of three cylinder(radius=3) each of which lie on x,y,z axis perpendicular to one another. I'm using a 5x5x5 cube

I can't use any of the more advanced matlab functions as we haven't been taught it yet.

I tried manipulating code for the area of a circle in a 2D square to fit my current problem.

Here's what I got:

N= 10000; % number of points generated
a = -5;
b = 5;
r=3;
hits=0;
x = a + (b-a).*rand(N,1);
y = a + (b-a).*rand(N,1);
z = a + (b-a).*rand(N,1);
radiixy = sqrt(z.^2+y.^2);
radiixz = sqrt(x.^2+z.^2);
radiizy = sqrt(x.^2+y.^2);
i = radiixy  && radiixz && radiizy <= r;
%count the hits
for j = 1:N
  hits=hits+i(j);
end
misses = N-hits;
disp('hits')
disp(hits)
disp('misses')
disp(misses)
plot(x(i),y(i), z(i),'.g');
hold;
plot(x(~i),y(~i),z(~i),'.r');
xlabel('x');
ylabel('y');
zlabel('z');
title('Intersection Volume');

ERROR: Operands to the and && operators must be convertible to logical scalar values

Best Answer

It looks like you packed all of your errors into one line of code. Everything else seems reasonable.

i = radiixy  && radiixz && radiizy <= r;

First of all, you use & as the operator here, not &&. Use && ONLY for scalar tests, essentially in an if statement. This is because the && operator is a short-circuited one. For example, in the logical statement:

(0 == 1) && (a == b)

MATLAB is smart enough to not bother to evaluate the second clause, because it sees that the first part is always false, therefore the statement MUST be false. This is why the && and operators were introduced.

Next, you may think that this tests if all of the radii are less than r, it does nothing of the sort, even if you used & instead of &&.

Do this instead:

i = (radiixy <= r) & (radiixz <= r) & (radiizy <= r);

Just because you might think in some short-hand notation, this does not mean you should your computer to understand what you mean.

Next, why have a loop to count the hits? Is sum really that advanced? How about this:

hits = ones(1,N)*i;

Related Solutions

MATLAB: Basic monte carlo question: area of circle

You are mixing things up. Your code is part vectorized, part scalar.

You have a loop over N, but then at each iteration you generate N points (so you're generating N^2 points overall). The variables x, y, and radii are vectors.

Note that

 [6 7 4 2 8 5] <=5

gives

 [0 0 1 1 0 1]

Therefore condition radii<=5 has a very low probability to be met ( all of the N radii should be less than 5).

As a result hits is (almost) never accumulated.

Your code is structured for a scalar radius. It would work if you substitute the lines

 x = a + (b-a).*rand(N,1);
 y = a + (b-a).*rand(N,1);

with

 x = a + (b-a).*rand;
 y = a + (b-a).*rand;

In this case hits will be a number between 0 and N, and hits/N is going to be proportional to the ratio of the areas of the circle and square, pi*(5/14)^2 (note that you can use simple power operator ^ instead of elementwise .^)

You can avoid the loop altogether using the vectorized variables.

 N= 1000; % number of points generated
 a = -7;
 b = 7;
 x = a + (b-a).*rand(N,1);
 y = a + (b-a).*rand(N,1);
 radii = sqrt(x.^2+y.^2);
 hits = sum(radii<=5);

MATLAB: [x,y,z] = sphere, except uniform or random distributed points

Here is a way to generate random points on a sphere that is statistically uniform.

   r = randn(3,n); % Use a large n
   r = bsxfun(@rdivide,r,sqrt(sum(r.^2,1)));
   x = r(1,:);
   y = r(2,:);
   z = r(3,:);

If you want a radius other than one, multiply the second line by the desired radius, R.

Note: If you use 'rand' rather than 'randn' in the above, the distribution will not be statistically uniform.

Best Answer

Related Solutions

MATLAB: Basic monte carlo question: area of circle

MATLAB: [x,y,z] = sphere, except uniform or random distributed points

Related Question