In Class Exercise: Monte Carlo Integration

Write a programme following the method of Monte Carlo Integration to calculate

$$ I = \int_0^2 \sin^2 [\frac{1}{x(2-x)}] dx. $$

As you will need to calculate $f(x) = \sin^2 [\frac{1}{x(2-x)}]$ many times please write a user defined function for this part of your programme.



In [6]:

    
import numpy as np
import matplotlib.pyplot as plt
#This just needed for the Notebook to show plots inline. 
%matplotlib inline 

#Define the function
def f(x):
    fx = (np.sin(1/(x*(2-x))))**2
    return fx

#Integrate the function from x=0-2
#Note that you need to know the maximum value of the function
#over this range (which is y=1), and therefore the area of the box
#from which we draw random number is A=2. 
N=100000
k=0
for i in range(N):
    x=2*np.random.random()
    y=np.random.random()
    if y<f(x):
        k+=1
A=2.        
I=A*k/N
print("The integral is equal to I = ",I)









    



The integral is equal to I =  1.4503

In Class Exercise

Returning to the programme you just wrote to calculate the integral of $f(x) = \sin^2 [\frac{1}{x(2-x)}]$, please now add an estimate of the accuracy of your integration. Using

$$\sigma_I = \frac{\sqrt{I(A-I)}}{\sqrt(N)} $$



In [7]:

    
#Calculate the error: 

sigmaI = np.sqrt(I*(A-I))/np.sqrt(N)
print("The integral is equal to I = ",I)
print("The error on the integral is equal to sigmaI = {0:.4f}".format(sigmaI))









    



The integral is equal to I =  1.4503
The error on the integral is equal to sigmaI = 0.0028

In Class Exercise

Now calculate the integral of $f(x) = \sin^2 [\frac{1}{x(2-x)}]$ using the mean value method, and compare the accuracy from 10000 sample to what you obtained previously.



In [8]:

    
#Draw N values from the distribution, and calculate their mean to 
#use the mean method for integration. 
xvalues=[]
fvalues=[]
for i in range(N):
    x=2*np.random.random()
    y=np.random.random()
    if y<f(x):
        fvalues.append(f(x))
        xvalues.append(x)

fmean=np.mean(fvalues)
Imean = 2*fmean
Imeansigma = 2*np.sqrt(np.var(fvalues))/np.sqrt(len(fvalues))

print("Mean Method: ")
print("The integral is equal to I = {0:.4f}".format(Imean))
print("The error on the integral is equal to sigmaI = {0:.4f}".format(Imeansigma))
print("The percent error is: {0:.2f} ".format(100*Imeansigma/Imean))
print("**********************")
print("Compare to the `hit or miss` Monte Carlo Method: ")
print("The integral is equal to I = {0:.4f}".format(I))
print("The error on the integral is equal to sigmaI = {0:.4f}".format(sigmaI))
print("The percent error is: {0:.2f} ".format(100*sigmaI/I))









    



Mean Method: 
The integral is equal to I = 1.6432
The error on the integral is equal to sigmaI = 0.0011
The percent error is: 0.07 
**********************
Compare to the `hit or miss` Monte Carlo Method: 
The integral is equal to I = 1.4503
The error on the integral is equal to sigmaI = 0.0028
The percent error is: 0.19



In [9]:

    
#Using mpl_style file.
import mpl_style
plt.style.use(mpl_style.style1)

plt.plot(xvalues,fvalues,'.')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.show()



In [18]:

    
np.mean?



In [ ]: