Java: How to approximate Pi with the Monte Carlo simulation

You can approximate the mathematical constant Pi with a Monte Carlo simulation which is a method based on lots of random experiments that allows to make statements about the simulated objects. [1]

Computation of Pi

Before I present my Java program for the Monte Carlo simulation, I would like to explain some mathematical basics. The area of a circle is computed with the radius r:

A = \pi r^2

If you use the unit circle with the radius r = 1, you can ignore the radius. Next, we want to be only focused on a quadrant of this unit circle and the unit square. That is the formula for Pi:

\pi = \frac{A}{4}

This chart displays a quadrant in a unit square. The green points are inside the quadrant and the red ones are outside.

All points are inside the quadrant if the following equation is fulfilled:

x^2 + y^2 \leq 1

The relative frequency is determined with the quotient of the area of the quadrant and the unit square:

\text{relative frequency (h)} = \frac{\text{points inside the quadrant (q)}}{\text{points inside the square (n)}}

q includes the amount of attempts in which a point is inside the quadrant, and n contains all points in the unit square which means both the points inside and outside the quadrant.

h = \frac{q}{n} \approx \frac{\pi}{4}

By the way, you can also determine the fraction

\frac{\pi}{4}

with the Leibniz formula for Pi. [2] This is the value for Pi with the first 24 digits:

\pi \approx 3.141592653589793238462643 \text{ [3]}

After this short excursion, all mathematical requirements are known for a Java program.

Java program for Pi approximation with the Monte Carlo method

It is a popular exercise in secondary schools and colleges to write a program that computes a value for n that approximates Pi in a good way. One approach is to compute Pi with different values for n:

// calculation for n
private void calculateForN(int n) {
 for (int i = 0; i < n; i++) {
   // method getRand() provides random values
   double x = this.getRand();
   double y = this.getRand();

   if (x * x + y * y <= 1) {
    // point is inside the circle
    this.in++;
   } else {
    // point is outside the circle
    this.out++;
   }
 }
 // print results in the console
 this.print(n);
}

The disadvantage of this approach is that you calculate Pi for a fixed n, e.g. for 1,000 random points. Depending on the current random numbers for the 1,000 points you get a approximation of Pi of for example 3.164.

\pi = 3.164 \text{ with n = 1,000}

The picture above shows the results for different n in one simulation. The bigger n, the better the approximation of Pi. This observation can be explained with the Law of large numbers (LLN). [4]

But a better approximation of Pi can be also reached with for example 904 random points although the value of n is small: 3.14159292.

\pi = 3.14159292 \text{ with n = 904}

And this insight moves us on to another approach in which an Epsilon is defined that sets an upper limit for an acceptable difference between the calculated value for Pi and Pi:

This second approach leads to a smaller n and more exact value for Pi because the computation stops when the difference between the calculated Pi and Pi is less than Epsilon.

// calculate PI for a certain Epsilon
 public void calcEpsilon() {
  // initiate values
  this.resetValues();
  // Boolean is used to determine whether the calculated value for PI is acceptable
  Boolean done = false;
  // compute until done is true
  while (!done) {
   double x = this.getRand();
   double y = this.getRand();

   if (x * x + y * y <= 1) {
    this.in++;
   } else {
    this.out++;
  }

  // calculate current PI
  this.piCalc = this.calcPi();
  // compare error with epsilon
  error = this.getError(this.piCalc);

  // epsilon is greater than error
  if (error == -1) {
   // error is less than epsilon
   done = true;
  }

  // exit if there cannot be determined a Pi with an error less than epsilon
  if ((this.in + this.out) > 100000000) {
   break;
  }
 }
 System.out.println("Error (" + this.pi.subtract(this.piCalc).abs().setScale(15, RoundingMode.CEILING) + ") < epsilon (" + this.epsilon.setScale(10, RoundingMode.CEILING) + ")? " + done);
 print((int) (this.in + this.out));
}

// compare error with epsilon
private int getError(BigDecimal currentPi) {
 // currently calculated PI minus PI with 24 digits
 BigDecimal diff = currentPi.subtract(this.pi).abs();
 int error = diff.compareTo(this.epsilon);
 // error = 1: diff > epsilon
 // error = 0: diff = epsilon
 // error = -1: diff < epsilon
 return error;
}

Helpful content

[1] Monte Carlo Estimate for Pi (wolfram.com),
Estimating Pi using the Monte Carlo Method (101computing.net, Python),
Estimating the value of Pi using Monte Carlo (www.geeksforgeeks.org, C++)
[2] Leibniz formula for Pi (Wikipedia)
[3] Pi (Wikipedia)
[4] Law of large numbers (Wikipedia)
[5] BigDecimal (oracle.com)
BigDecimal Class in java (geeksforgeeks.org)

