14.5 Monte Carlo, Buffon s Needle, and p

---

 

	Java Number Cruncher: The Java Programmer's Guide to Numerical Computing By Ronald  Mak
	Table of Contents
	Chapter  14.   Generating Random Numbers

---

14.5 Monte Carlo, Buffon's Needle, and p

Besides being the name of a famous European gambling resort, Monte Carlo refers to a type of algorithm that uses random numbers to compute a desired value. ^[2] Monte Carlo algorithms can succeed when there are no simple computational formulas.

^[2] This bit of naming whimsy is attributed to John von Neumann.

We'll examine a Monte Carlo algorithm for computing the value of p . In Chapter 13, we saw several elegant formulas for p that can compute hundreds of thousands, if not more, digits of p relatively efficiently . Our Monte Carlo algorithm does not do nearly as well, but it is fascinating, nevertheless.

Here is how the algorithm, known as Buffon's needle, works. ^[3] Suppose you dropped needles , each of length 1, at random onto a sheet of paper ruled with parallel lines 1 unit apart. What is the probability that a needle crosses a line?

^[3] It is named after the person who devised the algorithm, Georges-Louis Leclerc, Comte de Buffon (1707 §C1788), a French mathematician and scientist.

As shown in Figure 14-1, three uniformly distributed random values represent each needle. The first two random values are the x and y coordinates of the needle's midpoint , from which we measure the perpendicular distance d to the nearest line. The third random value is the needle's angle of rotation q about its midpoint. We can see that the needle crosses the line whenever .

Figure 14-1. A dropped Buffon needle of length 1 that crosses a line. The three uniformly distributed random values are the coordinates x and y of the needle's midpoint and its angle of rotation q .

The value of d can range from 0 to , and the value of q can range from 0 to p . This is the rectangle in Figure 14-2, and its area is . The curve inside the rectangle is the function , and the shaded region under the curve is where . The area of the shaded region is

Figure 14-2. The probability that a needle crosses a line is 2/ p , the ratio of the shaded area under the curve to the total area of the rectangle.

The probability that a needle crosses a line is the ratio of the shaded area under the curve to the total area of the rectangle:

Therefore, if we continue to drop needles, the ratio of the number of needles that cross a line to the total number of needles approaches , from which we can compute the value of p .

Listing 14-3a shows class Needles , which represents the dropped needles. Method dropNext() uses the translation and rotation formulas we saw in Chapter 9.

