brent7

Author: Dibbu-cell

maths/numerical_analysis/brents_method.py

@@ -5,17 +5,18 @@
 
 Brent's method combines bisection, secant, and inverse quadratic interpolation to efficiently and robustly find a root of a continuous function. It is guaranteed to converge as long as the root is bracketed.
 
-See: https://en.wikipedia.org/wiki/Brent%27s_method
+See:
+    https://en.wikipedia.org/wiki/Brent%27s_method
 
-Author: [Aryan Singh (2nd year LNMIIT)]
+Author: Aryan Singh (2nd year LNMIIT)
 """
 
 from collections.abc import Callable
 
 def brent_method(
     function: Callable[[float], float],
-    lower: float,
-    upper: float,
+    lower_bound: float,
+    upper_bound: float,
     tolerance: float = 1e-14,
     max_iterations: int = 100,
 ) -> float:
@@ -24,16 +25,16 @@ def brent_method(
 
     Args:
         function: A continuous function for which the root is sought.
-        lower: One end of the bracketing interval.
-        upper: The other end of the bracketing interval.
+        lower_bound: One end of the bracketing interval.
+        upper_bound: The other end of the bracketing interval.
         tolerance: The tolerance for convergence (default 1e-14).
         max_iterations: Maximum number of iterations (default 100).
 
     Returns:
         A float value approximating the root.
 
     Raises:
-        ValueError: If the root is not bracketed in [lower, upper].
+        ValueError: If the root is not bracketed in [lower_bound, upper_bound].
 
     Examples:
         >>> brent_method(lambda x: x**3 - 1, -5, 5)
@@ -47,66 +48,73 @@ def brent_method(
             ...
         ValueError: Root is not bracketed in the interval [4, 1000].
     """
-    fa = function(lower)
-    fb = function(upper)
-    if fa * fb >= 0:
-        raise ValueError(f"Root is not bracketed in the interval [{lower}, {upper}].")
-
-    if abs(fa) < abs(fb):
-        lower, upper = upper, lower
-        fa, fb = fb, fa
-
-    c = lower
-    fc = fa
-    d = upper - lower
-    mflag = True
+    function_lower = function(lower_bound)
+    function_upper = function(upper_bound)
+    if function_lower * function_upper >= 0:
+        error_message = (
+            "Root is not bracketed in the interval "
+            f"[{lower_bound}, {upper_bound}]."
+        )
+        raise ValueError(error_message)
+
+    if abs(function_lower) < abs(function_upper):
+        lower_bound, upper_bound = upper_bound, lower_bound
+        function_lower, function_upper = function_upper, function_lower
+
+    previous_bound = lower_bound
+    function_previous = function_lower
+    previous_step = upper_bound - lower_bound
+    bisect_flag = True
 
     for _ in range(max_iterations):
-        if fb == 0:
-            return upper
-        if fc not in {fa, fb}:
+        if function_upper == 0:
+            return upper_bound
+        if function_previous not in {function_lower, function_upper}:
             # Inverse quadratic interpolation
             s = (
-                lower * fb * fc / ((fa - fb) * (fa - fc))
-                + upper * fa * fc / ((fb - fa) * (fb - fc))
-                + c * fa * fb / ((fc - fa) * (fc - fb))
+                lower_bound * function_upper * function_previous
+                / ((function_lower - function_upper) * (function_lower - function_previous))
+                + upper_bound * function_lower * function_previous
+                / ((function_upper - function_lower) * (function_upper - function_previous))
+                + previous_bound * function_lower * function_upper
+                / ((function_previous - function_lower) * (function_previous - function_upper))
             )
         else:
             # Secant method
-            s = upper - fb * (upper - lower) / (fb - fa)
+            s = upper_bound - function_upper * (upper_bound - lower_bound) / (function_upper - function_lower)
 
         conditions = [
-            not ((3 * lower + upper) / 4 < s < upper if upper > lower else upper < s < (3 * lower + upper) / 4),
-            mflag and abs(s - upper) >= abs(upper - c) / 2,
-            not mflag and abs(s - upper) >= abs(c - d) / 2,
-            mflag and abs(upper - c) < tolerance,
-            not mflag and abs(c - d) < tolerance,
+            not ((3 * lower_bound + upper_bound) / 4 < s < upper_bound if upper_bound > lower_bound else upper_bound < s < (3 * lower_bound + upper_bound) / 4),
+            bisect_flag and abs(s - upper_bound) >= abs(upper_bound - previous_bound) / 2,
+            not bisect_flag and abs(s - upper_bound) >= abs(previous_bound - previous_step) / 2,
+            bisect_flag and abs(upper_bound - previous_bound) < tolerance,
+            not bisect_flag and abs(previous_bound - previous_step) < tolerance,
         ]
         if any(conditions):
-            s = (lower + upper) / 2
-            mflag = True
+            s = (lower_bound + upper_bound) / 2
+            bisect_flag = True
         else:
-            mflag = False
+            bisect_flag = False
 
-        fs = function(s)
-        d, c = c, upper
-        fc = fb
+        function_s = function(s)
+        previous_step, previous_bound = previous_bound, upper_bound
+        function_previous = function_upper
 
-        if fa * fs < 0:
-            upper = s
-            fb = fs
+        if function_lower * function_s < 0:
+            upper_bound = s
+            function_upper = function_s
         else:
-            lower = s
-            fa = fs
+            lower_bound = s
+            function_lower = function_s
 
-        if abs(fa) < abs(fb):
-            lower, upper = upper, lower
-            fa, fb = fb, fa
+        if abs(function_lower) < abs(function_upper):
+            lower_bound, upper_bound = upper_bound, lower_bound
+            function_lower, function_upper = function_upper, function_lower
 
-        if abs(upper - lower) < tolerance or fb == 0:
-            return upper
+        if abs(upper_bound - lower_bound) < tolerance or function_upper == 0:
+            return upper_bound
 
-    return upper
+    return upper_bound
 
 
 if __name__ == "__main__":