polykin.math.roots

fzero_secant

fzero_secant(
    f: Callable[[float], float],
    x0: float,
    x1: float,
    xtol: float = 1e-06,
    ftol: float = 1e-06,
    maxiter: int = 50,
) -> RootResult

Find the root of a scalar function using 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.

PARAMETER	DESCRIPTION
f	Function whose root is to be found. TYPE: Callable[[float], float]
x0	Initial guess. TYPE: float
x1	Second guess. 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_secant
>>> 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_secant(f, 0.3, 0.31)
>>> print(f"x= {sol.x:.3f}")
x= 0.213

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

def fzero_secant(f: Callable[[float], float],
                 x0: float,
                 x1: float,
                 xtol: float = 1e-6,
                 ftol: float = 1e-6,
                 maxiter: int = 50
                 ) -> RootResult:
    r"""Find the root of a scalar function using the secant method.

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

    Parameters
    ----------
    f : Callable[[float], float]
        Function whose root is to be found.
    x0 : float
        Initial guess.
    x1 : float
        Second guess.
    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_secant
    >>> 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_secant(f, 0.3, 0.31)
    >>> print(f"x= {sol.x:.3f}")
    x= 0.213
    """

    nfeval = 0
    message = ""

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

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

    success = False
    x2, f2 = np.nan, np.nan

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

        Δf = f1 - f0
        if abs(Δf) <= eps * max(abs(f0), abs(f1), 1.0):
            message = f"Nearly zero slope between x0={x0} and x1={x1} (Δf={Δf})."
            break

        x2 = x1 - f1*(x1 - x0) / Δf
        f2 = f(x2)
        nfeval += 1

        if (abs(x2 - x1) <= xtol):
            message = "|Δx| <= xtol"
            success = True
            break

        if (abs(f2) <= ftol):
            message = "|f| <= ftol"
            success = True
            break

        x0, f0 = x1, f1
        x1, f1 = x2, f2

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

    return RootResult(success, message, nfeval, niter, x2, f2)