Write a MATLAB script (M-File) to calculate pi number using MATLAB random number

Question

Write a MATLAB script (M-File) to calculate pi number using MATLAB random number generator rand(). Consider the circle and the square in the below figure. If N random points (with random X and Y) are generated inside the square, M random points will be inside the circle. The pi number can be estimated as: pi = 4 times M/N = 4 times AreaCircle/AreaSquare pi = 4 times AreaCircle/AreaSquare In your report, you should include the following: Problem statement: briefly explain the problem to be solved - Solution method: explain your solution method by writing mathematical equations (if required), physical relations (If required)... Flowchart Pseudocode MATLAB code Result Conclusion. Note You should email your MATLAB-file before 2.00PM Thursday 2/18.

Explanation / Answer

% this program calculates the value of "Pi" using sphere-cube method.
% sphere-cube method is analogous to circle-square-dart method
clear all;close all; %clear all residual variables and close all open figures.
tic % start clock.
PiValue=zeros(1,3); % 10 values of pi, each from one run is stored as elements in PiValue.
set(0,'DefaultFigureWindowStyle','docked') % dock the figures to matlab main window.
spherex=[];spherey=[];spherez=[];% initialise x y z coordinates of points inside sphere
annulusx=[];annulusy=[];annulusz=[]; % initialise x y z coordinates of points on and nearly inside the sphere surface
cubex=[];cubey=[];cubez=[];% initialise x y z coordinates of points inside cube and outside sphere
for count=1:size(PiValue,2) % do one trial after another till Pivalue vector is fully occupied.
    p3dm=zeros(16,32,64); % initialization of 3-D matrix storing pointers for x, y & z co-ordinates % for ease
    %of element access
    planes=size(p3dm,3);columns=size(p3dm,2);rows=size(p3dm,1);r=rows;c=columns; % extract variable values
    clear r c
    x=zeros(size(p3dm,1),size(p3dm,2),size(p3dm,3));y=x;z=y; % initialize 3-D matrices for co-ordinates
    x=rand(size(p3dm,1),size(p3dm,2),size(p3dm,3)); % generate random x-coordinates
    y=rand(size(p3dm,1),size(p3dm,2),size(p3dm,3)); % generate random y-coordinates
    z=rand(size(p3dm,1),size(p3dm,2),size(p3dm,3)); % generate random z-coordinates
    inside=0; % start counting, no of points inside the sphere of radius 0.5 units
    insideannulus=0; % start counting, no of points inside the sphere annulus
    insidecubeoutsidesphere=0; % start counting, no of points inside the cube and outside sphere
    for n=1:size(p3dm,1)*size(p3dm,2)*size(p3dm,3)
        distance=sqrt((x(n)-0.5)^2+(y(n)-0.5)^2+(z(n)-0.5)^2);
        if distance<=0.5 % distance criterion
            inside=inside+1; % increment by one of distance criterion satisfies
            if count==size(PiValue,2)
                spherex=[spherex x(n)];spherey=[spherey y(n)];spherez=[spherez z(n)];
            end
        end
        if distance<=0.5 & distance>0.49
            insideannulus=insideannulus+1; % increment by one of distance criterion satisfies
            if count==size(PiValue,2)
                annulusx=[annulusx x(n)];annulusy=[annulusy y(n)];annulusz=[annulusz z(n)];
            end
        end
        if distance>0.5
            insidecubeoutsidesphere=insidecubeoutsidesphere+1; % increment by one of distance criterion satisfies
            if count==size(PiValue,2)
                cubex=[cubex x(n)];cubey=[cubey y(n)];cubez=[cubez z(n)];
            end
        end
    end
   % pi_value=6*ratio of volume of sphere to volume of cube, which tends to
   % become equal to = 6*ratio of no of points inside sphere to total
   % points in cube at large random point density in cube.
    pi_value=6*(inside/(size(p3dm,1)*size(p3dm,2)*size(p3dm,3))); % calculate pi value
    PiValue(1,count)=pi_value; % fill the PiValue vector with the generated value of Pi; in the right place
    if count==1
        grid on;hold on % handle plot graphics
    end
    plot(count,PiValue(1,count),'ro','MarkerSize',5,'MarkerFaceColor','r','MarkerEdgeColor','k'); % plot points, each
    %representing value of obtained Pi vs. trial number
    pause(0.001) % pause to visualize the process
end
plot([1:size(PiValue,2)],PiValue(1,1:size(PiValue,2)),'k','LineWidth',1.5);grid on;axis tight % plot graph of different
%Pi values obtained and handle plot graphics
AverageValueOfPi=sum(PiValue)/size(PiValue,2) % calculate everage value
TrueValue=pi % true value of Pi
ErrorPercentage=sum((PiValue-pi)/pi*100)/size(PiValue,2) % percentage error in calculated average value
clear TrueValue columns rows planes x y z n distance % clear variables
count=1:size(PiValue,2); % create a vector of size same as PiValue - - > just for bar graph

figure;pause(0.25) % create new figure
bar(count,PiValue-pi);grid on;axis tight
xlabel('trial number');ylabel('obt. val - pi');title('deviation of obtained value from actual pi ');axis equal

pause(0.25);figure
plot3(spherex,spherey,spherez,'k.','MarkerSize',5);hold on;axis([0 1 0 1]);box on;grid on
xlabel('x-axis');ylabel('y-axis');zlabel('z-axis');title('Points inside sphere : of last run');axis equal

pause(0.25);figure
plot3(annulusx,annulusy,annulusz,'k.','MarkerSize',5);hold on;axis([0 1 0 1]);box on;grid on
xlabel('x-axis');ylabel('y-axis');zlabel('z-axis');title('Points on the sphere : of last run');axis equal

pause(0.25);figure
plot3(cubex,cubey,cubez,'k.','MarkerSize',5);hold on;axis([0 1 0 1]);box on;grid on
xlabel('x-axis');ylabel('y-axis');zlabel('z-axis');title('Points outside sphere : of last run');axis equal

clear count % clear variable
InAndTotal=[inside size(p3dm,1)*size(p3dm,2)*size(p3dm,3)]; % no of points inside the sphere an in the cube
clear inside % clear variable
toc % stop clock and show the 'elapsed time'

Related Questions

Navigate

• Browse (All)
• Browse W
• Subjects

---

---

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.