lab 7 Numerical Integration - FouadM0426/comp.physics-spring-2025 GitHub Wiki
Lab 7 Numerical Integration
- Goal In this lab we will be using different methods for numerical integration, numerical integration is a way approximate values for integrals that can be difficult to solve by hand. In this lab we will use our mathematical techniques to solve these situations, and compare to the python built-in code
- Part 1 Trapezoid We will write a code to implement trapezoid integral method code below
import numpy as np
import matplotlib.pyplot as plt
import math
# write your function here
def trap(y, x):
'''input y,x
output integral of f(x) between xmin and xmax'''
length = len(x)
integral = 0
for i in range(length - 1):
integral += (x[i + 1] - x[i]) * (y[i] + y[i + 1]) / 2
return integral
- Part 2 half circle In this part we will duplicate the concept above but with a semi-circle '''python
write your function mycircle here
def mycircle(x, r=2): '''input x, r output y value of a circle''' return np.sqrt(r2 - x2)
Now to show the integration using python
```python
array_x = np.linspace(-2, 2, 100)
array_y = mycircle(array_x)
area = trap(array_y, array_x)
print(f"The area under the half circle is equal = {area:.5f}")
The area under the half circle is equal = 6.27644
percent error of answer
accepted_value = np.pi * 2**2 / 2
percent_error = np.abs((area - accepted_value) / accepted_value) * 100
print(f"The percent error is {percent_error:.5f}%")
The percent error is 0.10742%
Modifying the spacing in the x-array
array_x = np.linspace(-2, 2, 1000)
array_y = mycircle(array_x)
area = trap(array_y, array_x)
print(f"The area under the half circle is equal = {area:.5f}")
The area under the half circle is equal = 6.28297
- Part 3 Integral of the Gaussian Distribution The gaussian function, 1/2pi(sigma)e^((x-x_avg)/sigma)^2 we will integrate this function to begin with.
def gauss(x,mean,sigma):
calc = 1/(np.sqrt(2*np.pi)*sigma)*np.exp(-0.5*((x-mean)/sigma)**2)
return calc
x =np.linspace(-1, 1, 100)
mean = 0
sigma = 1
y = gauss(x,mean,sigma)
intergal_mine = trap(y,x)
print(f"The integral of the gaussian function is equal to {intergal_mine:.5f}")
import scipy.integrate
integral = scipy.integrate.trapezoid(y,x)
print(f"The integral of the gaussian function is equal to {integral:.5f}")
print(math.isclose(.99729,1.2,abs_tol = (1.2-.99729)))
The integral of the gaussian function is equal to 0.68267
The integral of the gaussian function is equal to 0.68267
True
x =np.linspace(-2, 2, 100)
mean = 0
sigma = 1
y = gauss(x,mean,sigma)
integral_mine = trap(y,x)
print(f"The integral of the gaussian function is equal to {integral_mine:.5f}")
import scipy.integrate
integral = scipy.integrate.trapezoid(y,x)
print(f"The integral of the gaussian function is equal to {integral:.5f}")
print(math.isclose(.99729,1.2,abs_tol = (1.2-.99729)))
The integral of the gaussian function is equal to 0.95447
The integral of the gaussian function is equal to 0.95447
True
Now to do the integral in range −3σ≤x≤3σ
x =np.linspace(-3, 3, 100)
mean = 0
sigma = 1
y = gauss(x,mean,sigma)
integral_mine = trap(y,x)
print(f"The integral of the gaussian function is equal to {integral_mine:.5f}")
import scipy.integrate
integral = scipy.integrate.trapezoid(y,x)
print(f"The integral of the gaussian function is equal to {integral:.5f}")
print(math.isclose(.99729,1.2,abs_tol = (1.2-.99729)))
The integral of the gaussian function is equal to 0.99729
The integral of the gaussian function is equal to 0.99729
True
As you can see our answers were close but since we used the math.isclose function to compare it to the .99729 it says true meaning that were within approx of the actual.
-Part 4 built in functions In this part we will be using functions that were already created on part 1 questions of this lab
# 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
integral = scipy.integrate.trapezoid(y, x)
# print your result and the expected answer in a formatted print statement
print(f"The integral of the half circle is equal to {integral:.5f}")
#The integral of the half circle is equal to 6.27644
# import scipy.integrate.trapezoid
# set up x and y arrays for the gaussian
x = np.linspace(-1, 1, 100)
mean = 0
sigma = 1
y = gauss(x,mean,sigma)
# call scipy.integrate.trapezoid to perform the integration between $-1\sigma \le x \le 1\sigma$.
intergal = scipy.integrate.trapezoid(y, x)
# print your result and the expected answer in a formatted print statement
print(f"The integral of the gaussian function is equal to {integral:.5f}")
#The integral of the gaussian function is equal to 6.27644
Part 5 monte Carlo
we will use the function created above but to visualize what the half circle would look like in my example I used 250,000 points to plot trying to get as close as possible to the actual value
# your code to calculate the area of a half circle using monte carlo
z = mcintegrate(2,visualize=True)