1.3. Monte Carlo \(\pi\)#

For this Exercise, we import some useful libraries to deal with mathematical operations and random numbers:

import random as r
import math as m
import matplotlib.pyplot as plt

from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets

import sys
sys.path.append("..")

# The Random Seed is set here, you will focus on it in the last question.
# For the moment it is sufficient to know that this is used to initialize
# the generation of random numbers (if not set a default value is used, 
# but it can be useful to know what is the seed for reproducibility)
r.seed(42)

import helpers
helpers.set_style()

Modify the following piece of code to implement a numerical estimation of \(\pi\) via a Monte Carlo approach. Note that you only need to modify the code in the section explicitly highlighted).

You will need to generate demo points using random extractions. There are different ways to do so, one possibility is to use the uniform function of the random package that we have imported. Refer to the package documentation to understand how this function works, and how to obtain the desired outcome with it.

def estimatePI(n_trials=10):

    successfulHits = 0  
    circleRadius = 0.5
    squareLength = 1.0

    inside_x,inside_y=[],[]
    outside_x,outside_y=[],[]
    
    ## Begin code to modify ##

    # Iterate for the number of darts.
    for i in range(0, n_trials):
        x = 0 # generate a demo point here from random extraction
        y = 0
        # check if point in circle here and decide whether 
        # to append the point to inside or outside list    
        # also don't forget to count the number of succesful hits    
     
    # Estimate Pi 
    pi = 0 # Change code here to get an estimate of pi
    ## End Code to modify ##
     
    print(f"{pi:.10f}")

    fig, ax=plt.subplots(1, figsize=(5,5))

    ax.scatter(inside_x,inside_y)
    ax.scatter(outside_x,outside_y)
    circle = plt.Circle((0,0), circleRadius, fill=False)
    ax.add_patch(circle)

    plt.show()

Try to run a single test of your code with few trials (n=100). You should get a square plot with points inside the circle in blue and points outside the circle in orange.

estimatePI(n_trials=100)

Exercise 4

Perform the \(\pi\) estimation for 1000, 100000 and 1000000 trials. Take a screenshot of these estimations and include them in your report (e.g of the python plots). What happens to the accuracy of the \(\pi\) estimation when going from 1000 to 1000000 trials and why?

Everything working correctly?

Then, lets get a bit more fancy and allow changing the number of trials interactively.

The below code provides an interactive widget with a slider that will redraw the plot once you change the slider.

interact(estimatePI, n_trials=widgets.IntSlider(min=2, max=100000, step=1000, value=1000))

Exercise 5

What happens to the estimation of \(\pi\) when the circle origin is changed? Why?

Exercise 6

What happens to the accuracy of the estimation when you increase the square size, or decrease the circle size? Is there an optimal ratio between the square side (\(l\)) and the circle diameter (\(d\))?

Exercise 7

In our code we set the random seed (r.seed()) once in the first cell, but then it is unchanged. It means that the first time you run the MC code, the random number set as seed was used for the extraction, but then the later extractions do not start from that seed anymore. What happens if you use the same seed for the pseudo random number generator (PRNG) each time you start the MC code? And what would happen if you would use the same seed for each random number extraction? You can read up more details in Appendix: Random and Pseudo-random Numbers. Hint: To get an idea, you can either edit the MC code, or you can also write a piece of code to check how the numbers generated depend on the seed set.

Molecular Dynamics and Monte Carlo

Monte Carlo \pi

1.3. Monte Carlo \(\pi\)#