polykin.math.roots

root_secant

root_secant(
    f: Callable[[float], float],
    x0: float,
    x1: float,
    *,
    tolx: float = 1e-06,
    tolf: float = 1e-06,
    maxiter: int = 50,
    callback: (
        Callable[[int, float, float], tuple[bool, bool]]
        | None
    ) = None
) -> RootResult

Find the root of a scalar function using the secant method.

The secant method uses two initial guesses and approximates the derivative of the function to iteratively find the root according to the formula:

\[ x_{k+1} = x_k - f(x_k) \frac{x_k - x_{k-1}}{f(x_k) - f(x_{k-1})} \]

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	First initial guess. TYPE: float
x1	Second initial guess. 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
callback	Optional callback with signature callback(niter, x, fx)->(stop, success) called at each iteration. If stop is True, the iteration is terminated. If success is True, the optimization is considered successful. TYPE: Callable[[int, float, float], tuple[bool, bool]] | None DEFAULT: None

RETURNS	DESCRIPTION
RootResult	Dataclass with root solution results.

Examples:

Find a root of the Flory-Huggins equation.

>>> from polykin.math import root_secant
>>> 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 = root_secant(f, 0.3, 0.31)
>>> print(f"x = {sol.x:.3f}")
x = 0.213

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

def root_secant(
    f: Callable[[float], float],
    x0: float,
    x1: float,
    *,
    tolx: float = 1e-6,
    tolf: float = 1e-6,
    maxiter: int = 50,
    callback: Callable[[int, float, float], tuple[bool, bool]] | None = None,
) -> RootResult:
    r"""Find the root of a scalar function using the secant method.

    The secant method uses two initial guesses and approximates the derivative
    of the function to iteratively find the root according to the formula:

    $$ x_{k+1} = x_k - f(x_k) \frac{x_k - x_{k-1}}{f(x_k) - f(x_{k-1})} $$

    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
        First initial guess.
    x1 : float
        Second initial guess.
    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.
    callback : Callable[[int, float, float], tuple[bool, bool]] | None
        Optional callback with signature `callback(niter, x, fx)->(stop, success)` called
        at each iteration. If `stop` is `True`, the iteration is terminated. If `success`
        is `True`, the optimization is considered successful.

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

    Examples
    --------
    Find a root of the Flory-Huggins equation.
    >>> from polykin.math import root_secant
    >>> 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 = root_secant(f, 0.3, 0.31)
    >>> print(f"x = {sol.x:.3f}")
    x = 0.213
    """
    method = "Secant"
    success = False
    message = ""
    nfeval = 0

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

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

    x2, f2 = float("nan"), float("nan")
    niter = 0

    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 x[k-1]={x0} and x[k]={x1} (Δf={Δf})."
            break

        x2 = x1 - f1 * (x1 - x0) / Δf

        f2 = f(x2)
        nfeval += 1

        if callback is not None:
            stop, _success = callback(niter, x2, f2)
            if stop:
                message = "Terminated by user callback."
                success = _success
                break

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

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

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

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

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