Lab 07 - Curlss/Computational_physics GitHub Wiki

The purpose of this lab was to teach us how to do numerical integration so that we could write code to use the trapezoid method of integration. Then, we wanted to use this to create and run a Monte Carlo integrator that would predict the half-circle area and show the sampled spots. We would be able to understand half-circle area and gaussian distribution through this study.

In the lab, we need to write code for a simple  integrator so that we can test against scientific data and a library of data. We now know more about a statistical integration method and can better understand it. These methods work, and we learned how to use them. We also learned how their mistakes work and why one way will be better than another.

Lab

We wrote a function to use the trapezoid integration method. It ran over each set of x_i and x_i+1 to find the area of the trapezoid with heights y_i and y_i+1 and add them up. This set up a simple integrator that could be used for many things, which we then used in later parts of this lab.

# write your function here
def trap(y, x):
   integral = 0
   for i in range(len(x)-1):
      integral += (x[i+1] - x[i]) * (y[i] + y[i+1]) / 2
   return integral

To get the y values for our half circle, we then defined a function. Then, we used the x and y values to do the integral and find the area under the half circle.

# write your function mycircle here
def mycircle(x, r=2):
  # r^2 = x^2 + y^2
  y = np.sqrt(r**2 - x**2)
  return y

Calculate the area under the half circle using integration

  • In the next cell, set up an array of x values with $-2 \le x \le 2$, with 100 divisions between $-2$ and $2$.

  • Calculate the corresponding y values using your mycircle function.

  • Then use your trap function to calculate the area under the circle.

  • Print the value of your area using a formatted print statement. Print five decimal places.

x = np.linspace(-2, 2, 100)
y = mycircle(x)
area = trap(y, x)
print("The area under the half circle is {:.5f}".format(area))
The area under the half circle is 6.27644

We then defined a function that can integrate and measure the gaussian distribution. Our code and how well it works, show that even a simple trapezoid method can handle functions that change effortlessly, like the gaussian distribution.

def mygauss(x, mean, sigma):
  f = (1 / (np.sqrt(2 * np.pi) * sigma)) * np.exp(-0.5 * ((x - mean) / sigma)**2)
  return f

Integrate your gaussian

  • write code to integrate the gaussian between:

  • $ -\sigma \le x \le \sigma $

  • Compare your results with the expected answer, and calculate your percent error.

  • If you are not sure what the expected answer is, look here

mean = 0
sigma = 1
x = np.linspace(-sigma, sigma, 100)
y = mygauss(x, mean, sigma)
area = trap(y, x)
print("The area under the gaussian is {:.5f}".format(area))

accepted_value = 0.682689492137086
p = ((accepted_value - area) / accepted_value*100)
print(math.isclose(0.682689492137086,area, abs_tol=0.01))
print("The percent error is {:.2f}%".format(p))

As we did before, we use built-in methods in Python to finish the half-circle area integration and gaussian distribution integral. By doing that, we saw how possible and likely it is to use these functions we made instead of the known, optimized Python tools.

# import scipy.integrate.trapezoid
import scipy.integrate as trapezoid
# set up x and y arrays for the half circle
x = np.linspace(-2, 2, 100)
y = mycircle(x)
# call scipy.integrate.trapezoid to perform the integration
area = trapezoid.trapezoid(y, x)
# print your result and the expected answer in a formatted print statement
print("The area under the half circle is {:.5f}".format(area))
The area under the half circle is 6.27644

# import scipy.integrate.trapezoid
import scipy.integrate as trapezoid
# set up x and y arrays for the gaussian
x = np.linspace(-sigma, sigma, 100)
y = mygauss(x, mean, sigma)
# call scipy.integrate.trapezoid to perform the integration between $-1\sigma \le x \le 1\sigma$.
area2 = trapezoid.trapezoid(y, x)
# print your result and the expected answer in a formatted print statement
print("The area under the gaussian is {:.5f}".format(area2))
The area under the gaussian is 0.68267

Finally the code is used to show the true border curve along with randomly chosen places that are colored based on whether they are inside or outside the half-circle. This shows how chance can get close to an area and how the quality gets better as you take more samples.

#your monte carlo function here
#setting up arrays
xrandom = np.random.uniform(low=-2, high=2, size=100000)
yrandom = np.random.uniform(low=0, high=2, size=100000)

#setting up function
def mcintegrate (xrandom, yrandom, visualize = False):
  # get min and max values of xrandom
  xmin = np.min(xrandom)
  xmax = np.max(xrandom)
  # get min and max values of yrandom
  ymin = np.min(yrandom)
  ymax = np.max(yrandom)
  # test if yrandom < mycirc(xrandom)
  goodflag = yrandom < mycircle(xrandom)
  Ngood = np.sum(goodflag)
  Ntotal = len(xrandom)

  integral = Ngood/Ntotal * (xmax-xmin)*(ymax-ymin)

  if visualize == True:
    plt.plot(xrandom[goodflag],yrandom[goodflag],'b.')
    plt.plot(xrandom[~goodflag],yrandom[~goodflag],'r.')

    x = np.linspace(-2,2,100)
    y = mycircle(x, r=2)
    plt.plot(x, y, 'k')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title('Monte Carlo Integration')
    plt.text(0, 1.5, 'Area = {:.5f}'.format(integral))
    plt.show()

  return integral

# your code to calculate the area of a half circle using monte carlo
mcarea = mcintegrate(xrandom, yrandom, visualize=True)

image