EGR 103/Fall 2018/Lec 12

This page will go through the development of the program we wrote in class. Each step along the way, I will describe what functionality we wanted to add and you will be able to see the code that accomplished that goal.

Introduction

We will be writing a program that uses a Monte Carlo method to estimate \(\pi\). The program will simulate throwing some number of darts at a dart board of radius \(r\) mounted on a square back board that is \(2r\) on a side. If the \(x\) and </math>y</math> coordinates are each uniformly distributed between \(-r\) and \(r\), then the probability that a dart hits the dart board is the same as the ratio of the area of the dart board to the area of the back board - that is, \(p_{hit}=\frac{\pi r^2}{(2r)^2}=\frac{\pi}{4}\). We will thus estimate \(\pi\) by dividing the number of darts that hit the dart board by the total number of darts thrown and then multiplying that by 4.

Tasks

We broke the tasks up into several parts and wrote to code to accomplish those tasks. The main tasks were:

• Import modules
• Draw the board
• Throw a dart
• Throw lots of darts
• Plot the dart locations
• Determine how many darts land on the dart board
• Print how many darts land on the dart board and the estimates of \(\pi\)
• Refine how that is printed
• Plot the evolution of the estimates
• Have the value of a variable decide if we should plot while the simulation is running or only at the end

Importing modules

As always, start by importing the necessary modules. We will be doing math, using computational methods, and making graphs, so we will need:

1 import math as m
2 import numpy as np
3 import matplotlib.pyplot as plt

Drawing the board

Assuming a dart board of radius \(r\), the square that encloses it would be \(2r\) on a side. Since we will later be using the distance from the center of the circle to determine if a dart is on the dart board, it is mathematically easiest to have the middle of the board at (0, 0). Drawing the back board thus means drawing a square with vertices at (-r, -r), (r, -r), (r, r), and (-r, r). The the circle, we will need a list of x and y coordinates that defines a circle; the easiest way to do that is to use polar coordinates with a constant radius and an angle that goes from 0 to \(2\pi\). We need to know how many values to use in making the circle, so that will be another input to the function. We chose a default radius of 1 and a default number of points of 100:

 1 import math as m
 2 import numpy as np
 3 import matplotlib.pyplot as plt
 4 
 5 def draw_board(r=1, pts=100):
 6     plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
 7     theta = np.linspace(0, 2*m.pi, pts)
 8     x = r*np.cos(theta)
 9     y = r*np.sin(theta)
10     plt.plot(x, y, 'r-')

Test the code by calling it from the script:

 1 import math as m
 2 import numpy as np
 3 import matplotlib.pyplot as plt
 4 
 5 def draw_board(r=1, pts=100):
 6     plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
 7     theta = np.linspace(0, 2*m.pi, pts)
 8     x = r*np.cos(theta)
 9     y = r*np.sin(theta)
10     plt.plot(x, y, 'r-')
11 
12 if __name__ == '__main__':
13     plt.figure(1)
14     plt.clf()
15     draw_board()

The problem with this is that Python will use all the figure window and depending on the shape of the window, the square and the circle might look non-square and non-circular. To require that Python plot using the same scale along both axes, use the .axis('equal') method of plt:

 1 import math as m
 2 import numpy as np
 3 import matplotlib.pyplot as plt
 4 
 5 def draw_board(r=1, pts=100):
 6     plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
 7     theta = np.linspace(0, 2*m.pi, pts)
 8     x = r*np.cos(theta)
 9     y = r*np.sin(theta)
10     plt.plot(x, y, 'r-')
11     plt.axis('equal')
12 
13 if __name__ == '__main__':
14     plt.figure(1)
15     plt.clf()
16     draw_board()

You will now have a square back board and a circular dart board. As an aside, if you want to see something trippy, copy the following into the console after you have run your script:

r = 1
for k in range(10):
     draw_board(r)
     r = r*m.sqrt(2)

Throw a dart

Throwing a dart means figuring out the \(x\) and \(y\) coordinates for a dart. To determine the valid location of darts, we need to know the radius of the dart board. Given all that, adding a function to return the location of a thrown dart might look like:

 1 import math as m
 2 import numpy as np
 3 import matplotlib.pyplot as plt
 4 
 5 def draw_board(r=1, pts=100):
 6     plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
 7     theta = np.linspace(0, 2*m.pi, pts)
 8     x = r*np.cos(theta)
 9     y = r*np.sin(theta)
10     plt.plot(x, y, 'r-')
11     plt.axis('equal')
12     
13 def throw_darts(r=1):
14     x = np.random.uniform(-r, r)
15     y = np.random.uniform(-r, r)
16     return (x, y)
17 
18 if __name__ == '__main__':
19     plt.figure(1)
20     plt.clf()
21     draw_board()

You can test this in the console by throwing darts at a board of a given size (say, 10):

for k in range(10):
    print(throw_darts(10))

Make sure both coordinates are in the proper range.

Throw lots of darts

Now that you can throw a single dart, you might want to throw several darts. You should add variables to keep track of the total number of darts to throw and of the total number thrown:

 1 import math as m
 2 import numpy as np
 3 import matplotlib.pyplot as plt
 4 
 5 def draw_board(r=1, pts=100):
 6     plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
 7     theta = np.linspace(0, 2*m.pi, pts)
 8     x = r*np.cos(theta)
 9     y = r*np.sin(theta)
10     plt.plot(x, y, 'r-')
11     plt.axis('equal')
12     
13 def throw_darts(r=1):
14     x = np.random.uniform(-r, r)
15     y = np.random.uniform(-r, r)
16     return (x, y)
17     
18 if __name__ == '__main__':
19     darts_to_throw = 100
20     darts_thrown = 0
21     
22     plt.figure(1)
23     plt.clf()
24     draw_board()
25 
26     while darts_thrown < darts_to_throw:
27         darts_thrown += 1
28         x_d, y_d = throw_darts()

Plot the dart locations

Now that you are throwing darts, you can plot the locations at which they land. There are two fundamentally different ways to do this - you can plot each dart as it is thrown or you can wait until the end to plot all the darts. We are going to first look at how to plot each dart and watch the board fill up with darts. They key here will be a function that takes the dart coordinates and then updates the figure before moving on to the next command. Having the figure update right now requires drawing the plot and then adding a pause to let the plot...plot:

 1 import math as m
 2 import numpy as np
 3 import matplotlib.pyplot as plt
 4 
 5 def draw_board(r=1, pts=100):
 6     plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
 7     theta = np.linspace(0, 2*m.pi, pts)
 8     x = r*np.cos(theta)
 9     y = r*np.sin(theta)
10     plt.plot(x, y, 'r-')
11     plt.axis('equal')
12     
13 def throw_darts(r=1):
14     x = np.random.uniform(-r, r)
15     y = np.random.uniform(-r, r)
16     return (x, y)
17     
18 def plot_darts(x, y):
19     plt.figure(1)
20     plt.plot(x, y, 'b.')
21     plt.draw()
22     plt.pause(0.01)
23     
24 if __name__ == '__main__':
25     darts_to_throw = 100
26     darts_thrown = 0
27     
28     plt.figure(1)
29     plt.clf()
30     draw_board()
31 
32     while darts_thrown < darts_to_throw:
33         darts_thrown += 1
34         x_d, y_d = throw_darts()
35         plot_darts(x_d, y_d)