Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

chore: improve comments and add tests to trapezoidal rule #11640

Merged
merged 12 commits into from
Dec 31, 2024
80 changes: 68 additions & 12 deletions maths/trapezoidal_rule.py
Original file line number Diff line number Diff line change
@@ -1,47 +1,103 @@
"""
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'
This method is the classical approach of summing 'Equally Spaced Abscissas'

method 1:
"extended trapezoidal rule"
Method 1:
"Extended Trapezoidal Rule"

"""


def method_1(boundary, steps):
# "extended trapezoidal rule"
# int(f) = dx/2 * (f1 + 2f2 + ... + fn)
"""
This function implements the extended trapezoidal rule for numerical integration.
The function f(x) is provided below.

:param boundary: List containing the lower and upper bounds of integration [a, b]
:param steps: The number of steps (intervals) used in the approximation
:return: The numerical approximation of the integral

>>> abs(method_1([0, 1], 10) - 0.33333) < 0.01
True

>>> abs(method_1([0, 1], 100) - 0.33333) < 0.01
True

>>> abs(method_1([0, 2], 1000) - 2.66667) < 0.01
True

>>> abs(method_1([1, 2], 1000) - 2.33333) < 0.01
True
"""
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):
"""
Generates the points between a and b with spacing h for trapezoidal integration.

:param a: The lower bound of integration
:param b: The upper bound of integration
:param h: The step size
:yield: The next x-value in the range (a, b)

>>> list(make_points(0, 1, 0.1))
[0.1, 0.2, 0.30000000000000004, \
0.4, 0.5, 0.6, 0.7, \
0.7999999999999999, \
0.8999999999999999]
"""
x = a + h
while x < (b - h):
yield x
x = x + h


def f(x): # enter your function here
y = (x - 0) * (x - 0)
def f(x):
"""
This is the function to integrate, f(x) = (x - 0)^2 = x^2.

:param x: The input value
:return: The value of f(x)

>>> f(0)
0.0

>>> f(1)
1.0

>>> f(0.5)
0.25
"""
y = float((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
"""
Main function to test the trapezoidal rule.
:a: Lower bound of integration
:b: Upper bound of integration
:steps: define number of steps or resolution
:boundary: define boundary of integration

>>> main()
y = 0.3349999999999999
"""
a = 0.0
b = 1.0
steps = 10.0
boundary = [a, b]
y = method_1(boundary, steps)
print(f"y = {y}")

Expand Down