math/doc/roots/roots_without_derivatives.qbk

[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).
]