trapezoidal rule Algorithm
The trapezoidal rule algorithm is a numerical integration technique used to estimate the definite integral of a function over a given interval. It approximates the area under the curve by dividing the interval into smaller subintervals and approximating the function with a linear function (trapezoid) on each subinterval. The sum of the areas of all these trapezoids gives an approximation of the integral. The algorithm is widely used in mathematical applications, including engineering, physics, and economics, as a simple and efficient method for solving integrals that do not have an exact analytical solution.
The trapezoidal rule algorithm starts by dividing the interval of integration into n equal subintervals, each with width h = (b - a) / n, where a and b are the limits of integration. On each subinterval, the algorithm calculates the area of the trapezoid formed by connecting the endpoints of the function with a straight line. The area of each trapezoid is given by (h / 2) * (f(x_i) + f(x_{i+1})), where x_i and x_{i+1} are the endpoints of the subinterval, and f(x) is the function being integrated. The algorithm then sums the areas of all the trapezoids to obtain the estimated integral value. The trapezoidal rule provides better accuracy for smooth functions, and its accuracy can be increased by increasing the number of subintervals (n) or by using more advanced integration techniques like Simpson's rule or adaptive quadrature algorithms.
"""
Numerical integration or quadrature for a smooth function f with known values at x_i
This method is the classical approach of suming 'Equally Spaced Abscissas'
method 1:
"extended trapezoidal rule"
"""
def method_1(boundary, steps):
# "extended trapezoidal rule"
# int(f) = dx/2 * (f1 + 2f2 + ... + fn)
h = (boundary[1] - boundary[0]) / steps
a = boundary[0]
b = boundary[1]
x_i = make_points(a, b, h)
y = 0.0
y += (h / 2.0) * f(a)
for i in x_i:
# print(i)
y += h * f(i)
y += (h / 2.0) * f(b)
return y
def make_points(a, b, h):
x = a + h
while x < (b - h):
yield x
x = x + h
def f(x): # enter your function here
y = (x - 0) * (x - 0)
return y
def main():
a = 0.0 # Lower bound of integration
b = 1.0 # Upper bound of integration
steps = 10.0 # define number of steps or resolution
boundary = [a, b] # define boundary of integration
y = method_1(boundary, steps)
print(f"y = {y}")
if __name__ == "__main__":
main()
