polykin.math.roots

root_newton

root_newton(
    f: Callable[[complex], complex],
    x0: 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 Newton-Raphson method.

The Newton-Raphson method uses the first derivative of the function to iteratively find the root according to the formula:

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

Unlike the equivalent method in scipy, this method uses complex step differentiation to estimate the derivative of \(f(x)\) without loss of precision. Therefore, there is no need to provide \(f'(x)\). Its application is restricted to real functions that can be evaluated with complex inputs, but which per se do not implement complex arithmetic.

Note

The functions from the math module (e.g., math.log) do not support complex arguments. Use the equivalent functions from numpy instead.

PARAMETER	DESCRIPTION
f	Function whose root is to be found. TYPE: Callable[[complex], complex]
x0	Inital 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_newton
>>> 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_newton(f, 0.3)
>>> print(f"x = {sol.x:.3f}")
x = 0.213

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

def root_newton(f: Callable[[complex], complex],
                x0: 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 Newton-Raphson method.

    The Newton-Raphson method uses the first derivative of the function to
    iteratively find the root according to the formula:

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

    Unlike the equivalent method in [scipy](https://docs.scipy.org/doc/scipy/reference/optimize.root_scalar-newton.html),
    this method uses complex step differentiation to estimate the derivative of
    $f(x)$ without loss of precision. Therefore, there is no need to provide
    $f'(x)$. Its application is restricted to real functions that can be
    evaluated with complex inputs, but which per se do not implement complex
    arithmetic.

    !!! note

        The functions from the `math` module (e.g., `math.log`) do not support
        complex arguments. Use the equivalent functions from `numpy` instead.

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

    nfeval = 0
    message = ""
    success = False

    x = x0

    for k in range(maxiter):

        dfdx, fx = derivative_complex(f, x)
        nfeval += 1

        if verbose:
            print(f"Iteration {k+1}: x = {x}, f(x) = {fx}, df/dx = {dfdx}",
                  flush=True)

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

        if abs(dfdx) <= eps:
            message = f"Nearly zero derivative at x={x} (df/dx={dfdx})."
            break

        Δx = - fx / dfdx

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

        if k + 1 < maxiter:
            x += Δx

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

    return RootResult(success, message, nfeval, k+1, x, fx)