1 [section:roots_noderiv Root Finding Without Derivatives]
6 #include <boost/math/tools/roots.hpp>
9 namespace boost { namespace math {
10 namespace tools { // Note namespace boost::math::tools.
12 template <class F, class T, class Tol>
19 boost::uintmax_t& max_iter);
21 template <class F, class T, class Tol>
29 template <class F, class T, class Tol, class ``__Policy``>
36 boost::uintmax_t& max_iter,
39 // Bracket and Solve Root
40 template <class F, class T, class Tol>
42 bracket_and_solve_root(
48 boost::uintmax_t& max_iter);
50 template <class F, class T, class Tol, class ``__Policy``>
52 bracket_and_solve_root(
58 boost::uintmax_t& max_iter,
62 template <class F, class T, class Tol>
69 boost::uintmax_t& max_iter);
71 template <class F, class T, class Tol, class ``__Policy``>
78 boost::uintmax_t& max_iter,
81 template <class F, class T, class Tol>
90 boost::uintmax_t& max_iter);
92 template <class F, class T, class Tol, class ``__Policy``>
101 boost::uintmax_t& max_iter,
102 const ``__Policy``&);
104 // Termination conditions:
106 struct eps_tolerance;
110 struct equal_nearest_integer;
112 }}} // boost::math::tools namespaces
116 These functions solve the root of some function ['f(x)] -
117 ['without the need for any derivatives of ['f(x)]].
119 The `bracket_and_solve_root` functions use __root_finding_TOMS748
120 by Alefeld, Potra and Shi that is asymptotically the most efficient known,
121 and has been shown to be optimal for a certain classes of smooth functions.
122 Variants with and without __policy_section are provided.
124 Alternatively, __bisect is a simple __bisection_wikipedia routine which can be useful
125 in its own right in some situations, or alternatively for narrowing
126 down the range containing the root, prior to calling a more advanced
129 All the algorithms in this section reduce the diameter of the enclosing
130 interval with the same asymptotic efficiency with which they locate the
131 root. This is in contrast to the derivative based methods which may ['never]
132 significantly reduce the enclosing interval, even though they rapidly approach
133 the root. This is also in contrast to some other derivative-free methods
134 (for example, Brent's method described at
135 [@http://en.wikipedia.org/wiki/Brent%27s_method Brent-Dekker)]
136 which only reduces the enclosing interval on the final step.
137 Therefore these methods return a `std::pair` containing the enclosing interval found,
138 and accept a function object specifying the termination condition.
140 Three function objects are provided for ready-made termination conditions:
142 * ['eps_tolerance] causes termination when the relative error in the enclosing
143 interval is below a certain threshold.
144 * ['equal_floor] and ['equal_ceil] are useful for certain statistical applications
145 where the result is known to be an integer.
146 * Other user-defined termination conditions are likely to be used
147 only rarely, but may be useful in some specific circumstances.
149 [section:bisect Bisection]
151 template <class F, class T, class Tol>
153 bisect( // Unlimited iterations.
159 template <class F, class T, class Tol>
161 bisect( // Limited iterations.
166 boost::uintmax_t& max_iter);
168 template <class F, class T, class Tol, class ``__Policy``>
170 bisect( // Specified policy.
175 boost::uintmax_t& max_iter,
176 const ``__Policy``&);
178 These functions locate the root using __bisection_wikipedia.
180 `bisect` function arguments are:
183 [[f] [A unary functor (or C++ lambda) which is the function ['f(x)] whose root is to be found.]]
184 [[min] [The left bracket of the interval known to contain the root.]]
185 [[max] [The right bracket of the interval known to contain the root.[br]
186 It is a precondition that ['min < max] and ['f(min)*f(max) <= 0],
187 the function raises an __evaluation_error if these preconditions are violated.
188 The action taken on error is controlled by the __Policy template argument: the default behavior is to
189 throw a ['boost::math::evaluation_error]. If the __Policy is changed to not throw
190 then it returns ['std::pair<T>(min, min)].]]
191 [[tol] [A binary functor (or C++ lambda) that specifies the termination condition: the function
192 will return the current brackets enclosing the root when ['tol(min, max)] becomes true.
193 See also __root_termination.]]
194 [[max_iter][The maximum number of invocations of ['f(x)] to make while searching for the root. On exit, this is updated to the actual number of invocations performed.]]
199 [*Returns]: a pair of values ['r] that bracket the root so that:
201 [:f(r.first) * f(r.second) <= 0]
205 [:tol(r.first, r.second) == true]
211 where ['m] is the initial value of ['max_iter] passed to the function.
213 In other words, it's up to the caller to verify whether termination occurred
214 as a result of exceeding ['max_iter] function invocations (easily done by
215 checking the updated value of ['max_iter] when the function returns), rather than
216 because the termination condition ['tol] was satisfied.
218 [endsect] [/section:bisect Bisection]
220 [section:bracket_solve Bracket and Solve Root]
222 template <class F, class T, class Tol>
224 bracket_and_solve_root(
230 boost::uintmax_t& max_iter);
232 template <class F, class T, class Tol, class ``__Policy``>
234 bracket_and_solve_root(
240 boost::uintmax_t& max_iter,
241 const ``__Policy``&);
243 `bracket_and_solve_root` is a convenience function that calls __root_finding_TOMS748 internally
244 to find the root of ['f(x)]. It is generally much easier to use this function rather than __root_finding_TOMS748, since it
245 does the hard work of bracketing the root for you. It's bracketing routines are quite robust and will
246 usually be more foolproof than home-grown routines, unless the function can be analysed to yield tight
249 Note that this routine can only be used when:
251 * ['f(x)] is monotonic in the half of the real axis containing ['guess].
252 * The value of the inital guess must have the same sign as the root: the function
253 will ['never cross the origin] when searching for the root.
254 * The location of the root should be known at least approximately,
255 if the location of the root differs by many orders of magnitude
256 from ['guess] then many iterations will be needed to bracket the root in spite of
257 the special heuristics used to guard against this very situation. A typical example would be
258 setting the initial guess to 0.1, when the root is at 1e-300.
260 The `bracket_and_solve_root` parameters are:
263 [[f][A unary functor (or C++ lambda) that is the function whose root is to be solved.
264 ['f(x)] must be uniformly increasing or decreasing on ['x].]]
265 [[guess][An initial approximation to the root.]]
266 [[factor][A scaling factor that is used to bracket the root: the value
267 /guess/ is multiplied (or divided as appropriate) by /factor/
268 until two values are found that bracket the root. A value
269 such as 2 is a typical choice for ['factor].
270 In addition ['factor] will be multiplied by 2 every 32 iterations:
271 this is to guard against a really very bad initial guess, typically these occur
272 when it's known the result is very large or small, but not the exact order
274 [[rising][Set to ['true] if ['f(x)] is rising on /x/ and /false/ if ['f(x)]
275 is falling on /x/. This value is used along with the result
276 of /f(guess)/ to determine if /guess/ is
277 above or below the root.]]
278 [[tol] [A binary functor (or C++ lambda) that determines the termination condition for the search
279 for the root. /tol/ is passed the current brackets at each step,
280 when it returns true then the current brackets are returned as the pair result.
281 See also __root_termination.]]
282 [[max_iter] [The maximum number of function invocations to perform in the search
283 for the root. On exit is set to the actual number of invocations performed.]]
288 [*Returns]: a pair of values ['r] that bracket the root so that:
290 [:f(r.first) * f(r.second) <= 0]
294 [:tol(r.first, r.second) == true]
300 where ['m] is the initial value of ['max_iter] passed to the function.
302 In other words, it's up to the caller to verify whether termination occurred
303 as a result of exceeding ['max_iter] function invocations (easily done by
304 checking the value of ['max_iter] when the function returns), rather than
305 because the termination condition ['tol] was satisfied.
307 [endsect] [/section:bracket_solve Bracket and Solve Root]
309 [section:TOMS748 Algorithm TOMS 748: Alefeld, Potra and Shi: Enclosing zeros of continuous functions]
311 template <class F, class T, class Tol>
318 boost::uintmax_t& max_iter);
320 template <class F, class T, class Tol, class ``__Policy``>
327 boost::uintmax_t& max_iter,
328 const ``__Policy``&);
330 template <class F, class T, class Tol>
339 boost::uintmax_t& max_iter);
341 template <class F, class T, class Tol, class ``__Policy``>
350 boost::uintmax_t& max_iter,
351 const ``__Policy``&);
353 These functions implement TOMS Algorithm 748: it uses a mixture of
354 cubic, quadratic and linear (secant) interpolation to locate the root of
355 ['f(x)]. The two pairs of functions differ only by whether values for ['f(a)] and
356 ['f(b)] are already available.
358 Generally speaking it is easier (and often more efficient) to use __bracket_solve
359 rather than trying to bracket the root yourself as this function requires.
361 This function is provided rather than [@http://en.wikipedia.org/wiki/Brent%27s_method Brent's method] as it is known to be more
362 efficient in many cases (it is asymptotically the most efficient known,
363 and has been shown to be optimal for a certain classes of smooth functions).
364 It also has the useful property of decreasing the bracket size
365 with each step, unlike Brent's method which only shrinks the enclosing interval in the
366 final step. This makes it particularly useful when you need a result where the ends
367 of the interval round to the same integer: as often happens in statistical applications,
368 for example. In this situation the function is able to exit after a much smaller
369 number of iterations than would otherwise be possible.
371 The __root_finding_TOMS748 parameters are:
374 [[f] [A unary functor (or C++ lambda) that is the function whose root is to be solved.
375 f(x) need not be uniformly increasing or decreasing on ['x] and
376 may have multiple roots. However, the bounds given must bracket a single root.]]
377 [[a] [The lower bound for the initial bracket of the root.]]
378 [[b] [The upper bound for the initial bracket of the root.
379 It is a precondition that ['a < b] and that ['a] and ['b]
380 bracket the root to find so that ['f(a) * f(b) < 0].]]
381 [[fa] [Optional: the value of ['f(a)].]]
382 [[fb] [Optional: the value of ['f(b)].]]
383 [[tol] [A binary functor (or C++ lambda) that determines the termination condition for the search
384 for the root. ['tol] is passed the current brackets at each step,
385 when it returns true, then the current brackets are returned as the result.
386 See also __root_termination.]]
387 [[max_iter] [The maximum number of function invocations to perform in the search
388 for the root. On exit, ['max_iter] is set to actual number of function
394 `toms748_solve` returns: a pair of values ['r] that bracket the root so that:
396 [:['f(r.first) * f(r.second) <= 0]]
400 [:['tol(r.first, r.second) == true]]
406 where ['m] is the initial value of ['max_iter] passed to the function.
408 In other words, it's up to the caller to verify whether termination occurred
409 as a result of exceeding ['max_iter] function invocations (easily done by
410 checking the updated value of ['max_iter]
411 against its previous value passed as parameter),
412 rather than because the termination condition ['tol] was satisfied.
414 [endsect] [/section:TOMS748 Algorithm TOMS 748: Alefeld, Potra and Shi: Enclosing zeros of continuous functions]
416 [section:brent Brent-Decker Algorithm]
418 The [@http://en.wikipedia.org/wiki/Brent%27s_method Brent-Dekker algorithm], although very well know,
419 is not provided by this library as __root_finding_TOMS748 or
420 its slightly easier to use variant __bracket_solve are superior and provide equivalent functionality.
422 [endsect] [/section:brent Brent-Decker Algorithm]
424 [section:root_termination Termination Condition Functors]
430 eps_tolerance(int bits);
431 bool operator()(const T& a, const T& b)const;
434 `eps_tolerance` is the usual termination condition used with these root finding functions.
435 Its `operator()` will return true when the relative distance between ['a] and ['b]
436 is less than four times the machine epsilon for T, or 2[super 1-bits], whichever is
437 the larger. In other words, you set ['bits] to the number of bits of precision you
438 want in the result. The minimal tolerance of ['four times the machine epsilon of type T] is
439 required to ensure that we get back a bracketing interval, since this must clearly
440 be at greater than one epsilon in size. While in theory a maximum distance of twice
441 machine epsilon is possible to achieve, in practice this results in a great deal of "thrashing"
442 given that the function whose root is being found can only ever be accurate to 1 epsilon at best.
447 template <class T> bool operator()(const T& a, const T& b)const;
450 This termination condition is used when you want to find an integer result
451 that is the ['floor] of the true root. It will terminate as soon as both ends
452 of the interval have the same ['floor].
457 template <class T> bool operator()(const T& a, const T& b)const;
460 This termination condition is used when you want to find an integer result
461 that is the ['ceil] of the true root. It will terminate as soon as both ends
462 of the interval have the same ['ceil].
464 struct equal_nearest_integer
466 equal_nearest_integer();
467 template <class T> bool operator()(const T& a, const T& b)const;
470 This termination condition is used when you want to find an integer result
471 that is the /closest/ to the true root. It will terminate as soon as both ends
472 of the interval round to the same nearest integer.
474 [endsect] [/section:root_termination Termination Condition Functors]
476 [section:implementation Implementation]
478 The implementation of the bisection algorithm is extremely straightforward
479 and not detailed here.
481 __TOMS748 is described in detail in:
483 ['Algorithm 748: Enclosing Zeros of Continuous Functions,
484 G. E. Alefeld, F. A. Potra and Yixun Shi,
485 ACM Transactions on Mathematica1 Software, Vol. 21. No. 3. September 1995.
488 The implementation here is a faithful translation of this paper into C++.
490 [endsect] [/section:implementation Implementation]
492 [endsect] [/section:roots_noderiv Root Finding Without Derivatives]
495 Copyright 2006, 2010, 2015 John Maddock and Paul A. Bristow.
496 Distributed under the Boost Software License, Version 1.0.
497 (See accompanying file LICENSE_1_0.txt or copy at
498 http://www.boost.org/LICENSE_1_0.txt).