polykin.math

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/solvers.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
    """

    fa = f(xa)
    if (abs(fa) < ftol):
        return RootResult(True, 0, xa, fa)
    fb = f(xb)
    if (abs(fb) < ftol):
        return RootResult(True, 0, xb, fb)

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

    xc, fc = xb, fb
    success = False
    for iter 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
        tol1 = 2*eps*abs(xb) + 0.5*xtol
        xm = 0.5*(xc - xb)
        if (abs(xm) <= tol1) or (abs(fb) <= ftol):
            # return xb
            success = True
            break
        if (abs(e) >= tol1) and (abs(fa) > abs(fb)):
            s = fb/fa
            if xa == xc:
                p = 2*xm*s
                q = 1 - s
            else:
                q = fa/fc
                r = fb/fc
                p = s*(2*xm*q*(q - r) - (xb - xa)*(r - 1))
                q = (q - 1)*(r - 1)*(s - 1)
            if p > 0:
                q = -q
            p = abs(p)
            min1 = 3*xm*q - abs(tol1*q)
            min2 = abs(e*q)
            if 2*p < min(min1, min2):
                e = d
                d = p/q
            else:
                d = xm
                e = d
        else:
            d = xm
            e = d
        xa, fa = xb, fb
        if abs(d) > tol1:
            xb += d
        else:
            xb += math.copysign(tol1, xm)
        fb = f(xb)

    return RootResult(success, iter+1, xb, fb)