Introduction

The number \(\pi\) is a mathematical constant, approximately equal to 3.14159. It is defined in Euclidean geometry as the ratio of a circle’s circumference to its diameter.
Being an irrational number, \(pi\) cannot be expressed as a common fraction, although fractions such as \(\frac{22}{7}\) are commonly used to approximate it. Equivalently, its decimal representation never ends and never settles into a permanently repeating pattern. Its decimal (or other base) digits appear to be randomly distributed, and are conjectured to satisfy a specific kind of statistical randomness.
Here we will approximate the \(\pi\) using simulation.

Motivation

Let consider a circle of radios 1 which is inscribed under a square of side 2. Then the area of the circle is given by, \[ Area = \pi.(1)^2 = \pi \] And the area of the square is 4.

Lets throw \(N\) dirts into the square. Out these \(N\) dirts some of the dirts will be hit the circle and andd some will not. Let \(n\) be the numbers of dirtsout of the \(N\) dirts, which hit the circle. Therefore the proportion of of the dirts which hit the circle is an approximation of \(\pi /4\). \[i.e. \frac{n}{N}=\frac{\pi}{4}\] Hence, an approximation of \(\pi\) is given by, \[\hat{\pi} = 4 \times \frac {n}{N}\]. For large \(N\) like 10,000 ; 20,000 ; 50,000 etc this method gives satisfactory result.

Algorithm

Since, this method is appropriate with large \(N\), it is impossible to workout the simulation practically by hand, we need the help of programming here. To execute the above mentioned problem we may follow the following algorithm,
-Step 1:Draw \(N\) random samples from Uniform (-1,1) distribution for each of the variable \(x\) and \(y\).
-Step 2: Find the number of (\(x,y\)) which lie inside the circle \(x^2+y^2=1\). Let \(n\) be the number.
-Step 3: Then the approximation of \(\pi\) is \(4*(n/N)\).

Program

To execute the simulation we will use R here.

Approach 1

apx_pi<-function(N,draw=FALSE)
{
  x<-runif(N, min = -1, max = 1)
  y<-runif(N, min = -1, max = 1)
  
  m<- rep(0,N)
  for (i in 1:N) {
   if(x[i]^2+y[i]^2<=1)
     m[i]<-1
   else
     m[i]<-0
     }
  
sum(m)
p<-4* (sum(m)/N)

x1<-rep(0,N)
y1<-rep(0,N)
for (i in 1:N)
  {
  if(x[i]^2+y[i]^2<=1)
    {
    x1[i]<-x[i]
    y1[i]<-y[i]
    }
  else
    {
    x1[i]<-0
    y1[i]<-0
    }
}
library(ggplot2)
print( ggplot(aes(x,y),data=NULL)+
    geom_point(col="blue")+
    geom_point(aes(x1,y1),col="yellow"))
return(p)

}

apx_pi(100)

## [1] 3.12

apx_pi(20000)

## [1] 3.146

apx_pi(50000)

## [1] 3.13936

Approach 2

apx_pi2<-function(N)
{
  x<-runif(N, min = -1, max = 1)
  y<-runif(N, min = -1, max = 1)
  
  m<- rep(0,N)
  for (i in 1:N) {
   if(x[i]^2+y[i]^2<=1)
     m[i]<-1
   else
     m[i]<-0
     }
  
sum(m)
p<-4* (sum(m)/N)
x1<-rep(0,N)
y1<-x1
for (i in 1:N)
{
  if(x[i]^2+y[i]^2<=1)
  {
    x1[i]<- x[i]
    y1[i]<- y[i]
  }
}
plot(x1,y1,pch=".",col="blue")    ######giving the plot
abline(h=0,,v=0, col='yellow', lwd=2)
abline(h=c(1,-1), col=c('red','red'),lty=2, lwd=3)
abline(v=c(1,-1), col=c('red','red'),lty=2, lwd=3)

return(p)

}

apx_pi2(500)

## [1] 3.2

apx_pi2(10000)

## [1] 3.1664

apx_pi2(50000)

## [1] 3.14424

Mathematical Approach

Mathematically one can thought \(\pi\) as the area inside a circle of radios 1,
i.e. \[\pi = \int\int _A dx dy\] where, \(A=\{(x,y): |x|\leq1, |y|\leq1, x^2+y^2\leq1 \}\)

Search This Blog

Soumarya Basak

Pi Approximation

Approximation of Pi