chore: improve comments and add tests to trapezoidal rule

This commit is contained in:
Julien RICHARD 2024-10-01 19:28:17 +02:00
parent 0177ae1cd5
commit 2e67da3a7b
No known key found for this signature in database

View File

@ -1,17 +1,35 @@
"""
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]
@ -19,29 +37,64 @@ def method_1(boundary, steps):
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}")