Listing 14-3a Class Needles .

 package numbercruncher.program14_3; import java.util.Random; /**  * Implementation of Buffon's needles, which are randomly dropped  * onto a ruled sheet of paper.  */ class Needles {     /** x coord of one end of a needle */        private float x1;     /** x coord of the other end of a needle */  private float x2;     /** y coord of one end of a needle */        private float y1;     /** y coord of the other end of a needle */  private float y2;     /** angle of rotation about the midpoint */  private float theta;     /** paper's left edge */    private float xMin;     /** paper's width */        private float xWidth;     /** paper's bottom edge */  private float yMin;     /** paper's height */       private float yHeight;     /** count of dropped needles */      private int    count;     /** number that cross a line */      private int    crossings;     /** current computed value of pi */  private float  pi;     /** generator of uniformly-distributed random values */     private Random random = new Random(System.currentTimeMillis());     /**      * Constructor.      * @param xMin the paper's left edge      * @param xMax the paper's right edge      * @param yMin the paper's bottom edge      * @param yMax the paper's top edge      */     Needles(float xMin, float xMax, float yMin, float yMax)     {         this.xMin    = xMin;         this.xWidth  = xMax - xMin;         this.yMin    = yMin;         this.yHeight = yMax - yMin;     }     /**      * Return the x coordinate of one end of the needle.      * @return the coordinate      */     float getX1() { return x1; }     /**      * Return the x coordinate of one end of the needle.      * @return the coordinate      */     float getX2() { return x2; }     /**      * Return the y coordinate of one end of the needle.      * @return the coordinate      */     float getY1() { return y1; }     /**      * Return the y coordinate of the other end of the needle.      * @return the coordinate      */     float getY2() { return y2; }     /**      * Return the count of all the dropped needles.      * @return the count      */     int getCount() { return count; }     /**      * Return the number of needles that crossed a line.      * @return the number of crossings      */     int getCrossings() { return crossings; }     /**      * Return the current computed value of pi      * @return the value      */     float getPi() { return pi = (2f*count)/crossings; }     /**      * Return the error of the current computed value of pi      * @return the error      */     float getError() { return pi - ((float) Math.PI); }     /**      * Initialize.      */     void init()     {         count = crossings = 0;     }     /**      * Drop the next needle.      */     void dropNext()     {         ++count;         // Compute random values for the x and y coordinates of the         // needle's midpoint, and for the needle's angle of rotation.         float xCenter =  xWidth*random.nextFloat() + xMin;         float yCenter = yHeight*random.nextFloat() + yMin;         float theta   = (float) Math.PI*random.nextFloat();         float sin = (float) Math.sin(theta);         float cos = (float) Math.cos(theta);         // Rotate about the origin a 1-unit length needle on         // the x axis with endpoints at -0.5 and +0.5.         x1 = -0.5f*cos;         x2 = +0.5f*cos;         y1 = -0.5f*(-sin);         y2 = +0.5f*(-sin);         // Translate the needle to its location.         x1 += xCenter;         x2 += xCenter;         y1 += yCenter;         y2 += yCenter;         // Does the needle cross a line?         if (Math.floor(y1) != Math.floor(y2)) ++crossings;     } } 

Program 14 §C3 demonstrates the Buffon's needle algorithm. See Listing 14-3b. We see from the program's output that this is not really a good way to compute p ?ait converges very slowly and erratically.

Listing 14-3b A demonstration of Buffon's needle algorithm.

 package numbercruncher.program14_3; import numbercruncher.mathutils.AlignRight; /**  * PROGRAM 14-3: Buffon's Needle  *  * Demonstrate how we can calculate the value of pi using the  * Monte Carlo technique and Buffon's needle algorithm.  */ public class Buffon {     private static final int MAX_NEEDLES = 10000000;    // 10 million     private static final int GROUP_SIZE  = MAX_NEEDLES/20;     /**      * Main.      * @param args the array of program arguments      */     public static void main(String args[])     {         Needles    needles = new Needles(0, 0, 0, 10.5f);         AlignRight ar      = new AlignRight();         ar.print("Needles", 15);  ar.print("Crossings", 15);         ar.print("Pi", 15);       ar.print("Error", 15);         ar.underline();         do {             // Drop a group of needles.             for (int i = 0; i < GROUP_SIZE; ++i) needles.dropNext();             ar.print(needles.getCount(), 15);             ar.print(needles.getCrossings(), 15);             ar.print(needles.getPi(), 15);             ar.print(needles.getError(), 15);             ar.println();         } while (needles.getCount() < MAX_NEEDLES);     } } 

Output:

 Needles      Crossings             Pi          Error ------------------------------------------------------------          500000         318282      3.1418679   2.7513504E-4         1000000         636433      3.1425147   9.2196465E-4         1500000         954980       3.141427  -1.6570091E-4         2000000        1273603       3.140696   -8.966923E-4         2500000        1591847      3.1410053   -5.874634E-4         3000000        1910364      3.1407628   -8.299351E-4         3500000        2229064      3.1403315  -0.0012612343         4000000        2547078       3.140854   -7.388592E-4         4500000        2865942      3.1403286  -0.0012640953         5000000        3184564      3.1401472  -0.0014455318         5500000        3502476      3.1406353    -9.57489E-4         6000000        3820892      3.1406279     -9.6488E-4         6500000        4139230      3.1406808  -9.1195107E-4         7000000        4457471      3.1407945   -7.982254E-4         7500000        4775901      3.1407685    -8.24213E-4         8000000        5093656      3.1411622  -4.3058395E-4         8500000        5411821      3.1412716  -3.2114983E-4         9000000        5730361      3.1411633  -4.2939186E-4         9500000        6048652      3.1411958  -3.9696693E-4        10000000        6366532      3.1414278  -1.6498566E-4 

Screen 14-3 shows a screen shot of the interactive version of Program 14 §C3 in action.

Screen 14-3. Buffon's needle algorithm in action.

---

Top