polykin.math.roots

fzero_brent

fzero_brent(
    f: Callable[[float], float],
    xa: float,
    xb: float,
    tolx: float = 1e-06,
    tolf: float = 1e-06,
    maxiter: int = 50,
    verbose: bool = False,
) -> RootResult

Find the root of a scalar function using Brent's method.

Brent's method is a root-finding algorithm combining bisection, secant, and inverse quadratic interpolation. It is robust like the bisection method and fast like the secant method.

Unlike the equivalent method in scipy, this method also terminates based on the function value. This is sometimes a more meaningful stop criterion.

References

• William H. Press, Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery. 2007. "Numerical Recipes 3rd Edition: The Art of Scientific Computing" (3rd. ed.). Cambridge University Press, USA.

PARAMETER	DESCRIPTION
f	Function whose root is to be found. TYPE: Callable[[float], float]
xa	Lower bound of the bracketing interval. TYPE: float
xb	Upper bound of the bracketing interval. TYPE: float
tolx	Absolute tolerance for x value. The algorithm will terminate when the change in x between two iterations is less or equal than tolx. TYPE: float DEFAULT: 1e-06
tolf	Absolute tolerance for function value. The algorithm will terminate when |f(x)|<=tolf. TYPE: float DEFAULT: 1e-06
maxiter	Maximum number of iterations. TYPE: int DEFAULT: 50
verbose	Print iteration information. TYPE: bool DEFAULT: False

RETURNS	DESCRIPTION
RootResult	Dataclass with root solution results.

Examples:

Find a root of the Flory-Huggins equation.

>>> from polykin.math import fzero_brent
>>> from numpy import log
>>> def f(x, a=0.6, chi=0.4):
...     return log(x) + (1 - x) + chi*(1 - x)**2 - log(a)
>>> sol = fzero_brent(f, 0.1, 0.9)
>>> print(f"x = {sol.x:.3f}")
x = 0.213

Source code in src/polykin/math/roots/scalar.py

def fzero_brent(f: Callable[[float], float],
                xa: float,
                xb: float,
                tolx: float = 1e-6,
                tolf: float = 1e-6,
                maxiter: int = 50,
                verbose: bool = False
                ) -> RootResult:
    r"""Find the root of a scalar function using Brent's method.

    Brent's method is a root-finding algorithm combining bisection, secant,
    and inverse quadratic interpolation. It is robust like the bisection method
    and fast like the secant method.

    Unlike the equivalent method in [scipy](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.brentq.html),
    this method also terminates based on the function value. This is sometimes
    a more meaningful stop criterion.

    **References**

    * William H. Press, Saul A. Teukolsky, William T. Vetterling, and 
      Brian P. Flannery. 2007. "Numerical Recipes 3rd Edition: The Art of 
      Scientific Computing" (3rd. ed.). Cambridge University Press, USA.

    Parameters
    ----------
    f : Callable[[float], float]
        Function whose root is to be found.
    xa : float
        Lower bound of the bracketing interval.
    xb : float
        Upper bound of the bracketing interval.
    tolx : float
        Absolute tolerance for `x` value. The algorithm will terminate when the
        change in `x` between two iterations is less or equal than `tolx`.
    tolf : float
        Absolute tolerance for function value. The algorithm will terminate
        when `|f(x)|<=tolf`.
    maxiter : int
        Maximum number of iterations.
    verbose : bool
        Print iteration information.

    Returns
    -------
    RootResult
        Dataclass with root solution results.

    Examples
    --------
    Find a root of the Flory-Huggins equation.
    >>> from polykin.math import fzero_brent
    >>> from numpy import log
    >>> def f(x, a=0.6, chi=0.4):
    ...     return log(x) + (1 - x) + chi*(1 - x)**2 - log(a)
    >>> sol = fzero_brent(f, 0.1, 0.9)
    >>> print(f"x = {sol.x:.3f}")
    x = 0.213
    """

    nfeval = 0
    message = ""
    success = False

    fa = f(xa)
    nfeval += 1
    if abs(fa) <= tolf:
        message = "|f(xa)| ≤ tolf"
        return RootResult(True, message, nfeval, 0, xa, fa)

    fb = f(xb)
    nfeval += 1
    if abs(fb) <= tolf:
        message = "|f(xb)| ≤ tolf"
        return RootResult(True, message, nfeval, 0, xb, fb)

    if (fa*fb) > 0.0:
        raise ValueError("Root is not bracketed.")

    xc, fc = xb, fb

    for k in range(maxiter):

        if (fb*fc > 0.0):
            xc, fc = xa, fa
            d = xb - xa
            e = d

        if abs(fc) < abs(fb):
            xa, fa = xb, fb
            xb, fb = xc, fc
            xc, fc = xa, fa

        if verbose:
            print(f"Iteration {k+1}: x = {xb}, f(x) = {fb}", flush=True)

        tol1 = 2*eps*abs(xb) + 0.5*tolx
        m = 0.5*(xc - xb)

        if abs(m) <= tol1:
            message = "|Δx| ≤ tolx"
            success = True
            break

        if abs(fb) <= tolf:
            message = "|f(x)| ≤ tolf"
            success = True
            break

        if (abs(e) >= tol1) and (abs(fa) > abs(fb)):
            s = fb/fa
            if xa == xc:
                p = 2*m*s
                q = 1 - s
            else:
                q = fa/fc
                r = fb/fc
                p = s*(2*m*q*(q - r) - (xb - xa)*(r - 1))
                q = (q - 1)*(r - 1)*(s - 1)
            if p > 0:
                q = -q
            p = abs(p)
            min1 = 3*m*q - abs(tol1*q)
            min2 = abs(e*q)
            if 2*p < min(min1, min2):
                e = d
                d = p/q
            else:
                d = m
                e = d
        else:
            d = m
            e = d

        xa, fa = xb, fb
        if abs(d) > tol1:
            xb += d
        else:
            xb += np.copysign(tol1, m)

        fb = f(xb)
        nfeval += 1

    else:
        message = f"Maximum number of iterations ({maxiter}) reached."

    return RootResult(success, message, nfeval, k+1, xb, fb)