Skip to content

Commit 16e3b40

Browse files
debesh00800debeshpre-commit-ci[bot]cclauss
authored
Add Brent's Method for root finding (numerical analysis) (#13063)
* Add Brent's Method for root finding (numerical analysis) * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Fix parameter names and lint issues in Brent's Method * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Fix import name * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Rename maths/brent_method.py to maths/numerical_analysis/brent_method.py --------- Co-authored-by: debesh <debeshmaheshwari008@gmail.com> Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Christian Clauss <cclauss@me.com>
1 parent 80c39d0 commit 16e3b40

1 file changed

Lines changed: 121 additions & 0 deletions

File tree

Lines changed: 121 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,121 @@
1+
from collections.abc import Callable
2+
3+
4+
def brent_method(
5+
func: Callable[[float], float],
6+
left: float,
7+
right: float,
8+
tol: float = 1e-8,
9+
max_iter: int = 100,
10+
) -> float:
11+
"""
12+
Find the root of function func in the interval [left, right] using Brent's Method.
13+
14+
Brent's Method combines bisection, secant, and inverse quadratic interpolation.
15+
16+
17+
Parameters
18+
----------
19+
func : Callable[[float], float]
20+
Function for which to find the root.
21+
left : float
22+
Left endpoint of interval.
23+
right : float
24+
Right endpoint of interval.
25+
tol : float
26+
Tolerance for convergence (default 1e-8).
27+
max_iter : int
28+
Maximum number of iterations (default 100).
29+
30+
31+
Returns
32+
-------
33+
float
34+
Approximate root of func in [left, right].
35+
36+
Raises
37+
------
38+
ValueError
39+
If func(left) and func(right) do not have opposite signs.
40+
41+
Examples
42+
--------
43+
>>> def f(x): return x**3 - x - 2
44+
>>> round(brent_method(f, 1, 2), 5)
45+
1.52138
46+
47+
>>> def f2(x): return x**2 + 1
48+
>>> brent_method(f2, 0, 1)
49+
Traceback (most recent call last):
50+
...
51+
ValueError: func(left) and func(right) must have opposite signs
52+
"""
53+
fl = func(left)
54+
fr = func(right)
55+
56+
if fl * fr >= 0:
57+
raise ValueError("func(left) and func(right) must have opposite signs")
58+
59+
if abs(fl) < abs(fr):
60+
left, right = right, left
61+
fl, fr = fr, fl
62+
63+
c = left
64+
fc = fl
65+
d = right - left
66+
67+
for iteration in range(max_iter):
68+
if fr == 0:
69+
return right
70+
71+
if fc not in (fl, fr):
72+
# Inverse quadratic interpolation
73+
s = (
74+
left * fr * fc / ((fl - fr) * (fl - fc))
75+
+ right * fl * fc / ((fr - fl) * (fr - fc))
76+
+ c * fl * fr / ((fc - fl) * (fc - fr))
77+
)
78+
else:
79+
# Secant method
80+
s = right - fr * (right - left) / (fr - fl)
81+
82+
conditions = [
83+
not ((3 * left + right) / 4 < s < right)
84+
if right > left
85+
else not (right < s < (3 * left + right) / 4),
86+
iteration > 1 and abs(s - right) >= abs(right - c) / 2,
87+
iteration <= 1 and abs(s - right) >= abs(c - d) / 2,
88+
iteration > 1 and abs(right - c) < tol,
89+
iteration <= 1 and abs(c - d) < tol,
90+
]
91+
92+
if any(conditions):
93+
# Bisection fallback
94+
s = (left + right) / 2
95+
d = right - left
96+
97+
fs = func(s)
98+
d, c = c, right
99+
fc = fr
100+
101+
if fl * fs < 0:
102+
right = s
103+
fr = fs
104+
else:
105+
left = s
106+
fl = fs
107+
108+
if abs(fl) < abs(fr):
109+
left, right = right, left
110+
fl, fr = fr, fl
111+
112+
if abs(right - left) < tol:
113+
return right
114+
115+
return right
116+
117+
118+
if __name__ == "__main__":
119+
import doctest
120+
121+
doctest.testmod(verbose=True)

0 commit comments

Comments
 (0)