...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
#include <boost/math/tools/roots.hpp>
namespace boost{ namespace math{ namespace 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; }}} // namespaces
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 TOMS Algorithm 748 that is asymptotically the most efficient
known, and have been shown to be optimal for a certain classes of smooth
functions. Variants with and without __policies are provided.
Alternatively, there is a simple bisection 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 the methods of Brent or Dekker) which only reduce 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, while 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.
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&);
These functions locate the root using bisection: function arguments are:
A unary functor which is the function whose root is to be found.
The left bracket of the interval known to contain the root.
The right bracket of the interval known to contain the root. It is a precondition that min < max and f(min)*f(max) <= 0, the function signals evaluation error if these preconditions are violated. The action taken is controlled by the evaluation error policy. A best guess may be returned, perhaps significantly wrong.
A binary functor that specifies the termination condition: the function will return the current brackets enclosing the root when tol(min,max) becomes true.
The maximum number of invocations of f(x) to make while searching for the root.
The final Policy argument is optional and can be used to control the behaviour of the function: how it handles errors, what level of precision to use etc. Refer to the policy documentation for more details.
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.
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&);
This is a convenience function that calls toms748_solve internally to find the root of f(x). It's usable only when f(x) is a monotonic function, and the location of the root is known approximately, and in particular it is known whether the root is occurs for positive or negative x. The parameters are:
A unary functor that is the function whose root is to be solved. f(x) must be uniformly increasing or decreasing on x.
An initial approximation to the root
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.
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.
A binary functor 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.
The maximum number of function invocations to perform in the search for the root.
The final Policy argument is optional and can be used to control the behaviour of the function: how it handles errors, what level of precision to use etc. Refer to the policy documentation for more details.
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.
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 two 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 functions differ only by whether values for f(a) and f(b) are already available. The toms748_solve parameters are:
A unary functor 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.
The lower bound for the initial bracket of the root.
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.
Optional: the value of f(a).
Optional: the value of f(b).
A binary functor 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.
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.
The final Policy argument is optional and can be used to control the behaviour of the function: how it handles errors, what level of precision to use etc. Refer to the policy documentation for more details.
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.
template <class T> struct 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 twice the machine epsilon for T,
or 21-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 twice the machine epsilon of T is required to ensure that we get back
a bracketing interval: since this must clearly be at least 1 epsilon in size.
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.
The implementation of the bisection algorithm is extremely straightforward and not detailed here. TOMS Algorithm 748: enclosing zeros of continuous functions 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++.