Question

Volume of a cap (formula given) using Monte Carlo integration (MATLAB).

The volume of a cap in terms of the sphere's radius, r, and depth, d, is given by Vcap = (3r – d) Suppose we did not know this expression or maybe we just want to check this expression. Assume r = 2 m and d = 1.5 m. Use Monte- Carlo integration to find the volume of such a cap. To start the prob- lem take 1,000 random numbers. When you are sure the code is run- ning, use 100,000. The number you compute should be pretty close to the one given by this formula. Figure 1: Definitions of h, d and r.

Answer 1

Integral formulation for volume of shoes opherical cap having radius to and height d combe awritten as - V.2 tl la- y²)dy. . ford Here 222, de 15 Hence ve 12 Ch - y²) dy e n Actual volume - 4-4%dy Veap a to (Bord) » Ixlos? (2x2–1-5) TRADE x 4.50 2 Txig 210. 602875205-

% % Monte Carlo simulation for finding volume of spherical cap clear all close all %function for integral f=@ (y) pi.*(4-y.^2); $Radius and height of the spherical cap r=2; d=1.5; % number of points generated N = 1000; al = r-d; a 2=0; b1 = r; b2=12; xx = al + (bl-al).*rand (N,1); yy = a 2 + (b2-a2). *rand (N,1); hits = sum(yy<=f(xx)); yl=linspace (al, bl); fl=f(yl); figure(1) plot(yl,fi, 'Linewidth', 2) Plotting data k=0; hold on for i=1:N if yy (i)<=f(xx(i)) k=k+1; xx(k)=xx(i); yy (k)=f(i); plot(xx(i),yy(i), 'r*') else plot(xx(i),yy(i), 'b*') end end xlim( [0.5 2]) ylim( [012]) title (sprintf('Monte Carlo plot for N=%d',N)) xlabel('x') ylabel('y') ar=(hits/N)*(1.5*12); fprintf('\n\tFor N=%d volume of spherical cap using Monte Carlo method is f.\n',N, ar ) %Radius and height of the spherical cap r=2; d=1.5; % number of points generated N = 10000; al = r-d; a 2=0; b1 = r; b2=12; xx = al + (bl-al).*rand (1,1);

yy = a2 + (b2-a2).*rand (N,1); hits = sum(yy<=f(xx)); yl=linspace (al,bl); fl=f(yl); figure(2) plot(yl,fi, 'Linewidth', 2) Plotting data k=0; hold on for i=1:N if yy(i)<=f(xx(i)) k=k+1; xx(k)=xx(i); yy(k)=f(i); plot(xx(i),yy(i), 'r*') else plot(xx(i),yy(i), 'b*') end end xlim( [0.5 2]) ylim( [012]) title (sprintf('Monte Carlo plot for N=%d',N)) xlabel('x') ylabel('y') ar= (hits/N)*(1.5*12); fprintf('\n\tFor N=%d volume of spherical cap using Monte Carlo method is %f.\n',N, ar ) For N=1000 volume of spherical cap using Monte Carlo method is 10.548000. For N=10000 volume of spherical cap using Monte Carlo method is 10.580400.

Monte Carlo plot for N=1000 0.5 1.5 Monte Carlo plot for N=10000 12 to 12 * SZ SOOS * 49 Lt 2. MANNE we * ERE WH SchW YOU *** SW . . ZRIN * Ses # 399 0.5 1.5
%%Monte Carlo simulation for finding volume of spherical cap
clear all
close all
%function for integral
f=@(y) pi.*(4-y.^2);
%Radius and height of the spherical cap
r=2; d=1.5;
% number of points generated
N = 1000;
a1 = r-d; a2=0;
b1 = r;   b2=12;

xx = a1 + (b1-a1).*rand(N,1);
yy = a2 + (b2-a2).*rand(N,1);

hits = sum(yy<=f(xx));
y1=linspace(a1,b1);
f1=f(y1);
figure(1)
plot(y1,f1,'Linewidth',2)
%Plotting data
k=0;
hold on
for i=1:N
   
     if yy(i)<=f(xx(i))
         k=k+1;
         xx(k)=xx(i);
         yy(k)=f(i);
         plot(xx(i),yy(i),'r*')
     else
         plot(xx(i),yy(i),'b*')
     end
   
end
xlim([0.5 2])
ylim([0 12])
title(sprintf('Monte Carlo plot for N=%d',N))
xlabel('x')
ylabel('y')
ar=(hits/N)*(1.5*12);
fprintf('\n\tFor N=%d volume of spherical cap using Monte Carlo method is %f.\n',N,ar)

%Radius and height of the spherical cap
r=2; d=1.5;
% number of points generated
N = 10000;
a1 = r-d; a2=0;
b1 = r;   b2=12;

xx = a1 + (b1-a1).*rand(N,1);
yy = a2 + (b2-a2).*rand(N,1);

hits = sum(yy<=f(xx));
y1=linspace(a1,b1);
f1=f(y1);
figure(2)
plot(y1,f1,'Linewidth',2)
%Plotting data
k=0;
hold on
for i=1:N
   
     if yy(i)<=f(xx(i))
         k=k+1;
         xx(k)=xx(i);
         yy(k)=f(i);
         plot(xx(i),yy(i),'r*')
     else
         plot(xx(i),yy(i),'b*')
     end
   
end
xlim([0.5 2])
ylim([0 12])
title(sprintf('Monte Carlo plot for N=%d',N))
xlabel('x')
ylabel('y')
ar=(hits/N)*(1.5*12);
fprintf('\n\tFor N=%d volume of spherical cap using Monte Carlo method is %f.\n',N,ar)
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%