Coverage for /usr/lib/python3/dist-packages/scipy/optimize/_zeros_py.py: 11%
479 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1import warnings
2from collections import namedtuple
3import operator
4from . import _zeros
5from ._optimize import OptimizeResult
6import numpy as np
9_iter = 100
10_xtol = 2e-12
11_rtol = 4 * np.finfo(float).eps
13__all__ = ['newton', 'bisect', 'ridder', 'brentq', 'brenth', 'toms748',
14 'RootResults']
16# Must agree with CONVERGED, SIGNERR, CONVERR, ... in zeros.h
17_ECONVERGED = 0
18_ESIGNERR = -1
19_ECONVERR = -2
20_EVALUEERR = -3
21_EINPROGRESS = 1
23CONVERGED = 'converged'
24SIGNERR = 'sign error'
25CONVERR = 'convergence error'
26VALUEERR = 'value error'
27INPROGRESS = 'No error'
30flag_map = {_ECONVERGED: CONVERGED, _ESIGNERR: SIGNERR, _ECONVERR: CONVERR,
31 _EVALUEERR: VALUEERR, _EINPROGRESS: INPROGRESS}
34class RootResults(OptimizeResult):
35 """Represents the root finding result.
37 Attributes
38 ----------
39 root : float
40 Estimated root location.
41 iterations : int
42 Number of iterations needed to find the root.
43 function_calls : int
44 Number of times the function was called.
45 converged : bool
46 True if the routine converged.
47 flag : str
48 Description of the cause of termination.
50 """
52 def __init__(self, root, iterations, function_calls, flag):
53 self.root = root
54 self.iterations = iterations
55 self.function_calls = function_calls
56 self.converged = flag == _ECONVERGED
57 if flag in flag_map:
58 self.flag = flag_map[flag]
59 else:
60 self.flag = flag
63def results_c(full_output, r):
64 if full_output:
65 x, funcalls, iterations, flag = r
66 results = RootResults(root=x,
67 iterations=iterations,
68 function_calls=funcalls,
69 flag=flag)
70 return x, results
71 else:
72 return r
75def _results_select(full_output, r):
76 """Select from a tuple of (root, funccalls, iterations, flag)"""
77 x, funcalls, iterations, flag = r
78 if full_output:
79 results = RootResults(root=x,
80 iterations=iterations,
81 function_calls=funcalls,
82 flag=flag)
83 return x, results
84 return x
87def _wrap_nan_raise(f):
89 def f_raise(x, *args):
90 fx = f(x, *args)
91 f_raise._function_calls += 1
92 if np.isnan(fx):
93 msg = (f'The function value at x={x} is NaN; '
94 'solver cannot continue.')
95 err = ValueError(msg)
96 err._x = x
97 err._function_calls = f_raise._function_calls
98 raise err
99 return fx
101 f_raise._function_calls = 0
102 return f_raise
105def newton(func, x0, fprime=None, args=(), tol=1.48e-8, maxiter=50,
106 fprime2=None, x1=None, rtol=0.0,
107 full_output=False, disp=True):
108 """
109 Find a root of a real or complex function using the Newton-Raphson
110 (or secant or Halley's) method.
112 Find a root of the scalar-valued function `func` given a nearby scalar
113 starting point `x0`.
114 The Newton-Raphson method is used if the derivative `fprime` of `func`
115 is provided, otherwise the secant method is used. If the second order
116 derivative `fprime2` of `func` is also provided, then Halley's method is
117 used.
119 If `x0` is a sequence with more than one item, `newton` returns an array:
120 the roots of the function from each (scalar) starting point in `x0`.
121 In this case, `func` must be vectorized to return a sequence or array of
122 the same shape as its first argument. If `fprime` (`fprime2`) is given,
123 then its return must also have the same shape: each element is the first
124 (second) derivative of `func` with respect to its only variable evaluated
125 at each element of its first argument.
127 `newton` is for finding roots of a scalar-valued functions of a single
128 variable. For problems involving several variables, see `root`.
130 Parameters
131 ----------
132 func : callable
133 The function whose root is wanted. It must be a function of a
134 single variable of the form ``f(x,a,b,c...)``, where ``a,b,c...``
135 are extra arguments that can be passed in the `args` parameter.
136 x0 : float, sequence, or ndarray
137 An initial estimate of the root that should be somewhere near the
138 actual root. If not scalar, then `func` must be vectorized and return
139 a sequence or array of the same shape as its first argument.
140 fprime : callable, optional
141 The derivative of the function when available and convenient. If it
142 is None (default), then the secant method is used.
143 args : tuple, optional
144 Extra arguments to be used in the function call.
145 tol : float, optional
146 The allowable error of the root's value. If `func` is complex-valued,
147 a larger `tol` is recommended as both the real and imaginary parts
148 of `x` contribute to ``|x - x0|``.
149 maxiter : int, optional
150 Maximum number of iterations.
151 fprime2 : callable, optional
152 The second order derivative of the function when available and
153 convenient. If it is None (default), then the normal Newton-Raphson
154 or the secant method is used. If it is not None, then Halley's method
155 is used.
156 x1 : float, optional
157 Another estimate of the root that should be somewhere near the
158 actual root. Used if `fprime` is not provided.
159 rtol : float, optional
160 Tolerance (relative) for termination.
161 full_output : bool, optional
162 If `full_output` is False (default), the root is returned.
163 If True and `x0` is scalar, the return value is ``(x, r)``, where ``x``
164 is the root and ``r`` is a `RootResults` object.
165 If True and `x0` is non-scalar, the return value is ``(x, converged,
166 zero_der)`` (see Returns section for details).
167 disp : bool, optional
168 If True, raise a RuntimeError if the algorithm didn't converge, with
169 the error message containing the number of iterations and current
170 function value. Otherwise, the convergence status is recorded in a
171 `RootResults` return object.
172 Ignored if `x0` is not scalar.
173 *Note: this has little to do with displaying, however,
174 the `disp` keyword cannot be renamed for backwards compatibility.*
176 Returns
177 -------
178 root : float, sequence, or ndarray
179 Estimated location where function is zero.
180 r : `RootResults`, optional
181 Present if ``full_output=True`` and `x0` is scalar.
182 Object containing information about the convergence. In particular,
183 ``r.converged`` is True if the routine converged.
184 converged : ndarray of bool, optional
185 Present if ``full_output=True`` and `x0` is non-scalar.
186 For vector functions, indicates which elements converged successfully.
187 zero_der : ndarray of bool, optional
188 Present if ``full_output=True`` and `x0` is non-scalar.
189 For vector functions, indicates which elements had a zero derivative.
191 See Also
192 --------
193 root_scalar : interface to root solvers for scalar functions
194 root : interface to root solvers for multi-input, multi-output functions
196 Notes
197 -----
198 The convergence rate of the Newton-Raphson method is quadratic,
199 the Halley method is cubic, and the secant method is
200 sub-quadratic. This means that if the function is well-behaved
201 the actual error in the estimated root after the nth iteration
202 is approximately the square (cube for Halley) of the error
203 after the (n-1)th step. However, the stopping criterion used
204 here is the step size and there is no guarantee that a root
205 has been found. Consequently, the result should be verified.
206 Safer algorithms are brentq, brenth, ridder, and bisect,
207 but they all require that the root first be bracketed in an
208 interval where the function changes sign. The brentq algorithm
209 is recommended for general use in one dimensional problems
210 when such an interval has been found.
212 When `newton` is used with arrays, it is best suited for the following
213 types of problems:
215 * The initial guesses, `x0`, are all relatively the same distance from
216 the roots.
217 * Some or all of the extra arguments, `args`, are also arrays so that a
218 class of similar problems can be solved together.
219 * The size of the initial guesses, `x0`, is larger than O(100) elements.
220 Otherwise, a naive loop may perform as well or better than a vector.
222 Examples
223 --------
224 >>> import numpy as np
225 >>> import matplotlib.pyplot as plt
226 >>> from scipy import optimize
228 >>> def f(x):
229 ... return (x**3 - 1) # only one real root at x = 1
231 ``fprime`` is not provided, use the secant method:
233 >>> root = optimize.newton(f, 1.5)
234 >>> root
235 1.0000000000000016
236 >>> root = optimize.newton(f, 1.5, fprime2=lambda x: 6 * x)
237 >>> root
238 1.0000000000000016
240 Only ``fprime`` is provided, use the Newton-Raphson method:
242 >>> root = optimize.newton(f, 1.5, fprime=lambda x: 3 * x**2)
243 >>> root
244 1.0
246 Both ``fprime2`` and ``fprime`` are provided, use Halley's method:
248 >>> root = optimize.newton(f, 1.5, fprime=lambda x: 3 * x**2,
249 ... fprime2=lambda x: 6 * x)
250 >>> root
251 1.0
253 When we want to find roots for a set of related starting values and/or
254 function parameters, we can provide both of those as an array of inputs:
256 >>> f = lambda x, a: x**3 - a
257 >>> fder = lambda x, a: 3 * x**2
258 >>> rng = np.random.default_rng()
259 >>> x = rng.standard_normal(100)
260 >>> a = np.arange(-50, 50)
261 >>> vec_res = optimize.newton(f, x, fprime=fder, args=(a, ), maxiter=200)
263 The above is the equivalent of solving for each value in ``(x, a)``
264 separately in a for-loop, just faster:
266 >>> loop_res = [optimize.newton(f, x0, fprime=fder, args=(a0,),
267 ... maxiter=200)
268 ... for x0, a0 in zip(x, a)]
269 >>> np.allclose(vec_res, loop_res)
270 True
272 Plot the results found for all values of ``a``:
274 >>> analytical_result = np.sign(a) * np.abs(a)**(1/3)
275 >>> fig, ax = plt.subplots()
276 >>> ax.plot(a, analytical_result, 'o')
277 >>> ax.plot(a, vec_res, '.')
278 >>> ax.set_xlabel('$a$')
279 >>> ax.set_ylabel('$x$ where $f(x, a)=0$')
280 >>> plt.show()
282 """
283 if tol <= 0:
284 raise ValueError("tol too small (%g <= 0)" % tol)
285 maxiter = operator.index(maxiter)
286 if maxiter < 1:
287 raise ValueError("maxiter must be greater than 0")
288 if np.size(x0) > 1:
289 return _array_newton(func, x0, fprime, args, tol, maxiter, fprime2,
290 full_output)
292 # Convert to float (don't use float(x0); this works also for complex x0)
293 # Use np.asarray because we want x0 to be a numpy object, not a Python
294 # object. e.g. np.complex(1+1j) > 0 is possible, but (1 + 1j) > 0 raises
295 # a TypeError
296 x0 = np.asarray(x0)[()] * 1.0
297 p0 = x0
298 funcalls = 0
299 if fprime is not None:
300 # Newton-Raphson method
301 for itr in range(maxiter):
302 # first evaluate fval
303 fval = func(p0, *args)
304 funcalls += 1
305 # If fval is 0, a root has been found, then terminate
306 if fval == 0:
307 return _results_select(
308 full_output, (p0, funcalls, itr, _ECONVERGED))
309 fder = fprime(p0, *args)
310 funcalls += 1
311 if fder == 0:
312 msg = "Derivative was zero."
313 if disp:
314 msg += (
315 " Failed to converge after %d iterations, value is %s."
316 % (itr + 1, p0))
317 raise RuntimeError(msg)
318 warnings.warn(msg, RuntimeWarning)
319 return _results_select(
320 full_output, (p0, funcalls, itr + 1, _ECONVERR))
321 newton_step = fval / fder
322 if fprime2:
323 fder2 = fprime2(p0, *args)
324 funcalls += 1
325 # Halley's method:
326 # newton_step /= (1.0 - 0.5 * newton_step * fder2 / fder)
327 # Only do it if denominator stays close enough to 1
328 # Rationale: If 1-adj < 0, then Halley sends x in the
329 # opposite direction to Newton. Doesn't happen if x is close
330 # enough to root.
331 adj = newton_step * fder2 / fder / 2
332 if np.abs(adj) < 1:
333 newton_step /= 1.0 - adj
334 p = p0 - newton_step
335 if np.isclose(p, p0, rtol=rtol, atol=tol):
336 return _results_select(
337 full_output, (p, funcalls, itr + 1, _ECONVERGED))
338 p0 = p
339 else:
340 # Secant method
341 if x1 is not None:
342 if x1 == x0:
343 raise ValueError("x1 and x0 must be different")
344 p1 = x1
345 else:
346 eps = 1e-4
347 p1 = x0 * (1 + eps)
348 p1 += (eps if p1 >= 0 else -eps)
349 q0 = func(p0, *args)
350 funcalls += 1
351 q1 = func(p1, *args)
352 funcalls += 1
353 if abs(q1) < abs(q0):
354 p0, p1, q0, q1 = p1, p0, q1, q0
355 for itr in range(maxiter):
356 if q1 == q0:
357 if p1 != p0:
358 msg = "Tolerance of %s reached." % (p1 - p0)
359 if disp:
360 msg += (
361 " Failed to converge after %d iterations, value is %s."
362 % (itr + 1, p1))
363 raise RuntimeError(msg)
364 warnings.warn(msg, RuntimeWarning)
365 p = (p1 + p0) / 2.0
366 return _results_select(
367 full_output, (p, funcalls, itr + 1, _ECONVERR))
368 else:
369 if abs(q1) > abs(q0):
370 p = (-q0 / q1 * p1 + p0) / (1 - q0 / q1)
371 else:
372 p = (-q1 / q0 * p0 + p1) / (1 - q1 / q0)
373 if np.isclose(p, p1, rtol=rtol, atol=tol):
374 return _results_select(
375 full_output, (p, funcalls, itr + 1, _ECONVERGED))
376 p0, q0 = p1, q1
377 p1 = p
378 q1 = func(p1, *args)
379 funcalls += 1
381 if disp:
382 msg = ("Failed to converge after %d iterations, value is %s."
383 % (itr + 1, p))
384 raise RuntimeError(msg)
386 return _results_select(full_output, (p, funcalls, itr + 1, _ECONVERR))
389def _array_newton(func, x0, fprime, args, tol, maxiter, fprime2, full_output):
390 """
391 A vectorized version of Newton, Halley, and secant methods for arrays.
393 Do not use this method directly. This method is called from `newton`
394 when ``np.size(x0) > 1`` is ``True``. For docstring, see `newton`.
395 """
396 # Explicitly copy `x0` as `p` will be modified inplace, but the
397 # user's array should not be altered.
398 p = np.array(x0, copy=True)
400 failures = np.ones_like(p, dtype=bool)
401 nz_der = np.ones_like(failures)
402 if fprime is not None:
403 # Newton-Raphson method
404 for iteration in range(maxiter):
405 # first evaluate fval
406 fval = np.asarray(func(p, *args))
407 # If all fval are 0, all roots have been found, then terminate
408 if not fval.any():
409 failures = fval.astype(bool)
410 break
411 fder = np.asarray(fprime(p, *args))
412 nz_der = (fder != 0)
413 # stop iterating if all derivatives are zero
414 if not nz_der.any():
415 break
416 # Newton step
417 dp = fval[nz_der] / fder[nz_der]
418 if fprime2 is not None:
419 fder2 = np.asarray(fprime2(p, *args))
420 dp = dp / (1.0 - 0.5 * dp * fder2[nz_der] / fder[nz_der])
421 # only update nonzero derivatives
422 p = np.asarray(p, dtype=np.result_type(p, dp, np.float64))
423 p[nz_der] -= dp
424 failures[nz_der] = np.abs(dp) >= tol # items not yet converged
425 # stop iterating if there aren't any failures, not incl zero der
426 if not failures[nz_der].any():
427 break
428 else:
429 # Secant method
430 dx = np.finfo(float).eps**0.33
431 p1 = p * (1 + dx) + np.where(p >= 0, dx, -dx)
432 q0 = np.asarray(func(p, *args))
433 q1 = np.asarray(func(p1, *args))
434 active = np.ones_like(p, dtype=bool)
435 for iteration in range(maxiter):
436 nz_der = (q1 != q0)
437 # stop iterating if all derivatives are zero
438 if not nz_der.any():
439 p = (p1 + p) / 2.0
440 break
441 # Secant Step
442 dp = (q1 * (p1 - p))[nz_der] / (q1 - q0)[nz_der]
443 # only update nonzero derivatives
444 p = np.asarray(p, dtype=np.result_type(p, p1, dp, np.float64))
445 p[nz_der] = p1[nz_der] - dp
446 active_zero_der = ~nz_der & active
447 p[active_zero_der] = (p1 + p)[active_zero_der] / 2.0
448 active &= nz_der # don't assign zero derivatives again
449 failures[nz_der] = np.abs(dp) >= tol # not yet converged
450 # stop iterating if there aren't any failures, not incl zero der
451 if not failures[nz_der].any():
452 break
453 p1, p = p, p1
454 q0 = q1
455 q1 = np.asarray(func(p1, *args))
457 zero_der = ~nz_der & failures # don't include converged with zero-ders
458 if zero_der.any():
459 # Secant warnings
460 if fprime is None:
461 nonzero_dp = (p1 != p)
462 # non-zero dp, but infinite newton step
463 zero_der_nz_dp = (zero_der & nonzero_dp)
464 if zero_der_nz_dp.any():
465 rms = np.sqrt(
466 sum((p1[zero_der_nz_dp] - p[zero_der_nz_dp]) ** 2)
467 )
468 warnings.warn(
469 f'RMS of {rms:g} reached', RuntimeWarning)
470 # Newton or Halley warnings
471 else:
472 all_or_some = 'all' if zero_der.all() else 'some'
473 msg = f'{all_or_some:s} derivatives were zero'
474 warnings.warn(msg, RuntimeWarning)
475 elif failures.any():
476 all_or_some = 'all' if failures.all() else 'some'
477 msg = '{:s} failed to converge after {:d} iterations'.format(
478 all_or_some, maxiter
479 )
480 if failures.all():
481 raise RuntimeError(msg)
482 warnings.warn(msg, RuntimeWarning)
484 if full_output:
485 result = namedtuple('result', ('root', 'converged', 'zero_der'))
486 p = result(p, ~failures, zero_der)
488 return p
491def bisect(f, a, b, args=(),
492 xtol=_xtol, rtol=_rtol, maxiter=_iter,
493 full_output=False, disp=True):
494 """
495 Find root of a function within an interval using bisection.
497 Basic bisection routine to find a root of the function `f` between the
498 arguments `a` and `b`. `f(a)` and `f(b)` cannot have the same signs.
499 Slow but sure.
501 Parameters
502 ----------
503 f : function
504 Python function returning a number. `f` must be continuous, and
505 f(a) and f(b) must have opposite signs.
506 a : scalar
507 One end of the bracketing interval [a,b].
508 b : scalar
509 The other end of the bracketing interval [a,b].
510 xtol : number, optional
511 The computed root ``x0`` will satisfy ``np.allclose(x, x0,
512 atol=xtol, rtol=rtol)``, where ``x`` is the exact root. The
513 parameter must be positive.
514 rtol : number, optional
515 The computed root ``x0`` will satisfy ``np.allclose(x, x0,
516 atol=xtol, rtol=rtol)``, where ``x`` is the exact root. The
517 parameter cannot be smaller than its default value of
518 ``4*np.finfo(float).eps``.
519 maxiter : int, optional
520 If convergence is not achieved in `maxiter` iterations, an error is
521 raised. Must be >= 0.
522 args : tuple, optional
523 Containing extra arguments for the function `f`.
524 `f` is called by ``apply(f, (x)+args)``.
525 full_output : bool, optional
526 If `full_output` is False, the root is returned. If `full_output` is
527 True, the return value is ``(x, r)``, where x is the root, and r is
528 a `RootResults` object.
529 disp : bool, optional
530 If True, raise RuntimeError if the algorithm didn't converge.
531 Otherwise, the convergence status is recorded in a `RootResults`
532 return object.
534 Returns
535 -------
536 root : float
537 Root of `f` between `a` and `b`.
538 r : `RootResults` (present if ``full_output = True``)
539 Object containing information about the convergence. In particular,
540 ``r.converged`` is True if the routine converged.
542 Examples
543 --------
545 >>> def f(x):
546 ... return (x**2 - 1)
548 >>> from scipy import optimize
550 >>> root = optimize.bisect(f, 0, 2)
551 >>> root
552 1.0
554 >>> root = optimize.bisect(f, -2, 0)
555 >>> root
556 -1.0
558 See Also
559 --------
560 brentq, brenth, bisect, newton
561 fixed_point : scalar fixed-point finder
562 fsolve : n-dimensional root-finding
564 """
565 if not isinstance(args, tuple):
566 args = (args,)
567 maxiter = operator.index(maxiter)
568 if xtol <= 0:
569 raise ValueError("xtol too small (%g <= 0)" % xtol)
570 if rtol < _rtol:
571 raise ValueError(f"rtol too small ({rtol:g} < {_rtol:g})")
572 f = _wrap_nan_raise(f)
573 r = _zeros._bisect(f, a, b, xtol, rtol, maxiter, args, full_output, disp)
574 return results_c(full_output, r)
577def ridder(f, a, b, args=(),
578 xtol=_xtol, rtol=_rtol, maxiter=_iter,
579 full_output=False, disp=True):
580 """
581 Find a root of a function in an interval using Ridder's method.
583 Parameters
584 ----------
585 f : function
586 Python function returning a number. f must be continuous, and f(a) and
587 f(b) must have opposite signs.
588 a : scalar
589 One end of the bracketing interval [a,b].
590 b : scalar
591 The other end of the bracketing interval [a,b].
592 xtol : number, optional
593 The computed root ``x0`` will satisfy ``np.allclose(x, x0,
594 atol=xtol, rtol=rtol)``, where ``x`` is the exact root. The
595 parameter must be positive.
596 rtol : number, optional
597 The computed root ``x0`` will satisfy ``np.allclose(x, x0,
598 atol=xtol, rtol=rtol)``, where ``x`` is the exact root. The
599 parameter cannot be smaller than its default value of
600 ``4*np.finfo(float).eps``.
601 maxiter : int, optional
602 If convergence is not achieved in `maxiter` iterations, an error is
603 raised. Must be >= 0.
604 args : tuple, optional
605 Containing extra arguments for the function `f`.
606 `f` is called by ``apply(f, (x)+args)``.
607 full_output : bool, optional
608 If `full_output` is False, the root is returned. If `full_output` is
609 True, the return value is ``(x, r)``, where `x` is the root, and `r` is
610 a `RootResults` object.
611 disp : bool, optional
612 If True, raise RuntimeError if the algorithm didn't converge.
613 Otherwise, the convergence status is recorded in any `RootResults`
614 return object.
616 Returns
617 -------
618 root : float
619 Root of `f` between `a` and `b`.
620 r : `RootResults` (present if ``full_output = True``)
621 Object containing information about the convergence.
622 In particular, ``r.converged`` is True if the routine converged.
624 See Also
625 --------
626 brentq, brenth, bisect, newton : 1-D root-finding
627 fixed_point : scalar fixed-point finder
629 Notes
630 -----
631 Uses [Ridders1979]_ method to find a root of the function `f` between the
632 arguments `a` and `b`. Ridders' method is faster than bisection, but not
633 generally as fast as the Brent routines. [Ridders1979]_ provides the
634 classic description and source of the algorithm. A description can also be
635 found in any recent edition of Numerical Recipes.
637 The routine used here diverges slightly from standard presentations in
638 order to be a bit more careful of tolerance.
640 References
641 ----------
642 .. [Ridders1979]
643 Ridders, C. F. J. "A New Algorithm for Computing a
644 Single Root of a Real Continuous Function."
645 IEEE Trans. Circuits Systems 26, 979-980, 1979.
647 Examples
648 --------
650 >>> def f(x):
651 ... return (x**2 - 1)
653 >>> from scipy import optimize
655 >>> root = optimize.ridder(f, 0, 2)
656 >>> root
657 1.0
659 >>> root = optimize.ridder(f, -2, 0)
660 >>> root
661 -1.0
662 """
663 if not isinstance(args, tuple):
664 args = (args,)
665 maxiter = operator.index(maxiter)
666 if xtol <= 0:
667 raise ValueError("xtol too small (%g <= 0)" % xtol)
668 if rtol < _rtol:
669 raise ValueError(f"rtol too small ({rtol:g} < {_rtol:g})")
670 f = _wrap_nan_raise(f)
671 r = _zeros._ridder(f, a, b, xtol, rtol, maxiter, args, full_output, disp)
672 return results_c(full_output, r)
675def brentq(f, a, b, args=(),
676 xtol=_xtol, rtol=_rtol, maxiter=_iter,
677 full_output=False, disp=True):
678 """
679 Find a root of a function in a bracketing interval using Brent's method.
681 Uses the classic Brent's method to find a root of the function `f` on
682 the sign changing interval [a , b]. Generally considered the best of the
683 rootfinding routines here. It is a safe version of the secant method that
684 uses inverse quadratic extrapolation. Brent's method combines root
685 bracketing, interval bisection, and inverse quadratic interpolation. It is
686 sometimes known as the van Wijngaarden-Dekker-Brent method. Brent (1973)
687 claims convergence is guaranteed for functions computable within [a,b].
689 [Brent1973]_ provides the classic description of the algorithm. Another
690 description can be found in a recent edition of Numerical Recipes, including
691 [PressEtal1992]_. A third description is at
692 http://mathworld.wolfram.com/BrentsMethod.html. It should be easy to
693 understand the algorithm just by reading our code. Our code diverges a bit
694 from standard presentations: we choose a different formula for the
695 extrapolation step.
697 Parameters
698 ----------
699 f : function
700 Python function returning a number. The function :math:`f`
701 must be continuous, and :math:`f(a)` and :math:`f(b)` must
702 have opposite signs.
703 a : scalar
704 One end of the bracketing interval :math:`[a, b]`.
705 b : scalar
706 The other end of the bracketing interval :math:`[a, b]`.
707 xtol : number, optional
708 The computed root ``x0`` will satisfy ``np.allclose(x, x0,
709 atol=xtol, rtol=rtol)``, where ``x`` is the exact root. The
710 parameter must be positive. For nice functions, Brent's
711 method will often satisfy the above condition with ``xtol/2``
712 and ``rtol/2``. [Brent1973]_
713 rtol : number, optional
714 The computed root ``x0`` will satisfy ``np.allclose(x, x0,
715 atol=xtol, rtol=rtol)``, where ``x`` is the exact root. The
716 parameter cannot be smaller than its default value of
717 ``4*np.finfo(float).eps``. For nice functions, Brent's
718 method will often satisfy the above condition with ``xtol/2``
719 and ``rtol/2``. [Brent1973]_
720 maxiter : int, optional
721 If convergence is not achieved in `maxiter` iterations, an error is
722 raised. Must be >= 0.
723 args : tuple, optional
724 Containing extra arguments for the function `f`.
725 `f` is called by ``apply(f, (x)+args)``.
726 full_output : bool, optional
727 If `full_output` is False, the root is returned. If `full_output` is
728 True, the return value is ``(x, r)``, where `x` is the root, and `r` is
729 a `RootResults` object.
730 disp : bool, optional
731 If True, raise RuntimeError if the algorithm didn't converge.
732 Otherwise, the convergence status is recorded in any `RootResults`
733 return object.
735 Returns
736 -------
737 root : float
738 Root of `f` between `a` and `b`.
739 r : `RootResults` (present if ``full_output = True``)
740 Object containing information about the convergence. In particular,
741 ``r.converged`` is True if the routine converged.
743 Notes
744 -----
745 `f` must be continuous. f(a) and f(b) must have opposite signs.
747 Related functions fall into several classes:
749 multivariate local optimizers
750 `fmin`, `fmin_powell`, `fmin_cg`, `fmin_bfgs`, `fmin_ncg`
751 nonlinear least squares minimizer
752 `leastsq`
753 constrained multivariate optimizers
754 `fmin_l_bfgs_b`, `fmin_tnc`, `fmin_cobyla`
755 global optimizers
756 `basinhopping`, `brute`, `differential_evolution`
757 local scalar minimizers
758 `fminbound`, `brent`, `golden`, `bracket`
759 N-D root-finding
760 `fsolve`
761 1-D root-finding
762 `brenth`, `ridder`, `bisect`, `newton`
763 scalar fixed-point finder
764 `fixed_point`
766 References
767 ----------
768 .. [Brent1973]
769 Brent, R. P.,
770 *Algorithms for Minimization Without Derivatives*.
771 Englewood Cliffs, NJ: Prentice-Hall, 1973. Ch. 3-4.
773 .. [PressEtal1992]
774 Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; and Vetterling, W. T.
775 *Numerical Recipes in FORTRAN: The Art of Scientific Computing*, 2nd ed.
776 Cambridge, England: Cambridge University Press, pp. 352-355, 1992.
777 Section 9.3: "Van Wijngaarden-Dekker-Brent Method."
779 Examples
780 --------
781 >>> def f(x):
782 ... return (x**2 - 1)
784 >>> from scipy import optimize
786 >>> root = optimize.brentq(f, -2, 0)
787 >>> root
788 -1.0
790 >>> root = optimize.brentq(f, 0, 2)
791 >>> root
792 1.0
793 """
794 if not isinstance(args, tuple):
795 args = (args,)
796 maxiter = operator.index(maxiter)
797 if xtol <= 0:
798 raise ValueError("xtol too small (%g <= 0)" % xtol)
799 if rtol < _rtol:
800 raise ValueError(f"rtol too small ({rtol:g} < {_rtol:g})")
801 f = _wrap_nan_raise(f)
802 r = _zeros._brentq(f, a, b, xtol, rtol, maxiter, args, full_output, disp)
803 return results_c(full_output, r)
806def brenth(f, a, b, args=(),
807 xtol=_xtol, rtol=_rtol, maxiter=_iter,
808 full_output=False, disp=True):
809 """Find a root of a function in a bracketing interval using Brent's
810 method with hyperbolic extrapolation.
812 A variation on the classic Brent routine to find a root of the function f
813 between the arguments a and b that uses hyperbolic extrapolation instead of
814 inverse quadratic extrapolation. Bus & Dekker (1975) guarantee convergence
815 for this method, claiming that the upper bound of function evaluations here
816 is 4 or 5 times lesser than that for bisection.
817 f(a) and f(b) cannot have the same signs. Generally, on a par with the
818 brent routine, but not as heavily tested. It is a safe version of the
819 secant method that uses hyperbolic extrapolation.
820 The version here is by Chuck Harris, and implements Algorithm M of
821 [BusAndDekker1975]_, where further details (convergence properties,
822 additional remarks and such) can be found
824 Parameters
825 ----------
826 f : function
827 Python function returning a number. f must be continuous, and f(a) and
828 f(b) must have opposite signs.
829 a : scalar
830 One end of the bracketing interval [a,b].
831 b : scalar
832 The other end of the bracketing interval [a,b].
833 xtol : number, optional
834 The computed root ``x0`` will satisfy ``np.allclose(x, x0,
835 atol=xtol, rtol=rtol)``, where ``x`` is the exact root. The
836 parameter must be positive. As with `brentq`, for nice
837 functions the method will often satisfy the above condition
838 with ``xtol/2`` and ``rtol/2``.
839 rtol : number, optional
840 The computed root ``x0`` will satisfy ``np.allclose(x, x0,
841 atol=xtol, rtol=rtol)``, where ``x`` is the exact root. The
842 parameter cannot be smaller than its default value of
843 ``4*np.finfo(float).eps``. As with `brentq`, for nice functions
844 the method will often satisfy the above condition with
845 ``xtol/2`` and ``rtol/2``.
846 maxiter : int, optional
847 If convergence is not achieved in `maxiter` iterations, an error is
848 raised. Must be >= 0.
849 args : tuple, optional
850 Containing extra arguments for the function `f`.
851 `f` is called by ``apply(f, (x)+args)``.
852 full_output : bool, optional
853 If `full_output` is False, the root is returned. If `full_output` is
854 True, the return value is ``(x, r)``, where `x` is the root, and `r` is
855 a `RootResults` object.
856 disp : bool, optional
857 If True, raise RuntimeError if the algorithm didn't converge.
858 Otherwise, the convergence status is recorded in any `RootResults`
859 return object.
861 Returns
862 -------
863 root : float
864 Root of `f` between `a` and `b`.
865 r : `RootResults` (present if ``full_output = True``)
866 Object containing information about the convergence. In particular,
867 ``r.converged`` is True if the routine converged.
869 See Also
870 --------
871 fmin, fmin_powell, fmin_cg, fmin_bfgs, fmin_ncg : multivariate local optimizers
872 leastsq : nonlinear least squares minimizer
873 fmin_l_bfgs_b, fmin_tnc, fmin_cobyla : constrained multivariate optimizers
874 basinhopping, differential_evolution, brute : global optimizers
875 fminbound, brent, golden, bracket : local scalar minimizers
876 fsolve : N-D root-finding
877 brentq, brenth, ridder, bisect, newton : 1-D root-finding
878 fixed_point : scalar fixed-point finder
880 References
881 ----------
882 .. [BusAndDekker1975]
883 Bus, J. C. P., Dekker, T. J.,
884 "Two Efficient Algorithms with Guaranteed Convergence for Finding a Zero
885 of a Function", ACM Transactions on Mathematical Software, Vol. 1, Issue
886 4, Dec. 1975, pp. 330-345. Section 3: "Algorithm M".
887 :doi:`10.1145/355656.355659`
889 Examples
890 --------
891 >>> def f(x):
892 ... return (x**2 - 1)
894 >>> from scipy import optimize
896 >>> root = optimize.brenth(f, -2, 0)
897 >>> root
898 -1.0
900 >>> root = optimize.brenth(f, 0, 2)
901 >>> root
902 1.0
904 """
905 if not isinstance(args, tuple):
906 args = (args,)
907 maxiter = operator.index(maxiter)
908 if xtol <= 0:
909 raise ValueError("xtol too small (%g <= 0)" % xtol)
910 if rtol < _rtol:
911 raise ValueError(f"rtol too small ({rtol:g} < {_rtol:g})")
912 f = _wrap_nan_raise(f)
913 r = _zeros._brenth(f, a, b, xtol, rtol, maxiter, args, full_output, disp)
914 return results_c(full_output, r)
917################################
918# TOMS "Algorithm 748: Enclosing Zeros of Continuous Functions", by
919# Alefeld, G. E. and Potra, F. A. and Shi, Yixun,
920# See [1]
923def _notclose(fs, rtol=_rtol, atol=_xtol):
924 # Ensure not None, not 0, all finite, and not very close to each other
925 notclosefvals = (
926 all(fs) and all(np.isfinite(fs)) and
927 not any(any(np.isclose(_f, fs[i + 1:], rtol=rtol, atol=atol))
928 for i, _f in enumerate(fs[:-1])))
929 return notclosefvals
932def _secant(xvals, fvals):
933 """Perform a secant step, taking a little care"""
934 # Secant has many "mathematically" equivalent formulations
935 # x2 = x0 - (x1 - x0)/(f1 - f0) * f0
936 # = x1 - (x1 - x0)/(f1 - f0) * f1
937 # = (-x1 * f0 + x0 * f1) / (f1 - f0)
938 # = (-f0 / f1 * x1 + x0) / (1 - f0 / f1)
939 # = (-f1 / f0 * x0 + x1) / (1 - f1 / f0)
940 x0, x1 = xvals[:2]
941 f0, f1 = fvals[:2]
942 if f0 == f1:
943 return np.nan
944 if np.abs(f1) > np.abs(f0):
945 x2 = (-f0 / f1 * x1 + x0) / (1 - f0 / f1)
946 else:
947 x2 = (-f1 / f0 * x0 + x1) / (1 - f1 / f0)
948 return x2
951def _update_bracket(ab, fab, c, fc):
952 """Update a bracket given (c, fc), return the discarded endpoints."""
953 fa, fb = fab
954 idx = (0 if np.sign(fa) * np.sign(fc) > 0 else 1)
955 rx, rfx = ab[idx], fab[idx]
956 fab[idx] = fc
957 ab[idx] = c
958 return rx, rfx
961def _compute_divided_differences(xvals, fvals, N=None, full=True,
962 forward=True):
963 """Return a matrix of divided differences for the xvals, fvals pairs
965 DD[i, j] = f[x_{i-j}, ..., x_i] for 0 <= j <= i
967 If full is False, just return the main diagonal(or last row):
968 f[a], f[a, b] and f[a, b, c].
969 If forward is False, return f[c], f[b, c], f[a, b, c]."""
970 if full:
971 if forward:
972 xvals = np.asarray(xvals)
973 else:
974 xvals = np.array(xvals)[::-1]
975 M = len(xvals)
976 N = M if N is None else min(N, M)
977 DD = np.zeros([M, N])
978 DD[:, 0] = fvals[:]
979 for i in range(1, N):
980 DD[i:, i] = (np.diff(DD[i - 1:, i - 1]) /
981 (xvals[i:] - xvals[:M - i]))
982 return DD
984 xvals = np.asarray(xvals)
985 dd = np.array(fvals)
986 row = np.array(fvals)
987 idx2Use = (0 if forward else -1)
988 dd[0] = fvals[idx2Use]
989 for i in range(1, len(xvals)):
990 denom = xvals[i:i + len(row) - 1] - xvals[:len(row) - 1]
991 row = np.diff(row)[:] / denom
992 dd[i] = row[idx2Use]
993 return dd
996def _interpolated_poly(xvals, fvals, x):
997 """Compute p(x) for the polynomial passing through the specified locations.
999 Use Neville's algorithm to compute p(x) where p is the minimal degree
1000 polynomial passing through the points xvals, fvals"""
1001 xvals = np.asarray(xvals)
1002 N = len(xvals)
1003 Q = np.zeros([N, N])
1004 D = np.zeros([N, N])
1005 Q[:, 0] = fvals[:]
1006 D[:, 0] = fvals[:]
1007 for k in range(1, N):
1008 alpha = D[k:, k - 1] - Q[k - 1:N - 1, k - 1]
1009 diffik = xvals[0:N - k] - xvals[k:N]
1010 Q[k:, k] = (xvals[k:] - x) / diffik * alpha
1011 D[k:, k] = (xvals[:N - k] - x) / diffik * alpha
1012 # Expect Q[-1, 1:] to be small relative to Q[-1, 0] as x approaches a root
1013 return np.sum(Q[-1, 1:]) + Q[-1, 0]
1016def _inverse_poly_zero(a, b, c, d, fa, fb, fc, fd):
1017 """Inverse cubic interpolation f-values -> x-values
1019 Given four points (fa, a), (fb, b), (fc, c), (fd, d) with
1020 fa, fb, fc, fd all distinct, find poly IP(y) through the 4 points
1021 and compute x=IP(0).
1022 """
1023 return _interpolated_poly([fa, fb, fc, fd], [a, b, c, d], 0)
1026def _newton_quadratic(ab, fab, d, fd, k):
1027 """Apply Newton-Raphson like steps, using divided differences to approximate f'
1029 ab is a real interval [a, b] containing a root,
1030 fab holds the real values of f(a), f(b)
1031 d is a real number outside [ab, b]
1032 k is the number of steps to apply
1033 """
1034 a, b = ab
1035 fa, fb = fab
1036 _, B, A = _compute_divided_differences([a, b, d], [fa, fb, fd],
1037 forward=True, full=False)
1039 # _P is the quadratic polynomial through the 3 points
1040 def _P(x):
1041 # Horner evaluation of fa + B * (x - a) + A * (x - a) * (x - b)
1042 return (A * (x - b) + B) * (x - a) + fa
1044 if A == 0:
1045 r = a - fa / B
1046 else:
1047 r = (a if np.sign(A) * np.sign(fa) > 0 else b)
1048 # Apply k Newton-Raphson steps to _P(x), starting from x=r
1049 for i in range(k):
1050 r1 = r - _P(r) / (B + A * (2 * r - a - b))
1051 if not (ab[0] < r1 < ab[1]):
1052 if (ab[0] < r < ab[1]):
1053 return r
1054 r = sum(ab) / 2.0
1055 break
1056 r = r1
1058 return r
1061class TOMS748Solver:
1062 """Solve f(x, *args) == 0 using Algorithm748 of Alefeld, Potro & Shi.
1063 """
1064 _MU = 0.5
1065 _K_MIN = 1
1066 _K_MAX = 100 # A very high value for real usage. Expect 1, 2, maybe 3.
1068 def __init__(self):
1069 self.f = None
1070 self.args = None
1071 self.function_calls = 0
1072 self.iterations = 0
1073 self.k = 2
1074 # ab=[a,b] is a global interval containing a root
1075 self.ab = [np.nan, np.nan]
1076 # fab is function values at a, b
1077 self.fab = [np.nan, np.nan]
1078 self.d = None
1079 self.fd = None
1080 self.e = None
1081 self.fe = None
1082 self.disp = False
1083 self.xtol = _xtol
1084 self.rtol = _rtol
1085 self.maxiter = _iter
1087 def configure(self, xtol, rtol, maxiter, disp, k):
1088 self.disp = disp
1089 self.xtol = xtol
1090 self.rtol = rtol
1091 self.maxiter = maxiter
1092 # Silently replace a low value of k with 1
1093 self.k = max(k, self._K_MIN)
1094 # Noisily replace a high value of k with self._K_MAX
1095 if self.k > self._K_MAX:
1096 msg = "toms748: Overriding k: ->%d" % self._K_MAX
1097 warnings.warn(msg, RuntimeWarning)
1098 self.k = self._K_MAX
1100 def _callf(self, x, error=True):
1101 """Call the user-supplied function, update book-keeping"""
1102 fx = self.f(x, *self.args)
1103 self.function_calls += 1
1104 if not np.isfinite(fx) and error:
1105 raise ValueError(f"Invalid function value: f({x:f}) -> {fx} ")
1106 return fx
1108 def get_result(self, x, flag=_ECONVERGED):
1109 r"""Package the result and statistics into a tuple."""
1110 return (x, self.function_calls, self.iterations, flag)
1112 def _update_bracket(self, c, fc):
1113 return _update_bracket(self.ab, self.fab, c, fc)
1115 def start(self, f, a, b, args=()):
1116 r"""Prepare for the iterations."""
1117 self.function_calls = 0
1118 self.iterations = 0
1120 self.f = f
1121 self.args = args
1122 self.ab[:] = [a, b]
1123 if not np.isfinite(a) or np.imag(a) != 0:
1124 raise ValueError("Invalid x value: %s " % (a))
1125 if not np.isfinite(b) or np.imag(b) != 0:
1126 raise ValueError("Invalid x value: %s " % (b))
1128 fa = self._callf(a)
1129 if not np.isfinite(fa) or np.imag(fa) != 0:
1130 raise ValueError(f"Invalid function value: f({a:f}) -> {fa} ")
1131 if fa == 0:
1132 return _ECONVERGED, a
1133 fb = self._callf(b)
1134 if not np.isfinite(fb) or np.imag(fb) != 0:
1135 raise ValueError(f"Invalid function value: f({b:f}) -> {fb} ")
1136 if fb == 0:
1137 return _ECONVERGED, b
1139 if np.sign(fb) * np.sign(fa) > 0:
1140 raise ValueError("f(a) and f(b) must have different signs, but "
1141 "f({:e})={:e}, f({:e})={:e} ".format(a, fa, b, fb))
1142 self.fab[:] = [fa, fb]
1144 return _EINPROGRESS, sum(self.ab) / 2.0
1146 def get_status(self):
1147 """Determine the current status."""
1148 a, b = self.ab[:2]
1149 if np.isclose(a, b, rtol=self.rtol, atol=self.xtol):
1150 return _ECONVERGED, sum(self.ab) / 2.0
1151 if self.iterations >= self.maxiter:
1152 return _ECONVERR, sum(self.ab) / 2.0
1153 return _EINPROGRESS, sum(self.ab) / 2.0
1155 def iterate(self):
1156 """Perform one step in the algorithm.
1158 Implements Algorithm 4.1(k=1) or 4.2(k=2) in [APS1995]
1159 """
1160 self.iterations += 1
1161 eps = np.finfo(float).eps
1162 d, fd, e, fe = self.d, self.fd, self.e, self.fe
1163 ab_width = self.ab[1] - self.ab[0] # Need the start width below
1164 c = None
1166 for nsteps in range(2, self.k+2):
1167 # If the f-values are sufficiently separated, perform an inverse
1168 # polynomial interpolation step. Otherwise, nsteps repeats of
1169 # an approximate Newton-Raphson step.
1170 if _notclose(self.fab + [fd, fe], rtol=0, atol=32*eps):
1171 c0 = _inverse_poly_zero(self.ab[0], self.ab[1], d, e,
1172 self.fab[0], self.fab[1], fd, fe)
1173 if self.ab[0] < c0 < self.ab[1]:
1174 c = c0
1175 if c is None:
1176 c = _newton_quadratic(self.ab, self.fab, d, fd, nsteps)
1178 fc = self._callf(c)
1179 if fc == 0:
1180 return _ECONVERGED, c
1182 # re-bracket
1183 e, fe = d, fd
1184 d, fd = self._update_bracket(c, fc)
1186 # u is the endpoint with the smallest f-value
1187 uix = (0 if np.abs(self.fab[0]) < np.abs(self.fab[1]) else 1)
1188 u, fu = self.ab[uix], self.fab[uix]
1190 _, A = _compute_divided_differences(self.ab, self.fab,
1191 forward=(uix == 0), full=False)
1192 c = u - 2 * fu / A
1193 if np.abs(c - u) > 0.5 * (self.ab[1] - self.ab[0]):
1194 c = sum(self.ab) / 2.0
1195 else:
1196 if np.isclose(c, u, rtol=eps, atol=0):
1197 # c didn't change (much).
1198 # Either because the f-values at the endpoints have vastly
1199 # differing magnitudes, or because the root is very close to
1200 # that endpoint
1201 frs = np.frexp(self.fab)[1]
1202 if frs[uix] < frs[1 - uix] - 50: # Differ by more than 2**50
1203 c = (31 * self.ab[uix] + self.ab[1 - uix]) / 32
1204 else:
1205 # Make a bigger adjustment, about the
1206 # size of the requested tolerance.
1207 mm = (1 if uix == 0 else -1)
1208 adj = mm * np.abs(c) * self.rtol + mm * self.xtol
1209 c = u + adj
1210 if not self.ab[0] < c < self.ab[1]:
1211 c = sum(self.ab) / 2.0
1213 fc = self._callf(c)
1214 if fc == 0:
1215 return _ECONVERGED, c
1217 e, fe = d, fd
1218 d, fd = self._update_bracket(c, fc)
1220 # If the width of the new interval did not decrease enough, bisect
1221 if self.ab[1] - self.ab[0] > self._MU * ab_width:
1222 e, fe = d, fd
1223 z = sum(self.ab) / 2.0
1224 fz = self._callf(z)
1225 if fz == 0:
1226 return _ECONVERGED, z
1227 d, fd = self._update_bracket(z, fz)
1229 # Record d and e for next iteration
1230 self.d, self.fd = d, fd
1231 self.e, self.fe = e, fe
1233 status, xn = self.get_status()
1234 return status, xn
1236 def solve(self, f, a, b, args=(),
1237 xtol=_xtol, rtol=_rtol, k=2, maxiter=_iter, disp=True):
1238 r"""Solve f(x) = 0 given an interval containing a root."""
1239 self.configure(xtol=xtol, rtol=rtol, maxiter=maxiter, disp=disp, k=k)
1240 status, xn = self.start(f, a, b, args)
1241 if status == _ECONVERGED:
1242 return self.get_result(xn)
1244 # The first step only has two x-values.
1245 c = _secant(self.ab, self.fab)
1246 if not self.ab[0] < c < self.ab[1]:
1247 c = sum(self.ab) / 2.0
1248 fc = self._callf(c)
1249 if fc == 0:
1250 return self.get_result(c)
1252 self.d, self.fd = self._update_bracket(c, fc)
1253 self.e, self.fe = None, None
1254 self.iterations += 1
1256 while True:
1257 status, xn = self.iterate()
1258 if status == _ECONVERGED:
1259 return self.get_result(xn)
1260 if status == _ECONVERR:
1261 fmt = "Failed to converge after %d iterations, bracket is %s"
1262 if disp:
1263 msg = fmt % (self.iterations + 1, self.ab)
1264 raise RuntimeError(msg)
1265 return self.get_result(xn, _ECONVERR)
1268def toms748(f, a, b, args=(), k=1,
1269 xtol=_xtol, rtol=_rtol, maxiter=_iter,
1270 full_output=False, disp=True):
1271 """
1272 Find a root using TOMS Algorithm 748 method.
1274 Implements the Algorithm 748 method of Alefeld, Potro and Shi to find a
1275 root of the function `f` on the interval `[a , b]`, where `f(a)` and
1276 `f(b)` must have opposite signs.
1278 It uses a mixture of inverse cubic interpolation and
1279 "Newton-quadratic" steps. [APS1995].
1281 Parameters
1282 ----------
1283 f : function
1284 Python function returning a scalar. The function :math:`f`
1285 must be continuous, and :math:`f(a)` and :math:`f(b)`
1286 have opposite signs.
1287 a : scalar,
1288 lower boundary of the search interval
1289 b : scalar,
1290 upper boundary of the search interval
1291 args : tuple, optional
1292 containing extra arguments for the function `f`.
1293 `f` is called by ``f(x, *args)``.
1294 k : int, optional
1295 The number of Newton quadratic steps to perform each
1296 iteration. ``k>=1``.
1297 xtol : scalar, optional
1298 The computed root ``x0`` will satisfy ``np.allclose(x, x0,
1299 atol=xtol, rtol=rtol)``, where ``x`` is the exact root. The
1300 parameter must be positive.
1301 rtol : scalar, optional
1302 The computed root ``x0`` will satisfy ``np.allclose(x, x0,
1303 atol=xtol, rtol=rtol)``, where ``x`` is the exact root.
1304 maxiter : int, optional
1305 If convergence is not achieved in `maxiter` iterations, an error is
1306 raised. Must be >= 0.
1307 full_output : bool, optional
1308 If `full_output` is False, the root is returned. If `full_output` is
1309 True, the return value is ``(x, r)``, where `x` is the root, and `r` is
1310 a `RootResults` object.
1311 disp : bool, optional
1312 If True, raise RuntimeError if the algorithm didn't converge.
1313 Otherwise, the convergence status is recorded in the `RootResults`
1314 return object.
1316 Returns
1317 -------
1318 root : float
1319 Approximate root of `f`
1320 r : `RootResults` (present if ``full_output = True``)
1321 Object containing information about the convergence. In particular,
1322 ``r.converged`` is True if the routine converged.
1324 See Also
1325 --------
1326 brentq, brenth, ridder, bisect, newton
1327 fsolve : find roots in N dimensions.
1329 Notes
1330 -----
1331 `f` must be continuous.
1332 Algorithm 748 with ``k=2`` is asymptotically the most efficient
1333 algorithm known for finding roots of a four times continuously
1334 differentiable function.
1335 In contrast with Brent's algorithm, which may only decrease the length of
1336 the enclosing bracket on the last step, Algorithm 748 decreases it each
1337 iteration with the same asymptotic efficiency as it finds the root.
1339 For easy statement of efficiency indices, assume that `f` has 4
1340 continuouous deriviatives.
1341 For ``k=1``, the convergence order is at least 2.7, and with about
1342 asymptotically 2 function evaluations per iteration, the efficiency
1343 index is approximately 1.65.
1344 For ``k=2``, the order is about 4.6 with asymptotically 3 function
1345 evaluations per iteration, and the efficiency index 1.66.
1346 For higher values of `k`, the efficiency index approaches
1347 the kth root of ``(3k-2)``, hence ``k=1`` or ``k=2`` are
1348 usually appropriate.
1350 References
1351 ----------
1352 .. [APS1995]
1353 Alefeld, G. E. and Potra, F. A. and Shi, Yixun,
1354 *Algorithm 748: Enclosing Zeros of Continuous Functions*,
1355 ACM Trans. Math. Softw. Volume 221(1995)
1356 doi = {10.1145/210089.210111}
1358 Examples
1359 --------
1360 >>> def f(x):
1361 ... return (x**3 - 1) # only one real root at x = 1
1363 >>> from scipy import optimize
1364 >>> root, results = optimize.toms748(f, 0, 2, full_output=True)
1365 >>> root
1366 1.0
1367 >>> results
1368 converged: True
1369 flag: converged
1370 function_calls: 11
1371 iterations: 5
1372 root: 1.0
1373 """
1374 if xtol <= 0:
1375 raise ValueError("xtol too small (%g <= 0)" % xtol)
1376 if rtol < _rtol / 4:
1377 raise ValueError(f"rtol too small ({rtol:g} < {_rtol/4:g})")
1378 maxiter = operator.index(maxiter)
1379 if maxiter < 1:
1380 raise ValueError("maxiter must be greater than 0")
1381 if not np.isfinite(a):
1382 raise ValueError("a is not finite %s" % a)
1383 if not np.isfinite(b):
1384 raise ValueError("b is not finite %s" % b)
1385 if a >= b:
1386 raise ValueError(f"a and b are not an interval [{a}, {b}]")
1387 if not k >= 1:
1388 raise ValueError("k too small (%s < 1)" % k)
1390 if not isinstance(args, tuple):
1391 args = (args,)
1392 f = _wrap_nan_raise(f)
1393 solver = TOMS748Solver()
1394 result = solver.solve(f, a, b, args=args, k=k, xtol=xtol, rtol=rtol,
1395 maxiter=maxiter, disp=disp)
1396 x, function_calls, iterations, flag = result
1397 return _results_select(full_output, (x, function_calls, iterations, flag))