Skip to content

polykin.math¤

root_newton ¤

root_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)\). It's 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 smaller 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 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/solvers.py 
37
38
39
40
41
42
43
44
45
46
47
48
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
def root_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)$. It's 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, optional
        Absolute tolerance for `x` value. The algorithm will terminate when the
        change in `x` between two iterations is smaller than `xtol`.
    ftol : float, optional
        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 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
    """
    success = False
    niter = 0
    while niter < maxiter:
        dfdx, f0 = derivative_complex(f, x0)
        if (abs(f0) < ftol):
            success = True
            break
        x1 = x0 - f0 / dfdx
        niter += 1
        if (abs(x1 - x0) < xtol):
            success = True
            x0 = x1
            f0 = f(x0).real
            break
        x0 = x1

    return RootResult(success, niter, x0, f0)