Skip to content

polykin.math.roots¤

fzero_newton ¤

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

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

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.

PARAMETER DESCRIPTION
f

Function whose root is to be found.

TYPE: Callable[[complex], complex]

x0

Inital 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_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 = fzero_newton(f, 0.3)
>>> print(f"x= {sol.x:.3f}")
x= 0.213
Source code in src/polykin/math/roots.py 
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
def fzero_newton(f: Callable[[complex], complex],
                 x0: float,
                 xtol: float = 1e-6,
                 ftol: float = 1e-6,
                 maxiter: int = 50
                 ) -> RootResult:
    r"""Find the root of a scalar function using the newton method.

    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.

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

    nfeval = 0
    message = ""

    success = False
    for niter in range(1, maxiter+1):

        dfdx, f0 = derivative_complex(f, x0)
        nfeval += 1

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

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

        x1 = x0 - f0 / dfdx

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

        x0 = x1

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

    return RootResult(success, message, nfeval, niter, x0, f0)