polykin.math.roots

fzero_brent

fzero_brent(
    f: Callable[[float], float],
    xa: float,
    xb: float,
    xtol: float = 1e-06,
    ftol: float = 1e-06,
    maxiter: int = 50,
) -> RootResult

Find the root of a scalar function using Brent's 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
xtol	Absolute tolerance for x value. The algorithm will terminate when the change in x between two iterations is less or equal than xtol. TYPE: float DEFAULT: 1e-06
ftol	Absolute tolerance for function value. The algorithm will terminate when |f(x)|<=ftol. TYPE: float DEFAULT: 1e-06
maxiter	Maximum number of iterations. TYPE: int DEFAULT: 50

RETURNS	DESCRIPTION
RootResult	Dataclass with root solution results.

Examples:

Find a root of the Flory-Huggins equation.

>>> from polykin.math import fzero_brent
>>> from math 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.py

def fzero_brent(f: Callable[[float], float],
                xa: float,
                xb: float,
                xtol: float = 1e-6,
                ftol: float = 1e-6,
                maxiter: int = 50
                ) -> RootResult:
    r"""Find the root of a scalar function using Brent's 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.
    xtol : float
        Absolute tolerance for `x` value. The algorithm will terminate when the
        change in `x` between two iterations is less or equal than `xtol`.
    ftol : float
        Absolute tolerance for function value. The algorithm will terminate
        when `|f(x)|<=ftol`.
    maxiter : int
        Maximum number of iterations.

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

    Examples
    --------
    Find a root of the Flory-Huggins equation.
    >>> from polykin.math import fzero_brent
    >>> from math 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 = ""

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

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

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

    xc, fc = xb, fb
    success = False

    for niter in range(1, maxiter+1):

        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

        tol1 = 2*eps*abs(xb) + 0.5*xtol
        m = 0.5*(xc - xb)
        if abs(m) <= tol1:
            message = "|Δx| <= xtol"
            success = True
            break

        if abs(fb) <= ftol:
            message = "|f| <= ftol"
            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 += math.copysign(tol1, m)

        fb = f(xb)
        nfeval += 1

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

    return RootResult(success, message, nfeval, niter, xb, fb)