123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499 |
- [section:roots_noderiv Root Finding Without Derivatives]
- [h4 Synopsis]
- ``
- #include <boost/math/tools/roots.hpp>
- ``
- namespace boost { namespace math {
- namespace tools { // Note namespace boost::math::tools.
- // Bisection
- template <class F, class T, class Tol>
- std::pair<T, T>
- bisect(
- F f,
- T min,
- T max,
- Tol tol,
- boost::uintmax_t& max_iter);
- template <class F, class T, class Tol>
- std::pair<T, T>
- bisect(
- F f,
- T min,
- T max,
- Tol tol);
- template <class F, class T, class Tol, class ``__Policy``>
- std::pair<T, T>
- bisect(
- F f,
- T min,
- T max,
- Tol tol,
- boost::uintmax_t& max_iter,
- const ``__Policy``&);
- // Bracket and Solve Root
- template <class F, class T, class Tol>
- std::pair<T, T>
- bracket_and_solve_root(
- F f,
- const T& guess,
- const T& factor,
- bool rising,
- Tol tol,
- boost::uintmax_t& max_iter);
- template <class F, class T, class Tol, class ``__Policy``>
- std::pair<T, T>
- bracket_and_solve_root(
- F f,
- const T& guess,
- const T& factor,
- bool rising,
- Tol tol,
- boost::uintmax_t& max_iter,
- const ``__Policy``&);
- // TOMS 748 algorithm
- template <class F, class T, class Tol>
- std::pair<T, T>
- toms748_solve(
- F f,
- const T& a,
- const T& b,
- Tol tol,
- boost::uintmax_t& max_iter);
- template <class F, class T, class Tol, class ``__Policy``>
- std::pair<T, T>
- toms748_solve(
- F f,
- const T& a,
- const T& b,
- Tol tol,
- boost::uintmax_t& max_iter,
- const ``__Policy``&);
- template <class F, class T, class Tol>
- std::pair<T, T>
- toms748_solve(
- F f,
- const T& a,
- const T& b,
- const T& fa,
- const T& fb,
- Tol tol,
- boost::uintmax_t& max_iter);
- template <class F, class T, class Tol, class ``__Policy``>
- std::pair<T, T>
- toms748_solve(
- F f,
- const T& a,
- const T& b,
- const T& fa,
- const T& fb,
- Tol tol,
- boost::uintmax_t& max_iter,
- const ``__Policy``&);
- // Termination conditions:
- template <class T>
- struct eps_tolerance;
- struct equal_floor;
- struct equal_ceil;
- struct equal_nearest_integer;
- }}} // boost::math::tools namespaces
- [h4 Description]
- These functions solve the root of some function ['f(x)] -
- ['without the need for any derivatives of ['f(x)]].
- The `bracket_and_solve_root` functions use __root_finding_TOMS748
- by Alefeld, Potra and Shi that is asymptotically the most efficient known,
- and has been shown to be optimal for a certain classes of smooth functions.
- Variants with and without __policy_section are provided.
- Alternatively, __bisect is a simple __bisection_wikipedia routine which can be useful
- in its own right in some situations, or alternatively for narrowing
- down the range containing the root, prior to calling a more advanced
- algorithm.
- All the algorithms in this section reduce the diameter of the enclosing
- interval with the same asymptotic efficiency with which they locate the
- root. This is in contrast to the derivative based methods which may ['never]
- significantly reduce the enclosing interval, even though they rapidly approach
- the root. This is also in contrast to some other derivative-free methods
- (for example, Brent's method described at
- [@http://en.wikipedia.org/wiki/Brent%27s_method Brent-Dekker)]
- which only reduces the enclosing interval on the final step.
- Therefore these methods return a `std::pair` containing the enclosing interval found,
- and accept a function object specifying the termination condition.
- Three function objects are provided for ready-made termination conditions:
- * ['eps_tolerance] causes termination when the relative error in the enclosing
- interval is below a certain threshold.
- * ['equal_floor] and ['equal_ceil] are useful for certain statistical applications
- where the result is known to be an integer.
- * Other user-defined termination conditions are likely to be used
- only rarely, but may be useful in some specific circumstances.
- [section:bisect Bisection]
- template <class F, class T, class Tol>
- std::pair<T, T>
- bisect( // Unlimited iterations.
- F f,
- T min,
- T max,
- Tol tol);
- template <class F, class T, class Tol>
- std::pair<T, T>
- bisect( // Limited iterations.
- F f,
- T min,
- T max,
- Tol tol,
- boost::uintmax_t& max_iter);
- template <class F, class T, class Tol, class ``__Policy``>
- std::pair<T, T>
- bisect( // Specified policy.
- F f,
- T min,
- T max,
- Tol tol,
- boost::uintmax_t& max_iter,
- const ``__Policy``&);
- These functions locate the root using __bisection_wikipedia.
- `bisect` function arguments are:
- [variablelist
- [[f] [A unary functor (or C++ lambda) which is the function ['f(x)] whose root is to be found.]]
- [[min] [The left bracket of the interval known to contain the root.]]
- [[max] [The right bracket of the interval known to contain the root.[br]
- It is a precondition that ['min < max] and ['f(min)*f(max) <= 0],
- the function raises an __evaluation_error if these preconditions are violated.
- The action taken on error is controlled by the __Policy template argument: the default behavior is to
- throw a ['boost::math::evaluation_error]. If the __Policy is changed to not throw
- then it returns ['std::pair<T>(min, min)].]]
- [[tol] [A binary functor (or C++ lambda) that specifies the termination condition: the function
- will return the current brackets enclosing the root when ['tol(min, max)] becomes true.
- See also __root_termination.]]
- [[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.]]
- ]
- [optional_policy]
- [*Returns]: a pair of values ['r] that bracket the root so that:
- [:f(r.first) * f(r.second) <= 0]
- and either
- [:tol(r.first, r.second) == true]
- or
- [:max_iter >= m]
- where ['m] is the initial value of ['max_iter] passed to the function.
- In other words, it's up to the caller to verify whether termination occurred
- as a result of exceeding ['max_iter] function invocations (easily done by
- checking the updated value of ['max_iter] when the function returns), rather than
- because the termination condition ['tol] was satisfied.
- [endsect] [/section:bisect Bisection]
- [section:bracket_solve Bracket and Solve Root]
- template <class F, class T, class Tol>
- std::pair<T, T>
- bracket_and_solve_root(
- F f,
- const T& guess,
- const T& factor,
- bool rising,
- Tol tol,
- boost::uintmax_t& max_iter);
- template <class F, class T, class Tol, class ``__Policy``>
- std::pair<T, T>
- bracket_and_solve_root(
- F f,
- const T& guess,
- const T& factor,
- bool rising,
- Tol tol,
- boost::uintmax_t& max_iter,
- const ``__Policy``&);
- `bracket_and_solve_root` is a convenience function that calls __root_finding_TOMS748 internally
- to find the root of ['f(x)]. It is generally much easier to use this function rather than __root_finding_TOMS748, since it
- does the hard work of bracketing the root for you. It's bracketing routines are quite robust and will
- usually be more foolproof than home-grown routines, unless the function can be analysed to yield tight
- brackets.
- Note that this routine can only be used when:
- * ['f(x)] is monotonic in the half of the real axis containing ['guess].
- * The value of the inital guess must have the same sign as the root: the function
- will ['never cross the origin] when searching for the root.
- * The location of the root should be known at least approximately,
- if the location of the root differs by many orders of magnitude
- from ['guess] then many iterations will be needed to bracket the root in spite of
- the special heuristics used to guard against this very situation. A typical example would be
- setting the initial guess to 0.1, when the root is at 1e-300.
- The `bracket_and_solve_root` parameters are:
- [variablelist
- [[f][A unary functor (or C++ lambda) that is the function whose root is to be solved.
- ['f(x)] must be uniformly increasing or decreasing on ['x].]]
- [[guess][An initial approximation to the root.]]
- [[factor][A scaling factor that is used to bracket the root: the value
- /guess/ is multiplied (or divided as appropriate) by /factor/
- until two values are found that bracket the root. A value
- such as 2 is a typical choice for ['factor].
- In addition ['factor] will be multiplied by 2 every 32 iterations:
- this is to guard against a really very bad initial guess, typically these occur
- when it's known the result is very large or small, but not the exact order
- of magnitude.]]
- [[rising][Set to ['true] if ['f(x)] is rising on /x/ and /false/ if ['f(x)]
- is falling on /x/. This value is used along with the result
- of /f(guess)/ to determine if /guess/ is
- above or below the root.]]
- [[tol] [A binary functor (or C++ lambda) that determines the termination condition for the search
- for the root. /tol/ is passed the current brackets at each step,
- when it returns true then the current brackets are returned as the pair result.
- See also __root_termination.]]
- [[max_iter] [The maximum number of function invocations to perform in the search
- for the root. On exit is set to the actual number of invocations performed.]]
- ]
- [optional_policy]
- [*Returns]: a pair of values ['r] that bracket the root so that:
- [:f(r.first) * f(r.second) <= 0]
- and either
- [:tol(r.first, r.second) == true]
- or
- [:max_iter >= m]
- where ['m] is the initial value of ['max_iter] passed to the function.
- In other words, it's up to the caller to verify whether termination occurred
- as a result of exceeding ['max_iter] function invocations (easily done by
- checking the value of ['max_iter] when the function returns), rather than
- because the termination condition ['tol] was satisfied.
- [endsect] [/section:bracket_solve Bracket and Solve Root]
- [section:TOMS748 Algorithm TOMS 748: Alefeld, Potra and Shi: Enclosing zeros of continuous functions]
- template <class F, class T, class Tol>
- std::pair<T, T>
- toms748_solve(
- F f,
- const T& a,
- const T& b,
- Tol tol,
- boost::uintmax_t& max_iter);
- template <class F, class T, class Tol, class ``__Policy``>
- std::pair<T, T>
- toms748_solve(
- F f,
- const T& a,
- const T& b,
- Tol tol,
- boost::uintmax_t& max_iter,
- const ``__Policy``&);
- template <class F, class T, class Tol>
- std::pair<T, T>
- toms748_solve(
- F f,
- const T& a,
- const T& b,
- const T& fa,
- const T& fb,
- Tol tol,
- boost::uintmax_t& max_iter);
- template <class F, class T, class Tol, class ``__Policy``>
- std::pair<T, T>
- toms748_solve(
- F f,
- const T& a,
- const T& b,
- const T& fa,
- const T& fb,
- Tol tol,
- boost::uintmax_t& max_iter,
- const ``__Policy``&);
- These functions implement TOMS Algorithm 748: it uses a mixture of
- cubic, quadratic and linear (secant) interpolation to locate the root of
- ['f(x)]. The two pairs of functions differ only by whether values for ['f(a)] and
- ['f(b)] are already available.
- Generally speaking it is easier (and often more efficient) to use __bracket_solve
- rather than trying to bracket the root yourself as this function requires.
- This function is provided rather than [@http://en.wikipedia.org/wiki/Brent%27s_method Brent's method] as it is known to be more
- efficient in many cases (it is asymptotically the most efficient known,
- and has been shown to be optimal for a certain classes of smooth functions).
- It also has the useful property of decreasing the bracket size
- with each step, unlike Brent's method which only shrinks the enclosing interval in the
- final step. This makes it particularly useful when you need a result where the ends
- of the interval round to the same integer: as often happens in statistical applications,
- for example. In this situation the function is able to exit after a much smaller
- number of iterations than would otherwise be possible.
- The __root_finding_TOMS748 parameters are:
- [variablelist
- [[f] [A unary functor (or C++ lambda) that is the function whose root is to be solved.
- f(x) need not be uniformly increasing or decreasing on ['x] and
- may have multiple roots. However, the bounds given must bracket a single root.]]
- [[a] [The lower bound for the initial bracket of the root.]]
- [[b] [The upper bound for the initial bracket of the root.
- It is a precondition that ['a < b] and that ['a] and ['b]
- bracket the root to find so that ['f(a) * f(b) < 0].]]
- [[fa] [Optional: the value of ['f(a)].]]
- [[fb] [Optional: the value of ['f(b)].]]
- [[tol] [A binary functor (or C++ lambda) that determines the termination condition for the search
- for the root. ['tol] is passed the current brackets at each step,
- when it returns true, then the current brackets are returned as the result.
- See also __root_termination.]]
- [[max_iter] [The maximum number of function invocations to perform in the search
- for the root. On exit, ['max_iter] is set to actual number of function
- invocations used.]]
- ]
- [optional_policy]
- `toms748_solve` returns: a pair of values ['r] that bracket the root so that:
- [:['f(r.first) * f(r.second) <= 0]]
- and either
- [:['tol(r.first, r.second) == true]]
- or
- [:['max_iter >= m]]
- where ['m] is the initial value of ['max_iter] passed to the function.
- In other words, it's up to the caller to verify whether termination occurred
- as a result of exceeding ['max_iter] function invocations (easily done by
- checking the updated value of ['max_iter]
- against its previous value passed as parameter),
- rather than because the termination condition ['tol] was satisfied.
- [endsect] [/section:TOMS748 Algorithm TOMS 748: Alefeld, Potra and Shi: Enclosing zeros of continuous functions]
- [section:brent Brent-Decker Algorithm]
- The [@http://en.wikipedia.org/wiki/Brent%27s_method Brent-Dekker algorithm], although very well know,
- is not provided by this library as __root_finding_TOMS748 or
- its slightly easier to use variant __bracket_solve are superior and provide equivalent functionality.
- [endsect] [/section:brent Brent-Decker Algorithm]
- [section:root_termination Termination Condition Functors]
- template <class T>
- struct eps_tolerance
- {
- eps_tolerance();
- eps_tolerance(int bits);
- bool operator()(const T& a, const T& b)const;
- };
- `eps_tolerance` is the usual termination condition used with these root finding functions.
- Its `operator()` will return true when the relative distance between ['a] and ['b]
- is less than four times the machine epsilon for T, or 2[super 1-bits], whichever is
- the larger. In other words, you set ['bits] to the number of bits of precision you
- want in the result. The minimal tolerance of ['four times the machine epsilon of type T] is
- required to ensure that we get back a bracketing interval, since this must clearly
- be at greater than one epsilon in size. While in theory a maximum distance of twice
- machine epsilon is possible to achieve, in practice this results in a great deal of "thrashing"
- given that the function whose root is being found can only ever be accurate to 1 epsilon at best.
- struct equal_floor
- {
- equal_floor();
- template <class T> bool operator()(const T& a, const T& b)const;
- };
- This termination condition is used when you want to find an integer result
- that is the ['floor] of the true root. It will terminate as soon as both ends
- of the interval have the same ['floor].
- struct equal_ceil
- {
- equal_ceil();
- template <class T> bool operator()(const T& a, const T& b)const;
- };
- This termination condition is used when you want to find an integer result
- that is the ['ceil] of the true root. It will terminate as soon as both ends
- of the interval have the same ['ceil].
- struct equal_nearest_integer
- {
- equal_nearest_integer();
- template <class T> bool operator()(const T& a, const T& b)const;
- };
- This termination condition is used when you want to find an integer result
- that is the /closest/ to the true root. It will terminate as soon as both ends
- of the interval round to the same nearest integer.
- [endsect] [/section:root_termination Termination Condition Functors]
- [section:implementation Implementation]
- The implementation of the bisection algorithm is extremely straightforward
- and not detailed here.
- __TOMS748 is described in detail in:
- ['Algorithm 748: Enclosing Zeros of Continuous Functions,
- G. E. Alefeld, F. A. Potra and Yixun Shi,
- ACM Transactions on Mathematica1 Software, Vol. 21. No. 3. September 1995.
- Pages 327-344.]
- The implementation here is a faithful translation of this paper into C++.
- [endsect] [/section:implementation Implementation]
- [endsect] [/section:roots_noderiv Root Finding Without Derivatives]
- [/
- Copyright 2006, 2010, 2015 John Maddock and Paul A. Bristow.
- Distributed under the Boost Software License, Version 1.0.
- (See accompanying file LICENSE_1_0.txt or copy at
- http://www.boost.org/LICENSE_1_0.txt).
- ]
|