libs/math/doc/roots/roots.qbk

   1 [section:roots_deriv Root Finding With Derivatives: Newton-Raphson, Halley & Schr'''&#xf6;'''der]
   2
   3 [h4 Synopsis]
   4
   5 ``
   6 #include <boost/math/tools/roots.hpp>
   7 ``
   8
   9    namespace boost { namespace math {
  10    namespace tools { // Note namespace boost::math::tools.
  11    // Newton-Raphson
  12    template <class F, class T>
  13    T newton_raphson_iterate(F f, T guess, T min, T max, int digits);
  14
  15    template <class F, class T>
  16    T newton_raphson_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter);
  17
  18    // Halley
  19    template <class F, class T>
  20    T halley_iterate(F f, T guess, T min, T max, int digits);
  21
  22    template <class F, class T>
  23    T halley_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter);
  24
  25    // Schr'''&#xf6;'''der
  26    template <class F, class T>
  27    T schroder_iterate(F f, T guess, T min, T max, int digits);
  28
  29    template <class F, class T>
  30    T schroder_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter);
  31
  32    template <class F, class Complex>
  33    Complex complex_newton(F f, Complex guess, int max_iterations = std::numeric_limits<typename Complex::value_type>::digits);
  34
  35    template<class T>
  36    auto quadratic_roots(T const & a, T const & b, T const & c);
  37
  38    }}} // namespaces boost::math::tools.
  39
  40 [h4 Description]
  41
  42 These functions all perform iterative root-finding [*using derivatives]:
  43
  44 * `newton_raphson_iterate` performs second-order __newton.
  45
  46 * `halley_iterate` and `schroder_iterate` perform third-order
  47 __halley and __schroder iteration.
  48
  49 * `complex_newton` performs Newton's method on complex-analytic functions.
  50
  51 * `solve_quadratic` solves quadratic equations using various tricks to keep catastrophic cancellation from occurring in computation of the discriminant.
  52
  53
  54 [variablelist Parameters of the real-valued root finding functions
  55 [[F f] [Type F must be a callable function object (or C++ lambda) that accepts one parameter and
  56         returns a __tuple_type:
  57
  58 For second-order iterative method ([@http://en.wikipedia.org/wiki/Newton_Raphson Newton Raphson])
  59         the `tuple` should have [*two] elements containing the evaluation
  60         of the function and its first derivative.
  61
  62 For the third-order methods
  63 ([@http://en.wikipedia.org/wiki/Halley%27s_method Halley] and
  64 Schr'''&#xf6;'''der)
  65         the `tuple` should have [*three] elements containing the evaluation of
  66         the function and its first and second derivatives.]]
  67 [[T guess] [The initial starting value. A good guess is crucial to quick convergence!]]
  68 [[T min] [The minimum possible value for the result, this is used as an initial lower bracket.]]
  69 [[T max] [The maximum possible value for the result, this is used as an initial upper bracket.]]
  70 [[int digits] [The desired number of binary digits precision.]]
  71 [[uintmax_t& max_iter] [An optional maximum number of iterations to perform.  On exit, this is updated to the actual number of iterations performed.]]
  72 ]
  73
  74 When using these functions you should note that:
  75
  76 * Default `max_iter = (std::numeric_limits<boost::uintmax_t>::max)()` is effectively 'iterate for ever'.
  77 * They may be very sensitive to the initial guess, typically they converge very rapidly
  78 if the initial guess has two or three decimal digits correct.  However convergence
  79 can be no better than __bisect, or in some rare cases, even worse than __bisect if the
  80 initial guess is a long way from the correct value and the derivatives are close to zero.
  81 * These functions include special cases to handle zero first (and second where appropriate)
  82 derivatives, and fall back to __bisect in this case.  However, it is helpful
  83 if functor F is defined to return an arbitrarily small value ['of the correct sign] rather
  84 than zero.
  85 * The functions will raise an __evaluation_error if arguments `min` and `max` are the wrong way around
  86 or if they converge to a local minima.
  87 * If the derivative at the current best guess for the result is infinite (or
  88 very close to being infinite) then these functions may terminate prematurely.
  89 A large first derivative leads to a very small next step, triggering the termination
  90 condition.  Derivative based iteration may not be appropriate in such cases.
  91 * If the function is 'Really Well Behaved' (is monotonic and has only one root)
  92 the bracket bounds ['min] and ['max] may as well be set to the widest limits
  93 like zero and `numeric_limits<T>::max()`.
  94 *But if the function more complex and may have more than one root or a pole,
  95 the choice of bounds is protection against jumping out to seek the 'wrong' root.
  96 * These functions fall back to __bisect if the next computed step would take the
  97 next value out of bounds.  The bounds are updated after each step to ensure this leads
  98 to convergence.  However, a good initial guess backed up by asymptotically-tight
  99 bounds will improve performance no end - rather than relying on __bisection.
 100 * The value of ['digits] is crucial to good performance of these functions,
 101 if it is set too high then at best you will get one extra (unnecessary)
 102 iteration, and at worst the last few steps will proceed by __bisection.
 103 Remember that the returned value can never be more accurate than ['f(x)] can be
 104 evaluated, and that if ['f(x)] suffers from cancellation errors as it
 105 tends to zero then the computed steps will be effectively random.  The
 106 value of ['digits] should be set so that iteration terminates before this point:
 107 remember that for second and third order methods the number of correct
 108 digits in the result is increasing quite
 109 substantially with each iteration, ['digits] should be set by experiment so that the final
 110 iteration just takes the next value into the zone where ['f(x)] becomes inaccurate.
 111 A good starting point for ['digits] would be 0.6*D for Newton and 0.4*D for Halley or Shr'''&#xf6;'''der
 112 iteration, where D is `std::numeric_limits<T>::digits`.
 113 * If you need some diagnostic output to see what is going on, you can
 114 `#define BOOST_MATH_INSTRUMENT` before the `#include <boost/math/tools/roots.hpp>`,
 115 and also ensure that display of all the significant digits with
 116 ` cout.precision(std::numeric_limits<double>::digits10)`:
 117 or even possibly significant digits with
 118 ` cout.precision(std::numeric_limits<double>::max_digits10)`:
 119 but be warned, this may produce copious output!
 120 * Finally: you may well be able to do better than these functions by hand-coding
 121 the heuristics used so that they are tailored to a specific function.  You may also
 122 be able to compute the ratio of derivatives used by these methods more efficiently
 123 than computing the derivatives themselves.  As ever, algebraic simplification can
 124 be a big win.
 125
 126 [h4:newton Newton Raphson Method]
 127
 128 Given an initial guess ['x0] the subsequent values are computed using:
 129
 130 [equation roots1]
 131
 132 Out-of-bounds steps revert to __bisection of the current bounds.
 133
 134 Under ideal conditions, the number of correct digits doubles with each iteration.
 135
 136 [h4:halley Halley's Method]
 137
 138 Given an initial guess ['x0] the subsequent values are computed using:
 139
 140 [equation roots2]
 141
 142 Over-compensation by the second derivative (one which would proceed
 143 in the wrong direction) causes the method to
 144 revert to a Newton-Raphson step.
 145
 146 Out of bounds steps revert to bisection of the current bounds.
 147
 148 Under ideal conditions, the number of correct digits trebles with each iteration.
 149
 150 [h4:schroder Schr'''&#xf6;'''der's Method]
 151
 152 Given an initial guess x0 the subsequent values are computed using:
 153
 154 [equation roots3]
 155
 156 Over-compensation by the second derivative (one which would proceed
 157 in the wrong direction) causes the method to
 158 revert to a Newton-Raphson step.  Likewise a Newton step is used
 159 whenever that Newton step would change the next value by more than 10%.
 160
 161 Out of bounds steps revert to __bisection_wikipedia of the current bounds.
 162
 163 Under ideal conditions, the number of correct digits trebles with each iteration.
 164
 165 This is Schr'''&#xf6;'''der's general result (equation 18 from [@http://drum.lib.umd.edu/handle/1903/577 Stewart, G. W.
 166 "On Infinitely Many Algorithms for Solving Equations." English translation of Schr'''&#xf6;'''der's original paper.
 167 College Park, MD: University of Maryland, Institute for Advanced Computer Studies, Department of Computer Science, 1993].)
 168
 169 This method guarantees at least quadratic convergence (the same as Newton's method), and is known to work well in the presence of multiple roots:
 170 something that neither Newton nor Halley can do.
 171
 172 The complex Newton method works slightly differently than the rest of the methods:
 173 Since there is no way to bracket roots in the complex plane,
 174 the `min` and `max` arguments are not accepted.
 175 Failure to reach a root is communicated by returning `nan`s.
 176 Remember that if a function has many roots,
 177 then which root the complex Newton's method converges to is essentially impossible to predict a priori; see the Newton's fractal for more information.
 178
 179 Finally, the derivative of /f/ must be continuous at the root or else non-roots can be found; see [@https://math.stackexchange.com/questions/3017766/constructing-newton-iteration-converging-to-non-root here] for an example.
 180
 181 An example usage of `complex_newton` is given in `examples/daubechies_coefficients.cpp`.
 182
 183 [h4 Quadratics]
 184
 185 To solve a quadratic /ax/[super 2] + /bx/ + /c/ = 0, we may use
 186
 187     auto [x0, x1] = boost::math::tools::quadratic_roots(a, b, c);
 188
 189 If the roots are real, they are arranged so that `x0` \u2264 `x1`.
 190 If the roots are complex and the inputs are real, `x0` and `x1` are both `std::numeric_limits<Real>::quiet_NaN()`.
 191 In this case we must cast `a`, `b` and `c` to a complex type to extract the complex roots.
 192 If `a`, `b` and `c` are integral, then the roots are of type double.
 193 The routine is much faster if the fused-multiply-add instruction is available on your architecture.
 194 If the fma is not available, the function resorts to slow emulation.
 195 Finally, speed is improved if you compile for your particular architecture.
 196 For instance, if you compile without any architecture flags, then the `std::fma` call compiles down to `call _fma`,
 197 which dynamically chooses to emulate or execute the `vfmadd132sd` instruction based on the capabilities of the architecture.
 198 If instead, you compile with (say) `-march=native` then no dynamic choice is made:
 199 The `vfmadd132sd` instruction is always executed if available and emulation is used if not.
 200
 201
 202 [h4 Examples]
 203
 204 See __root_finding_examples.
 205
 206 [endsect] [/section:roots_deriv Root Finding With Derivatives]
 207
 208 [/
 209   Copyright 2006, 2010, 2012 John Maddock and Paul A. Bristow.
 210   Distributed under the Boost Software License, Version 1.0.
 211   (See accompanying file LICENSE_1_0.txt or copy at
 212   http://www.boost.org/LICENSE_1_0.txt).
 213 ]