Calculating "pi" by Monte Carlo Integration

Calculating π by Monte Carlo Integration

The Monte Carlo method is used to solve many complex problems for which an exact solution is not known. The name "Monte Carlo" refers to gambling, and in Monte Carlo methods many random outcomes, or guesses, are made and evaluated. π =3.14159..., the ratio of a circles diameter to its circumference, is numerically not exactly known. While there are far better ways to estimate π, I wrote this program to start learning C/C++ programming and play with Monte Carlo methods. It uses the relationship between π and the area of a circle, A=πr². An elegant way to visualize this relationship is to imagine slicing a circle into smaller and smaller pieces, then rearranging these pieces into a rectangle. The area of the rectangle is equal to the area of the circle and is equal to its width multiplied by its length, or A=πr².

If we set the radius to 1 then the area becomes equal to π. The trick with Monte Carlo integration is to find a function that includes your unknown function of interest, then guess values and evaluate if they are included in your function of interest. The fraction that are contained are proportional to the fraction of the area of the unknown function in your known functions area. I considered a quarter of a circle with a radius=1 that is included in a square with a side=1. Then random x, y coordinates are guessed and the distance from the center of the circle is determined by Pythagoreans' theorem, d=sqrt(x²+y²). If this distance is less than 1 they are in the quarter circle. The fraction in the circle vs. the square, or total, is multiplied by 4 to get the area of the whole circle and therefore π.

The following graph illustrates the diminishing returns in accuracy achived by increasing sample size. 12 independent points are plotted for each number of repititions. As sample size inceases the scatter diminishes and the π estimates rapidly approach a single value. But it becomes increasingly more difficult to "iron out" the remaining errors.

Beckman, Petr. A History of π (PI). 1973. St. Martin's Press. New York. is an entertaining read on math history pertaining to π as well as a dose of the authors political views. amazon.com Barnes&Noble

Download a ziped microsoft DOS/windows executable circle_pi.zip 26.8 KB
I also found that this executable, compiled for microsoft, also runs with wine on linux. Just download, unzip it, and type "wine circle_pi.exe" at the command line.

C++ Source Code:

//
// by Floyd Reed July 10, 2002
// after writing this I found a website 
// http://www.daimi.aau.dk/~u951581/pi/MonteCarlo/pimc.html
// that uses almost identical reasoning.  The only real difference
// is I used a quarter of a circle to make it easier to program.
//

#include <stdio.h>
#include <iostream.h>
#include <stdlib.h>	// required for rand function
#include <time.h>	// required for srand seed with time
#include <math.h>	// required for square root calculation

int main(int nNumberofArgs, char* pszArgs[])
{
	srand (time(NULL));		// uses time to seed random generator
	int nCArea = 0;
	int nGuess = 0;
	double dPi;
	cout << "CirclePi, version 2\n";
	cout << "This program attempts to calculate pi by a\n";
	cout << "  very inefficient 'guessing'  or 'Monte Carlo'\n";
	cout << "  method using the relationship between pi and\n";
	cout << "  the area of a circle.\n";
	cout << "Enter number of guesses to make, the more the\n";
	cout << "  better the estimate: ";
	cin >> nGuess;
	for (int nI=1; nI<=nGuess; nI++)
	{
		// generate x coordinate
		float fRandomx = rand();		// creates a random number, rand only returns integers so must convert to float to get 0 - 1 range
		float fTotalx = RAND_MAX;		// maximum value of random number
		float fX = fRandomx/fTotalx;	// crates a random zero to 1 fraction
		// generate y coordinate
		float fRandomy = rand();		// creates a random number, rand only returns integers so must convert to float to get 0 - 1 range
		float fTotaly = RAND_MAX;		// maximum value of random number
		float fY = fRandomy/fTotaly;	// crates a random zero to 1 fraction
		// calculate distance
		double dDistance = sqrt((fX*fX)+(fY*fY));	// distance from (0,0) equation
		//cout << "x=";   // you can "turn these on" to see all the calculations
		//cout << fX;
		//cout << "\ty=";
		//cout << fY;
		//cout << " \td=";
		//cout << fDistance;
		if (dDistance <= 1) {nCArea++;}		// counts distances within the quarter circle
		double dCArea = nCArea;
		double dI = nI;
		dPi = 4*dCArea/dI;
		//cout << " \tpi=";
		//cout << fPi;
		//cout << "\n";
	}
	cout << " = pi = \n";			// for some reason this prints after the printf command
	printf("\n%.6f", dPi);
	cout << "\n(I tried it for 1 billion - it took 4 minutes\n  and estimated pi as 3.141526)\nType 1 and return to close.";
	cin >> nGuess;
	return 0;
}

C Programming | Floyd's Home | MountainSmoke.com