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,
    verbose: bool = False,
) -> 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
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 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,
                verbose: bool = False
                ) -> 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.
    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 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
    """

    nfeval = 0
    message = ""
    success = False

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

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

    x2, f2 = np.nan, np.nan

    for k in range(maxiter):

        Δ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 verbose:
            print(f"Iteration {k+1}: x = {x2}, f(x) = {f2}", flush=True)

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

        if (abs(f2) <= tolf):
            message = "|f(x)| ≤ tolf"
            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, k+1, x2, f2)