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:
Solving for pi yields:
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.
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 %
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.