scijit.optimize.newton

scijit.optimize.newton(f, x0, fprime=None, args=(), tol=1.48e-08, maxiter=50, fprime2=None, x1=None, rtol=0.0, full_output=False, disp=True)

Newton-Raphson, secant and Halley root-finder.

Which derivatives are supplied selects the method: no fprime runs the secant method, fprime alone runs Newton-Raphson, fprime with fprime2 runs Halley’s method.

Callback style A for every callback: plain @njit functions of one float, passed as first-class arguments:

@njit
def f(x):
    return x * x - 2.0

@njit
def fp(x):
    return 2.0 * x

res = newton(f, 1.0, fp)
res.root, res.converged

Parameters:

f@njit function f(x) -> float

Function whose root is wanted.

x0float or complex

Starting guess. A complex x0 runs the same iteration over complex128 and returns a complex root, on all three methods. An ARRAY of any shape runs a different algorithm over every element at once, described under Notes; the callbacks then receive the whole array and return one of the same shape. A 0-d array is a scalar start.

fprime@njit function fprime(x) -> float, or None, optional

First derivative. None (default) runs the secant method with x1 as the second starting guess. Positional third.

argstuple, optional

Extra arguments, unpacked into every call as f(x, *args) and into fprime and fprime2 the same way. Must be a tuple; anything else raises TypeError. Default ().

tolfloat, optional

Absolute step tolerance; convergence when |p - p0| <= tol + rtol * |p0|. Default 1.48e-8.

maxiterint, optional

Iteration cap. Default 50. maxiter < 1 raises ValueError.

fprime2@njit function fprime2(x) -> float, or None, optional

Second derivative. Supplying it alongside fprime runs Halley’s method. Supplying it without fprime raises ValueError.

x1float, complex or None, optional

Second starting guess for the secant method, used only when fprime is None. The default None seeds it from x0, x0 * (1 + 1e-4) moved a further 1e-4 away from zero. Any float is taken as the seed, NaN included. x1 == x0 raises ValueError.

rtolfloat, optional

Relative step tolerance. Default 0.0.

full_outputbool, optional

False (default) returns the root alone. True returns (x, RootResults). Inside @njit it must be a compile-time constant; see Notes.

dispbool, optional

True (default) raises RuntimeError when the iteration limit is reached. False returns converged=False instead.

Returns:

xfloat or complex

The estimated root, when full_output is False. Complex when x0 is.

(x, res)tuple of (float or complex, RootResults)

When full_output is True. res fields are reached by attribute, by index or by unpacking.

rootfloat or complex

Best estimate.

iterationsint

Iterations used. An immediate exact hit, where f(x0) == 0, reports 0.

function_callsint

Evaluations of f, plus the derivative evaluations, counted on the same counter.

convergedbool

True if a tolerance test was met before maxiter.

flagstr

'converged', or 'convergence error'.

methodstr

The method that produced the result, by name.

resArrayNewtonResult

When x0 is an array of one dimension or more, and full_output is True. root, converged and zero_der are each shaped like x0. Under full_output=False that route returns the root member alone.

Raises:

ValueError

If fprime2 is given without fprime; if tol <= 0; if maxiter < 1; or if x1 == x0.

RuntimeError

If maxiter is reached, or the derivative is zero so the next step is undefined, or two secant iterates coincide. All three are suppressed by disp=False. On an array x0 the iteration limit raises only when EVERY element failed, and disp does not suppress it.

Warns:

RuntimeWarning

Under disp=False, on two exits: "Derivative was zero." and "Tolerance of {p1 - p0} reached.". Plain non-convergence does not warn.

An array x0 warns from four more, whatever disp is: "all derivatives were zero" and "some derivatives were zero" when a derivative was supplied, "RMS of {rms:g} reached" on the secant, and "some failed to converge after {maxiter} iterations" when part of the array reached the limit.

See also

scipy.optimize.newton

The scipy routine this mirrors.

scijit.optimize.root_scalar

The same methods, by name.

scijit.optimize.brentq

Bracketed, and guaranteed to converge.

Notes

full_output selects the RETURN SHAPE, and a compiled function has one return type per signature, so inside @njit the flag has to be readable when the call compiles. A literal, an omitted default and a module-level constant all are; a variable is not, and raises TypingError naming the constraint. From Python a runtime value is fine.

The result is a namedtuple, where scipy’s is a dict subclass. See RootResults.

Which of the three methods runs follows the VALUE of fprime and fprime2 from Python, as it does in scipy: one variable holding either a function or None reaches both methods and reports the one it ran. Inside @njit the choice follows the TYPE instead, so it is settled when the call compiles, and a variable that may hold either does not type. Both spellings work there written on their own.

fprime2 without fprime raises. scipy’s Halley block sits inside if fprime is not None, so scipy ignores a lone fprime2, runs the secant loop and reports method='secant'.

An array x0 runs a different algorithm, not the scalar loop repeated. Every element advances together, an element whose derivative is zero stops being updated while the others continue, and the outcome is reported per element in converged and zero_der. x1, rtol and disp are not read on that route, and the result carries no iterations or function_calls. Its secant seed is eps ** 0.33, where the scalar secant seeds with 1e-4.

A LENGTH-1 array takes the array route here. scipy switches on np.size(x0) > 1, so a length-1 array takes its scalar route instead. Under full_output scipy returns (array, RootResults) for that one size and this returns ArrayNewtonResult. The size of an array is a run-time value while its number of dimensions is not, and a compiled function has one return type per signature. A 0-d array is the scalar route on both sides.

An integer array x0 returns float64. scipy returns the input integer dtype in the single case where no element is ever updated, which is an x0 whose every element is already an exact root.

Newton’s method is not bracketed, so a bad start can diverge. There is no warning on divergence, only the iteration limit.

Pure @njit, safe to call from a numba.prange loop.

https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.newton.html

Examples

>>> from numba import njit
>>> from scijit.optimize import newton
>>> @njit
... def f(x):
...     return x * x - 2.0
>>> @njit
... def fp(x):
...     return 2.0 * x
>>> @njit
... def run():
...     return newton(f, 1.0, fp)
>>> round(run(), 12)
1.414213562373
>>> @njit
... def run_full():
...     return newton(f, 1.0, fp, full_output=True)
>>> x, res = run_full()
>>> round(x, 12), res.converged, res.iterations
(1.414213562373, True, 5)

With no derivative, which runs the secant method:

>>> @njit
... def run_secant():
...     return newton(f, 1.0)
>>> round(run_secant(), 12)
1.414213562373

An array of starting points, solved together. The callbacks receive the whole array:

>>> import numpy as np
>>> @njit
... def run_vec():
...     return newton(f, np.array([1.0, 10.0, -3.0]), fp)
>>> np.round(run_vec(), 12)
array([ 1.41421356,  1.41421356, -1.41421356])
>>> @njit
... def run_vec_full():
...     return newton(f, np.array([1.0, 10.0, -3.0]), fp,
...                   full_output=True)
>>> res = run_vec_full()
>>> res.converged, res.zero_der
(array([ True,  True,  True]), array([False, False, False]))