simpson rule Algorithm
The Simpson's Rule algorithm is a numerical integration technique used to approximate definite integrals of functions. Named after the British mathematician Thomas Simpson, this method is part of a broader family of techniques called Newton-Cotes formulas, which are used for approximating the definite integral of a function using equally spaced sample points. Simpson's Rule specifically employs a quadratic interpolating function to approximate the area under the curve of a given function, providing a more accurate estimate than other simpler methods, such as the midpoint or trapezoidal rule.
Simpson's Rule works by dividing the area under the curve into an even number of equally spaced intervals, then fitting a parabola (a quadratic polynomial) to each pair of adjacent intervals. The algorithm then computes the area under each parabola and sums these areas to obtain an approximation of the total area under the curve, which represents the definite integral of the given function. This method is particularly effective for functions that have a smooth curve without abrupt changes, and its accuracy generally increases with the number of intervals used. However, it is important to note that, like other numerical integration methods, Simpson's Rule is still an approximation and may not always provide an exact solution, especially for functions with complex or irregular behavior.
"""
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 2:
"Simpson Rule"
"""
def method_2(boundary, steps):
# "Simpson Rule"
# int(f) = delta_x/2 * (b-a)/3*(f1 + 4f2 + 2f_3 + ... + 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 / 3.0) * f(a)
cnt = 2
for i in x_i:
y += (h / 3) * (4 - 2 * (cnt % 2)) * f(i)
cnt += 1
y += (h / 3.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_2(boundary, steps)
print(f"y = {y}")
if __name__ == "__main__":
main()
