Estimating value of pi using Monte Carlo Method in Scilab

 Approximate Pi in Scilab 

Monte Carlo Method is a spectrum of computational algorithms which use randomness to solve problems which are deterministic in nature.

Here, we are estimating the value of  pi in Scilab by using the above stated method.

Suppose we have a square of side 'r' with quarter circle of radius 'r' inscribed in it.

Consider a game of darts with a dartboard just like the above one. You have been given N darts and of those N darts you have managed to hit n darts inside(including the edge of circle). Now, what is the probability that a dart hitting inside circle would be 

This probability is also equal to ratio of area of circle to area of square i.e,

The value of  N should be very large for better accuracy of  pi.

Now you might be wondering as to how probability is equal to ratio of areas and to answer that one might consider a similar situation of a square of side 2 unit with 4 squares each of side unit in it. 

Probability of choosing any unit square randomly is

 If we take a look at the ratio of areas i.e.,

from this we can infer that probability is equal to ratio of area i.e.,

Algorithm 

For Estimating Pi in Scilab

1.  input N(No. of sampling points), define r(radius of circle or side of square),ptcir(No. of point in circle),ptsq(No. of point in square)
2. initialize array x, y having dimension(1,N) (for storing random numbers from 0 to r) using function rand() 
3. initialize loop from 1 to N,inside loop 
• define 
  
• condition if  
• increment ptcir
• define array val i.e.,
4.  mean(val)

For Plotting in Scilab

1. define array theta  having equally spaced value from 0 to 2pi.
2. plot quadrant of circle using polarplot() function.
3.  initialize loop from 1 to N,inside loop 
• define 
  
• condition if   use plot() function to plot point giving a color 

• else plot point using some other color.

 Scilab Program

 //To find value of pi using Monte Carlo simulations

clc;clf;clear;
N=900//No. of pts to be generated
x=rand(1,N)//array of random no. in x-axis
y=rand(1,N)//array of random no. in y-axis
r=1//Radius of circle 
ptcir=0//points on circle
ptsq=N//points on square
for i=1:N
    d=sqrt(x(i)**2+y(i)**2)
    if d<=1
        ptcir=ptcir+1
    end   
    val(i)=4*(ptcir/ptsq)
end
avgval=mean(val)
theta = 0:.01:%pi/2;
polarplot(theta,r)
for i=1:N
    d=sqrt(x(i)**2+y(i)**2)
    xtitle("Approximate value of pi: "+string(val(i))+"    
No. of pts to be generated,N: "+string(i),"X-axis","Y-axis")
    if d<=1
        plot(x(i),y(i),'ro')
    else
        plot(x(i),y(i),'go')
    end   
end

Scilab Output

Given below give a little understanding to how by increase number of sampling would give out better approximate value of pi.

Bibliography

1. GeeksforGeeks
2. 101computing  
3. Wikipedia 
4. Scilab Help