Goal

We want to approximate the value of an integral of a function. To do so, we utilize the trapezoid integration method. The integral can also be calculated using the Monte Carlo method, and it provides a more visual depiction as well.

Overview

The trapezoidal integration method involves summing the trapezoidal area under a function for each $\Delta x$ between the range: $x_i$ to $x_f$. Each area is calculated via the following formula:

$$A = \frac{f(x_i)+f(x_i+\Delta x)}{2}\Delta x$$

We can create a function that utilizes this method, and we use it to calculate the area of a half-circle and the integral of the Gaussian distribution (Normal distribution). Later, we apply the Monte Carlo method to our half-circle example as an alternate way to determine its area. The Monte Carlo method involves generating random values for $x$ and $y$; those can then be plotted, and we can estimate the area based on the number of points that fall within the shape. This is given by:

$$\frac{N_{within}}{N_{total}}A_{total}$$

$A_{total}$ in this case refers to the total area of the plot.

This lab follows the following sequence:

• Function for Trapezoid Integration Method
• Area of Half-Circle
• Integral of Gaussian Distribution
• Compared to Built-In Functions
• Monte Carlo

Function for Trapezoid Integration Method

Given arrays of x and y, we can define a function to calculate the area using the formula detailed above.

def trap(x,y):
  z=0
  for i in range(1,len(x)):
    z+=(y[i-1]+y[i])/2*(x[i]-x[i-1])
  return z

Area of Half-Circle

Given an array of x values, we define a function to calculate exact y values corresponding to a circle of specified radius.

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

To test our trapezoid function, we make an array for x for a number of values between -r and +r, and use the mycircle function to get an array for y.

x = np.linspace(-2,2,100)
y = mycircle(x)

a = trap(x,y)
print(f"Area of half of a circle of 2 radius: {a:.5f}")

This returns the value: 6.27644. This can be compared to the actual value: $2\pi$, which evaluates to about 6.283: about a 0.1% error between the two. If we increase the number of steps for x, the accuracy improves.

Gaussian Distribution

The Gaussian function is given by:

$$f(x) = \frac{1}{\sqrt{2\pi}\sigma}e^{\frac{1}{2}(\frac{x-\bar{x}}{\sigma})^2}$$

First, we define a function for this function:

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

We can use our function alongside our trapezoid function to get the area and compare it to known values.

x = np.linspace(-1,1,1000)
y = mygauss(x,0,1)
a = trap(x,y)
a_theo = 0.682
error = (a-a_theo)/a_theo*100
print(f"percent error = {error}%")

This returns a %error of 0.101%.

x = np.linspace(-2,2,100)
y = mygauss(x,0,1)
a = trap(x,y)
a_theo = 0.272+0.682
error = (a-a_theo)/a_theo*100
print(f"percent error = {error}%")

This returns a %error of 0.049%.

x = np.linspace(-3,3,100)
y = mygauss(x,0,1)
a = trap(x,y)
a_theo = 0.042+0.272+0.682
error = (a-a_theo)/a_theo*100
print(f"percent error = {error}%")

This returns a %error of 0.130%.

Compared to Built-In Functions

Numpy has its own trapezoid integration function so we can compare our function to theirs.

# import scipy.integrate.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
a=trap(x,y)
a_t=np.trapezoid(y,x)
# print your result and the expected answer in a formatted print statement
print(f"res: {a}, exp: {a_t}")

This returns: res: 6.276436078403621, exp: 6.276436078403619

# import scipy.integrate.trapezoid

# set up x and y arrays for the gaussian
x = np.linspace(-1,1,100)
y = mygauss(x,0,1)
# call scipy.integrate.trapezoid to perform the integration between $-1\sigma \le x \le 1\sigma$.
a=trap(x,y)
a_t=np.trapezoid(y,x)
# print your result and the expected answer in a formatted print statement
print(f"res: {a}, exp: {a_t}")

This returns: res: 0.6826730329991482, exp: 0.6826730329991479.

Monte Carlo

For x and y arrays containing random values within the plot area of our half-circle, we define a function to plot them and calculate the area.

def mcintegrate(x,y,visualize=False):
  if visualize==True:
    plt.figure()
    fx = mycircle(x)

    plt.xlabel("x")
    plt.ylabel("y")

    goodx=[]
    goody=[]
    badx=[]
    bady=[]

    for i in range(0,len(x)):
      if y[i]<=fx[i]:
        goodx.append(x[i])
        goody.append(y[i])
      else:
        badx.append(x[i])
        bady.append(y[i])

    plt.plot(goodx,goody,"bo",label="good")
    plt.plot(badx,bady,"ro",label="bad")



    title=len(goodx)/(len(goodx)+len(badx))*8
    plt.title(f"area={title}")

    x1=np.linspace(-2,2,100)
    y1=mycircle(x1)
    plt.plot(x1,y1,"k-",label="function",linewidth=4)

    plt.legend()
  else:
    return None

Then we run our function:

x=np.random.uniform(-2,2,1000)
y=np.random.uniform(0,2,1000)
mcintegrate(x,y,visualize=True)

Resultant plot:

The title of the plot (area) is calculated using the formula that was displayed above in the Overview.