Learntofish's Blog

A blog about math, physics and computer science

Calculating Pi with the Monte Carlo method

Posted by Ed on October 13, 2010

We can calculate an approximate value for pi by using the Monte Carlo method. This probabilistic method relies on a random number generator and is described below. At the end of the post there is an excellent video by Kevin Wallenstein.


The idea is to approximate an area by counting dots. As an example see the image above. We want to determine the area of the smaller square compared to the bigger one. Obviously, the ratio between both areas is 1/4. Thus, if we randomly distribute 1000 dots then we expect about 250 dots to lie within the smaller square.

Now, consider the image with the circle. We know that the area of the circle is pi*r^2 and that of the square (2r)^2. The ratio of both areas is:

\frac{A_{circle}}{A_{square}} = \frac{\pi r^2} {(2r)^2} = \frac{\pi r^2}{4r^2} = \frac{\pi}{4}

Solving for pi yields: \pi = 4 * \frac{A_{circle}}{A_{square}}

This means if we can calculate the ratio A_circle / A_square on the right hand side of the equation we also get the value of pi. What we will do is approximate this ratio by choosing points (x,y) randomly in the square. (x,y) will either lie within or outside the circle.

circle_dots

The algorithm in pseudo-code is given below where for simplicity we have chosen the radius to be r=1 and examine only the first quadrant. Choosing only the first quadrant does not change the ratio. The code below does the following:
– Choose 1000 dots randomly.
– If the dot lands within the circle count the dot by increasing the variable “circleArea” by 1.
– At the end compare the number of dots within the circle with the total number of dots.

circleArea = 0
squareArea = 0
for(i=1 to 1000):
            x = random value from [0,1]
            y = random value from [0,1]
            if (x,y) within the circle:
                     circleArea = circleArea + 1
            squareArea = squareArea + 1
pi = 4.0*circleArea/squareArea
print pi

To check whether (x,y) lies within the circle (including the circumference) use the following function:

withinCircle(x,y)
    if(x^2+y^2<=1):
        return True
    else:
        return False

Below is an actual code in Python:

import random
import math

def withinCircle(x,y):
    if(x**2+y**2<1):
        return True
    else:
        return False

def main():
    circleArea = 0
    squareArea = 0
    pi = 0
    for i in range(0,100000):
        x = random.random()
        y = random.random()
        if(withinCircle(x,y)==1):
                   circleArea=circleArea+1
        squareArea=squareArea+1
    pi = 4.0*circleArea/squareArea
    print "Approximate value for pi: ", pi
    print "Difference to exact value of pi: ", pi-math.pi
    print "Error: (approx-exact)/exact=", (pi-math.pi)/math.pi*100, "%"

Run the code in Python by typing “main()” in the console and as an example output you will get:

Approximate value for pi:  3.056
Difference to exact value of pi:  -0.0855926535898
Error: (approx-exact)/exact= -2.72449878223 %

Type “main()” several times in the console. Each time you will get another value, but this is no surprise since our algorithm depends on a random number generator for x and y. You can also experiment by increasing the number of dots from 1000 to 10000. You will notice that the approximation becomes better:

Approximate value for pi:  3.152
Difference to exact value of pi:  0.0104073464102
Error: (approx-exact)/exact= 0.331276125131 %

References:
The video below by Kevin Wallenstein is an excellent explanation on Monte Carlo simulations. The beginning of the video is black but after 35 seconds you will see the graphics.

The graphics with the circles and dots were created with the excellent freeware program Inkscape.

Advertisements

6 Responses to “Calculating Pi with the Monte Carlo method”

  1. Younes said

    Thank you for sharing this useful information.
    Here is shortened Python code:

    #---------------------------------------------
    import random as r
    a = 0; c = 0
    for i in range(10000):
    x,y = r.random(), r.random()
    c+= int((x**2.0+y**2.0)<1)
    a+=1
    print 4.0*c/a-3.14159265358979
    #---------------------------------------------

  2. […] Click here to read more about this topic. […]

  3. […] from learntofish has a post with his implementation of the algorithm in […]

  4. […] have written a example of Mesos Framework by python. It simply calculate “Pi” by using mento-carlo algorithm The whole source code is at […]

  5. […] Monte Carlo algorithm for Pi value. […]

  6. […] Monte Carlo algorithm for Pi value. […]

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

 
%d bloggers like this: