Python/maths/brent_method.py at b6902519f962cf315d62953f4077e5682e324777 · TheAlgorithms/Python · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
from collections.abc import Callable


def brent_method(
    func: Callable[[float], float],
    left: float,
    right: float,
    tol: float = 1e-8,
    max_iter: int = 100,
) -> float:
    """
    Find the root of function func in the interval [left, right] using Brent's Method.

    Brent's Method combines bisection, secant, and inverse quadratic interpolation.


    Parameters
    ----------
    func : Callable[[float], float]
        Function for which to find the root.
    left : float
        Left endpoint of interval.
    right : float
        Right endpoint of interval.
    tol : float
        Tolerance for convergence (default 1e-8).
    max_iter : int
        Maximum number of iterations (default 100).


    Returns
    -------
    float
        Approximate root of func in [left, right].

    Raises
    ------
    ValueError
        If func(left) and func(right) do not have opposite signs.

    Examples
    --------
    >>> def f(x): return x**3 - x - 2
    >>> round(brent_method(f, 1, 2), 5)
    1.52138

    >>> def f2(x): return x**2 + 1
    >>> brent_method(f2, 0, 1)
    Traceback (most recent call last):
    ...
    ValueError: func(left) and func(right) must have opposite signs
    """
    fl = func(left)
    fr = func(right)

    if fl * fr >= 0:
        raise ValueError("func(left) and func(right) must have opposite signs")

    if abs(fl) < abs(fr):
        left, right = right, left
        fl, fr = fr, fl

    c = left
    fc = fl
    d = right - left

    for iteration in range(max_iter):
        if fr == 0:
            return right

        if fc not in (fl, fr):
            # Inverse quadratic interpolation
            s = (
                left * fr * fc / ((fl - fr) * (fl - fc))
                + right * fl * fc / ((fr - fl) * (fr - fc))
                + c * fl * fr / ((fc - fl) * (fc - fr))
            )
        else:
            # Secant method
            s = right - fr * (right - left) / (fr - fl)

        conditions = [
            not ((3 * left + right) / 4 < s < right)
            if right > left
            else not (right < s < (3 * left + right) / 4),
            iteration > 1 and abs(s - right) >= abs(right - c) / 2,
            iteration <= 1 and abs(s - right) >= abs(c - d) / 2,
            iteration > 1 and abs(right - c) < tol,
            iteration <= 1 and abs(c - d) < tol,
        ]

        if any(conditions):
            # Bisection fallback
            s = (left + right) / 2
            d = right - left

        fs = func(s)
        d, c = c, right
        fc = fr

        if fl * fs < 0:
            right = s
            fr = fs
        else:
            left = s
            fl = fs

        if abs(fl) < abs(fr):
            left, right = right, left
            fl, fr = fr, fl

        if abs(right - left) < tol:
            return right

    return right


if __name__ == "__main__":
    import doctest

    doctest.testmod(verbose=True)