For a large number of throws, we see that:

hits area_of_quadrant

------ = ----------------

throws area_of_square

Some half-remembered geometry tells us that:

area_of_quadrant (1 / 4) * pi * r^2 pi

---------------- = ------------------ = --

area_of_square r^2 4

Or:

area_of_quadrant hits

pi = 4 * ---------------- = 4 * ------

area_of_square throws

I first solved this problem as an undergraduate sometime in 1994 as part of a Computational Physics module. Using FORTRAN.

The full explanation can be found here

#!/usr/bin/env python

import random

import math

# class representing a single throw of a dart

class Throw:

def __init__(self):

# generate two random coordinates and work out how far away we are from

# the origin

self.x = random.random()

self.y = random.random()

self.distance = self.distance()

def distance(self):

# the distance from the origin is the hypotenuse of a right-angled

# triangle with sides of length and x and y. Pythagoras told us that:

# distance = sqrt((x^2) + (y^2))

# which looks like this in python

return math.sqrt(self.x**2 + self.y**2)

# did we land inside the quadrant?

def is_a_hit(self):

return self.distance <= 1.0

# main class

class MonteCarlo:

def __init__(self):

# initialise everything

self.hits = 0

self.throws = 0

self.pi = 0

# this method is called on every throw

def increment(self, throw):

# we always clock up another throw

self.throws += 1

# and accumulate a hit if we scored

if throw.is_a_hit():

self.hits += 1

# then get a new value of pi

self.calculate_pi()

# explanation can be found here: http://icanhaz.com/calculatingpi

def calculate_pi(self):

self.pi = 4 * (float(self.hits) / float(self.throws))

# we use this in determining whether to print our status

def divides_by(self, number):

return float(self.throws) % number == 0

# represent the current state as a string

def __repr__(self):

return "Throws: %10d, Hits: %10d, Pi: %10f, Actual Pi: %10f" % (self.throws, self.hits, self.pi, math.pi)

# if we're called on the command line

if __name__ == '__main__':

# construct a new instance

m = MonteCarlo()

# loop forever

while 1:

# keep throwing darts

m.increment(Throw())

# only print on every 100000th iteration

if m.divides_by(100000): print m

Get the code here:

`svn co http://pikesley.homedns.org/svn/code/MonteCarlo`

## No comments:

## Post a Comment