From 2e67da3a7b0f009f0c41c0f350a2a9207d3b5bba Mon Sep 17 00:00:00 2001 From: Julien RICHARD Date: Tue, 1 Oct 2024 19:28:17 +0200 Subject: [PATCH] chore: improve comments and add tests to trapezoidal rule --- maths/trapezoidal_rule.py | 77 +++++++++++++++++++++++++++++++++------ 1 file changed, 65 insertions(+), 12 deletions(-) diff --git a/maths/trapezoidal_rule.py b/maths/trapezoidal_rule.py index 9a4ddc8af..7ac0771f1 100644 --- a/maths/trapezoidal_rule.py +++ b/maths/trapezoidal_rule.py @@ -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}")