Approximation Of Pi Using The Monte Carlo Method

i-Maze An Academic Festival Organized by IIPS, DAVV, Indore

Approximation of pi using the Monte Carlo Method

Dr. W.U.Khan

Palash Kar

Department of Computer Science Shri Govindram Institute of Technology and Science Indore

Abstract

Numerical methods that are known as Monte Carlo methods can be loosely described as statistical simulation methods, where statistical simulation is defined in quite general terms to be any method that utilizes sequences of random numbers to perform the simulation. This paper finds a way to approximate the value of pi using the Monte Carlo method. We have used a unit circle (radius=1) within a square with sides equal to 2.Using the probability of a point to lie inside the circle when picked randomly from the points inside the square we can approximate the value of pi. A c-program is also provided which models the method and provides a value of pi correct upto 3 to 4 places of decimal, the only limiting factor being the discrete representation of real Nos. in a computer rather than the method itself.

Introduction to Monte Carlo Methods Numerical methods that are known as Monte Carlo methods can be loosely described as statistical simulation methods, where statistical simulation is defined in quite general terms to be any method that utilizes sequences of random numbers to perform the simulation. Monte Carlo methods have been used for centuries, but only in the past several decades has the technique gained the status of a full-fledged numerical method capable of addressing the most complex applications. The name ``Monte Carlo'' was coined by Metropolis (inspired by Ulam's interest in poker) during the Manhattan Project of World War II, because of the similarity of statistical simulation to games of chance, and because the capital of Monaco was a center for gambling and similar pursuits. Monte Carlo is now used routinely in many diverse fields, from the simulation of complex physical phenomena such as radiation transport in the earth's atmosphere and the simulation of the esoteric subnuclear processes in high energy physics experiments, to the mundane, such as the simulation of a Bingo game or the outcome of Monty Hall's vexing offer to the contestant in ``Let's Make a Deal.'' The analogy of Monte Carlo methods to games of chance is a good one, but the ``game'' is a physical system, and the outcome of the game is not a pot of money or stack of chips (unless simulated) but rather a solution to some problem. The ``winner'' is the scientist, who judges the value of his

results on their intrinsic worth, rather than the extrinsic worth of his holdings. Statistical simulation methods may be contrasted to conventional numerical discretization methods, which typically are applied to ordinary or partial differential equations that describe some underlying physical or mathematical system. In many applications of Monte Carlo, the physical process is simulated directly, and there is no need to even write down the differential equations that describe the behavior of the system. The only requirement is that the physical (or mathematical) system be described by probability density functions (pdf's), which will be discussed in more detail later in this chapter. For now, we will assume that the behavior of a system can be described by pdf's. Once the pdf's are known, the Monte Carlo simulation can proceed by random sampling from the pdf's. Many simulations are then performed (multiple ``trials'' or ``histories'') and the desired result is taken as an average over the number of observations (which may be a single observation or perhaps millions of observations). In many practical applications, one can predict the statistical error (the ``variance'') in this average result, and hence an estimate of the number of Monte Carlo trials that are needed to achieve a given error. Figure 1 illustrates the idea of Monte Carlo, or statistical, simulation as applied to an arbitrary physical system. Assuming that the evolution of the physical system can be described by probability density functions (pdf's), then the Monte Carlo simulation can proceed by sampling from these pdf's, which necessitates a fast and effective way to generate random numbers uniformly distributed on the interval [0,1]. The outcomes of these random samplings, or trials, must be accumulated or tallied in an appropriate

manner to produce the desired result, but the essential characteristic of Monte Carlo is the use of random sampling techniques (and perhaps other algebra to manipulate the outcomes) to arrive at a solution of the physical problem. In contrast, a conventional numerical solution approach would start with the mathematical model of the physical system, discretizing the differential equations and then solving a set of algebraic equations for the unknown state of the system.

It should be kept in mind though that this general description of Monte Carlo methods

may not directly apply to some applications. It is natural to think that Monte Carlo methods are used to simulate random, or stochastic, processes, since these can be described by pdf's. However, this coupling is actually too restrictive because many Monte Carlo applications have no apparent stochastic content, such as the evaluation of a definite integral or the inversion of a system of linear equations. However, in these cases and others, one can pose the desired solution in terms of pdf's, and while this transformation may seem artificial, this step allows the system to be treated as a stochastic process for the purpose of simulation and hence Monte Carlo methods can be applied to simulate the system. Therefore, we take a broad view of the definition of Monte Carlo methods and include in the Monte Carlo rubric all methods that involve statistical simulation of some underlying system, whether or not the system represents a real physical process. To illustrate the diversity of Monte Carlo methods, Figure 2 lists applications that have been addressed with statistical simulation techniques. As can be seen, the range of applications is enormous, from the simulation of galactic formation to quantum chromodynamics to the solution of systems of linear equations.

Major Components of a Monte Carlo Algorithm Given our definition of Monte Carlo, let us now describe briefly the major components of a Monte Carlo method. These components comprise the foundation of most Monte Carlo applications, and the following sections will explore them in more detail. An understanding of these major components will provide a sound foundation for the reader to construct his or her own Monte Carlo method, although of course the physics and

mathematics of the specific application are well beyond the scope of this chapter. The primary components of a Monte Carlo simulation method include the following: • Probability distribution functions (pdf's) - the physical (or mathematical) system must be described by a set of pdf's. • Random number generator - a source of random numbers uniformly distributed on the unit interval must be available. • Sampling rule - a prescription for sampling from the specified pdf's, assuming the availability of random numbers on the unit interval, must be given. • Scoring (or tallying) - the outcomes must be accumulated into overall tallies or scores for the quantities of interest. • Error estimation - an estimate of the statistical error (variance) as a function of the number of trials and other quantities must be determined. • Variance reduction techniques - methods for reducing the variance in the estimated solution to reduce the computational time for Monte Carlo simulation • Parallelization and vectorization - algorithms to allow Monte Carlo methods to be implemented efficiently on advanced computer architectures.

History of Monte Carlo method The method is called after the city in the Monaco principality, because of a roulette, a simple random number generator. The name and the systematic development of Monte

Carlo methods date from about 1944. There are however a number of isolated and undeveloped instances on much earlier occasions. For example, in the second half of the nineteenth century a number of people performed experiments, in which they threw a needle in a haphazard manner onto a board ruled with parallel straight lines and inferred the value of PI = 3.14… from observations of the number of intersections between needle and lines. In 1899 Lord Rayleigh showed that a one-dimensional random walk without absorbing barriers could provide an approximate solution to a parabolic differential equation. In 1931 Kolmogorov showed the relationship between Markov stochastic processes and certain integro-differential equations. In early part of the twentieth century, British statistical schools indulged in a fair amount of unsophisticated Monte Carlo work. Most of this seems to have been of didactic character and rarely used for research or discovery. Only on a few rare occasions was the emphasis on original discovery rather than comforting verification. In 1908 Student (W.S. Gosset) used experimental sampling to help him towards his discovery of the distribution of the correlation coefficient. In the same year Student also used sampling to bolster his faith in his so-called t-distribution, which he had derived by a somewhat shaky and incomplete theoretical analysis. The real use of Monte Carlo methods as a research tool stems from work on the atomic bomb during the second world war. This work involved a direct simulation of the probabilistic problems concerned with random neutron diffusion in fissile material; but

even at an early stage of these investigations, von Neumann and Ulam refined this particular " Russian roulette" and "splitting" methods. However, the systematic development of these ideas had to await the work of Harris and Herman Kahn in 1948. About 1948 Fermi, Metropolis, and Ulam obtained Monte Carlo estimates for the eigenvalues of Schrodinger equation. In about 1970, the newly developing theory of computational complexity began to provide a more precise and persuasive rationale for employing the Monte Carlo method. The theory identified a class of problems for which the time to evaluate the exact solution to a problem within the class grows, at least, exponentially with M. The question to be resolved was whether or not the Monte Carlo method could estimate the solution to a problem in this intractable class to within a specified statistical accuracy in time bounded above by a polynomial in M. Numerous examples now support this contention. Karp (1985) shows this property for estimating reliability in a planar multiterminal network with randomly failing edges. Dyer (1989) establish it for estimating the volume of a convex body in M-dimensional Euclidean space. Broder (1986) and Jerrum and Sinclair (1988) establish the property for estimating the permanent of a matrix or, equivalently, the number of perfect matchings in a bipartite graph.

How to calculate

We have to look at the unit circle (radius=1) within a square with sides equal to 2 (see figure 1). Now if we pick a random point (x, y) were both x and y are between -1..1, the probability of that

this random point lies inside the unit circle is given as the proportion between the area of the unit circle and the square:

If we pick a random point N times and M of those times the point lies inside the unit circle, the probability of that a random point lies inside the unit circle is given as:

Where the dot indicates that this is a discreet distribution (because M and N are integers). But if N becomes very large (theoretically infinite), the two probabilities will become equal and we can write:

Programming As you can see in the final formula for Pi, the precision (number of digits) depends on how many times you pick a point. The greater N and M are the more digits you get. However, the discrete representation of real numbers in a computer is also a limiting factor for the precision of this method. Because we only have a finite number of points (1019 for 32 bits) to choose from, we actually approximate a rational number close to pi, rather than pi itself!

An actual c-program is provided in which we have modeled the method to approximate the value of pi. The value of pi obtained is correct upto 3 places of decimal with n=1000000 and with n=100000000 it is obtained correct upto 4 to 5 places of decimal. Higher accuracy may be obtained by taking a higher value of n constrained by the architecture of the computer. For statistics sake the value of pi obtained with n=100000 was 3.139640, that with n=1000000 was 3.140545 and the one obtained with n=100000000 was 3.141554 correct upto 4 decimal places.

The Program: #include<stdio.h> #include<math.h> #include #include #include<stdlib.h> #define LIMIT 100000000

//max value of n

void main(){ double mypi,x,y ; int minus ; // 1 for pos., 2 for neg. long int m,n ; time_t t; n = 0; m = 0;

0..1

0..1

srand((unsigned) time(&t)); clrscr(); do{ n++; x = (rand()/(double)RAND_MAX); //pick a random X between minus = 1+(rand()/(double)RAND_MAX);//pick a sign for X switch (minus) { case 1 : x = x; case 2 : x = -1*x; } y = rand()/(double)RAND_MAX; //pick a random Y between minus = 1+rand()/(double)RAND_MAX; //pick a sign for Y switch (minus) { case 1 : x = x; case 2 : x = -1*x; } if (pow(x,2)+pow(y,2) < 1) m++; //the point is inside the unit circle mypi = 4.0*m/n; //calculate pi if (n % (LIMIT / 1000) == 0) { gotoxy(1,1); printf("X:%f ",x); printf("Y:%f ",y); printf("n:%ld ",n);

}

}

printf("m:%ld ",m); printf(" - "); printf("PI:%f ",mypi);

} while (n != LIMIT); getch();