Java program for Monte Carlo simulation for Pi

This is the entire source code of my small Java program for a Monte Carlo simulation for Pi. I decided to use the data type BigDecimal instead of double because BigDecimal represents numbers better. [5]

import java.math.BigDecimal;
import java.math.RoundingMode;
import java.util.Random;

public class MonteCarlo {
 private int in = 0;
 private int out = 0;
 // Pi with 24 digits
 private BigDecimal pi = new BigDecimal(3.141592653589793238462643);
 private BigDecimal piCalc = new BigDecimal(0.0);
 private BigDecimal epsilon = new BigDecimal(0.00001);
 private int error = 1;

 // get random value
 private double getRand() {
  Random rd = new Random();
  return rd.nextDouble();
 }

 // get amount of points inside the circle
 private int getIn() {
  return this.in;
 }

 // get amount of points outside the circle
 private int getOut() {
  return this.out;
 }

 // calculate PI
 private BigDecimal calcPi() {
  // multiply the result with 4 because the simulation was only
  BigDecimal piTmp = BigDecimal.valueOf(4.0 * ((double) this.getIn() / ((double) this.getIn() + (double) this.getOut())));
  return piTmp;
 }

 // get Pi with 24 digits
 private BigDecimal getPi() {
  return this.pi;
 }

 // reset values to ensure valid results
 private void resetValues() {
  this.in = 0;
  this.out = 0;
  this.piCalc = new BigDecimal(0.0);
  this.error = 1;
 }

 // print the results to the console
 private void print(int n) {
  System.out.println("Calculation with: " + n);
  System.out.println("Points in: " + this.getIn() + ", points out: " + this.getOut());
  System.out.println("Pi calculated: " + this.calcPi() + ", In+out: " + (this.getOut() + this.getIn()));
  System.out.println("Difference: Pi calculated (" + this.calcPi() + ") - Pi (" + this.pi.setScale(10, RoundingMode.CEILING)	+ "): " + this.calcPi().subtract(this.getPi()).abs().setScale(10, RoundingMode.CEILING));
  System.out.println("--------------- " + "\n");
 }

 // start calculation with a certain value n
 public void startCalc(int n) {
  System.out.println("Monte Carlo simulation has started with n = " + n);
  // reset the values because of a previous calculation
  this.resetValues();
  // start calculation
  this.calculateForN(n);
 }

 // calculation for n
 private void calculateForN(int n) {
  for (int i = 0; i < n; i++) {
   double x = this.getRand();
   double y = this.getRand();
   if (x * x + y * y <= 1) {
    // point is inside the circle
    this.in++;
   } else {
    // point is outside the circle
    this.out++;
   }
  }
  this.print(n);
 }

 // compare error with epsilon
 private int getError(BigDecimal currentPi) {
  // currently calculated PI minus PI with 24 digits
  BigDecimal diff = currentPi.subtract(this.pi).abs();
  int error = diff.compareTo(this.epsilon);

  // error = 1: diff > epsilon
  // error = 0: diff = epsilon
  // error = -1: diff < epsilon
  return error;
 }

 // calculate PI for a certain Epsilon
 public void calcEpsilon() {
  // initiate values
  this.resetValues();
  // Boolean is used to determine whether the calculated value for PI is acceptable
  Boolean done = false;

  // compute until done is true
  while (!done) {
   double x = this.getRand();
   double y = this.getRand();
   if (x * x + y * y <= 1) {
    this.in++;
    } else {
    this.out++;
   }
   // calculate current PI
   this.piCalc = this.calcPi();
   // compare error with epsilon
   error = this.getError(this.piCalc);

   // epsilon is greater than error
   if (error == -1) {
    // error is less than epsilon
    done = true;
   }
   // exit if there cannot be determined a Pi with an error less than epsilon
   if ((this.in + this.out) > 100000000) {
    break;
   }
  }
   System.out.println("Error (" + this.pi.subtract(this.piCalc).abs().setScale(15, RoundingMode.CEILING)	+ ") < epsilon (" + this.epsilon.setScale(10, RoundingMode.CEILING) + ")? " + done);
  print((int) (this.in + this.out));
 }

 public static void main(String[] args) {
  // create an object
  MonteCarlo mc = new MonteCarlo();
  // calculate PI with n random points
  mc.startCalc(1000);
  mc.startCalc(10000);
  mc.startCalc(1000000);
  // calculate PI for a certain epsilon
  mc.calcEpsilon();
 }
}