Goal

The goal of this lab is to effectively model and implement the trapezoid method of integration.

Overview

We create the function to calculate integrals using the trapezoid method and apply them to a few physical situations. We then compare them to some python built-in numerical integrators.

Functions

The first function we made was to implement the Trapezoidal Integration Method which takes in arrays of x and y values and returns the integral between xmin and xmax.

def trap(x,y):
  A=0
  for i in range(1,len(x)):
    A+=.5*(y[i-1]+y[i])*(-x[i-1]+x[i])
  return A

The second equation is to calculate the area under the half circle using the trapezoid method. We want to return the y values.

def mycircle(x,r=2):
  y=np.sqrt(r**2-x**2)
  return y

We then used both of these functions to calculate the area under the half circle.

x=np.linspace(-2,2,100)
y=mycircle(x,r=2)
A=trap(x,y)
print(f"The area under the curve is {A:.5f}")

Then we determined the percent error by using the following.

actual= (1/2)*(np.pi)*(2**2)

pError=(abs(actual-A)/actual)*100
print(f"Our percent error is {pError:.4f}")

It is import to note that the more steps you use in the linspace function the more accurate your number will be.

We then make a function for the Gaussian distribution which takes in x values, mean, and $\sigma$.

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

And then we integrate it between $-\sigma$ and $\sigma$, and use our previous equations using the following code.

sigma=1
x=np.linspace(-sigma,sigma,100)
fx=mygauss(x,mean=0,sigma=1)
#print(fx)

A=trap(x,fx)
print(f"{A:.4f}")
error=((abs(.682-A))/.682)*100
print(f"The percent error is {error:.4f}")

Using Python Built-In Functions

For this we use scipy.integrate.trapezoid function to perform our integration.

# import scipy.integrate.trapezoid
import scipy
# set up x and y arrays for the half circle
x=np.linspace(-2,2,100)
y=mycircle(x,r=2)
# call scipy.integrate.trapezoid to perform the integration
A=scipy.integrate.trapezoid(y,x)
# print your result and the expected answer in a formatted print statement
actual= (1/2)*(np.pi)*(2**2)
print(f" Our result is {A:.4f} while our expected value was {actual:.4f}")

These answers are very similar to those that we calculated making our own functions.

Calculating Area using the Monte Carlo Method and Plot

This function takes random points around the Area under the circle to calculate the area using the following code.

# your monte carlo function here
def mcintegrate(radius = 1, num = 1000, visualize = False):

  xpt = np.random.uniform(-radius, radius, num)
  ypt = np.random.uniform(0, radius, num)

  ptBelow = ypt <= mycircle(xpt,radius)
  ptAbove = ypt > mycircle(xpt,radius)

  areaCircle = np.sum(ptBelow) / num * (2 * radius ** 2)

  if visualize == True:
    plt.figure(figsize=(8, 8))

    # Plot the half circle
    x = np.linspace(-radius, radius, 100)
    plt.plot(x, mycircle(x), 'k-', label='Half Circle')

    # Plot points below and above the half circle
    plt.scatter(xpt[ptBelow], ypt[ptBelow], color='blue', label='Points Below')
    plt.scatter(xpt[ptAbove], ypt[ptAbove], color='red', label='Points Above')

  return areaCircle

And to call and plot it you use this code.

xval=np.linspace(-2,2,100)
yval=mycircle(xval)

areaval=trap(xval,yval)

areaValMonte = mcintegrate(2, visualize = True)

print(f"This is my calculated area value using my trap function, {areaVal}")

print(f"This is my calculated area value using my mcintegrate function, {areaValMonte}")

Which gives you this plot: