ceph/src/boost/libs/math/doc/html/math_toolkit/roots/roots_deriv.html

   1 <html>
   2 <head>
   3 <meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
   4 <title>Root Finding With Derivatives: Newton-Raphson, Halley &amp; Schr&#246;der</title>
   5 <link rel="stylesheet" href="../../math.css" type="text/css">
   6 <meta name="generator" content="DocBook XSL Stylesheets V1.77.1">
   7 <link rel="home" href="../../index.html" title="Math Toolkit 2.5.1">
   8 <link rel="up" href="../roots.html" title="Root finding">
   9 <link rel="prev" href="roots_noderiv/implementation.html" title="Implementation">
  10 <link rel="next" href="root_finding_examples.html" title="Examples of Root-Finding (with and without derivatives)">
  11 </head>
  12 <body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
  13 <table cellpadding="2" width="100%"><tr>
  14 <td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../../boost.png"></td>
  15 <td align="center"><a href="../../../../../../index.html">Home</a></td>
  16 <td align="center"><a href="../../../../../../libs/libraries.htm">Libraries</a></td>
  17 <td align="center"><a href="http://www.boost.org/users/people.html">People</a></td>
  18 <td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td>
  19 <td align="center"><a href="../../../../../../more/index.htm">More</a></td>
  20 </tr></table>
  21 <hr>
  22 <div class="spirit-nav">
  23 <a accesskey="p" href="roots_noderiv/implementation.html"><img src="../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../roots.html"><img src="../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="root_finding_examples.html"><img src="../../../../../../doc/src/images/next.png" alt="Next"></a>
  24 </div>
  25 <div class="section">
  26 <div class="titlepage"><div><div><h3 class="title">
  27 <a name="math_toolkit.roots.roots_deriv"></a><a class="link" href="roots_deriv.html" title="Root Finding With Derivatives: Newton-Raphson, Halley &amp; Schr&#246;der">Root Finding With Derivatives:
  28       Newton-Raphson, Halley &amp; Schr&#246;der</a>
  29 </h3></div></div></div>
  30 <h5>
  31 <a name="math_toolkit.roots.roots_deriv.h0"></a>
  32         <span class="phrase"><a name="math_toolkit.roots.roots_deriv.synopsis"></a></span><a class="link" href="roots_deriv.html#math_toolkit.roots.roots_deriv.synopsis">Synopsis</a>
  33       </h5>
  34 <pre class="programlisting"><span class="preprocessor">#include</span> <span class="special">&lt;</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">tools</span><span class="special">/</span><span class="identifier">roots</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">&gt;</span>
  35 </pre>
  36 <pre class="programlisting"><span class="keyword">namespace</span> <span class="identifier">boost</span> <span class="special">{</span> <span class="keyword">namespace</span> <span class="identifier">math</span> <span class="special">{</span>
  37 <span class="keyword">namespace</span> <span class="identifier">tools</span> <span class="special">{</span> <span class="comment">// Note namespace boost::math::tools.</span>
  38 <span class="comment">// Newton-Raphson</span>
  39 <span class="keyword">template</span> <span class="special">&lt;</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">T</span><span class="special">&gt;</span>
  40 <span class="identifier">T</span> <span class="identifier">newton_raphson_iterate</span><span class="special">(</span><span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">guess</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">min</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">max</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">digits</span><span class="special">);</span>
  41
  42 <span class="keyword">template</span> <span class="special">&lt;</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">T</span><span class="special">&gt;</span>
  43 <span class="identifier">T</span> <span class="identifier">newton_raphson_iterate</span><span class="special">(</span><span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">guess</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">min</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">max</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">digits</span><span class="special">,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span><span class="special">&amp;</span> <span class="identifier">max_iter</span><span class="special">);</span>
  44
  45 <span class="comment">// Halley</span>
  46 <span class="keyword">template</span> <span class="special">&lt;</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">T</span><span class="special">&gt;</span>
  47 <span class="identifier">T</span> <span class="identifier">halley_iterate</span><span class="special">(</span><span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">guess</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">min</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">max</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">digits</span><span class="special">);</span>
  48
  49 <span class="keyword">template</span> <span class="special">&lt;</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">T</span><span class="special">&gt;</span>
  50 <span class="identifier">T</span> <span class="identifier">halley_iterate</span><span class="special">(</span><span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">guess</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">min</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">max</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">digits</span><span class="special">,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span><span class="special">&amp;</span> <span class="identifier">max_iter</span><span class="special">);</span>
  51
  52 <span class="comment">// Schr'''&amp;#xf6;'''der</span>
  53 <span class="keyword">template</span> <span class="special">&lt;</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">T</span><span class="special">&gt;</span>
  54 <span class="identifier">T</span> <span class="identifier">schroder_iterate</span><span class="special">(</span><span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">guess</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">min</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">max</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">digits</span><span class="special">);</span>
  55
  56 <span class="keyword">template</span> <span class="special">&lt;</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">T</span><span class="special">&gt;</span>
  57 <span class="identifier">T</span> <span class="identifier">schroder_iterate</span><span class="special">(</span><span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">guess</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">min</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">max</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">digits</span><span class="special">,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span><span class="special">&amp;</span> <span class="identifier">max_iter</span><span class="special">);</span>
  58
  59 <span class="special">}}}</span> <span class="comment">// namespaces boost::math::tools.</span>
  60 </pre>
  61 <h5>
  62 <a name="math_toolkit.roots.roots_deriv.h1"></a>
  63         <span class="phrase"><a name="math_toolkit.roots.roots_deriv.description"></a></span><a class="link" href="roots_deriv.html#math_toolkit.roots.roots_deriv.description">Description</a>
  64       </h5>
  65 <p>
  66         These functions all perform iterative root-finding <span class="bold"><strong>using
  67         derivatives</strong></span>:
  68       </p>
  69 <div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; ">
  70 <li class="listitem">
  71             <code class="computeroutput"><span class="identifier">newton_raphson_iterate</span></code>
  72             performs second-order <a class="link" href="roots_deriv.html#math_toolkit.roots.roots_deriv.newton">Newton-Raphson
  73             iteration</a>.
  74           </li>
  75 <li class="listitem">
  76             <code class="computeroutput"><span class="identifier">halley_iterate</span></code> and <code class="computeroutput"><span class="identifier">schroder_iterate</span></code> perform third-order
  77             <a class="link" href="roots_deriv.html#math_toolkit.roots.roots_deriv.halley">Halley</a> and
  78             <a class="link" href="roots_deriv.html#math_toolkit.roots.roots_deriv.schroder">Schr&#246;der</a>
  79             iteration.
  80           </li>
  81 </ul></div>
  82 <p>
  83         The functions all take the same parameters:
  84       </p>
  85 <div class="variablelist">
  86 <p class="title"><b>Parameters of the root finding functions</b></p>
  87 <dl class="variablelist">
  88 <dt><span class="term">F f</span></dt>
  89 <dd>
  90 <p>
  91               Type F must be a callable function object that accepts one parameter
  92               and returns a <a class="link" href="../internals/tuples.html" title="Tuples">std::pair,
  93               std::tuple, boost::tuple or boost::fusion::tuple</a>:
  94             </p>
  95 <p>
  96               For second-order iterative method (<a href="http://en.wikipedia.org/wiki/Newton_Raphson" target="_top">Newton
  97               Raphson</a>) the <code class="computeroutput"><span class="identifier">tuple</span></code>
  98               should have <span class="bold"><strong>two</strong></span> elements containing
  99               the evaluation of the function and its first derivative.
 100             </p>
 101 <p>
 102               For the third-order methods (<a href="http://en.wikipedia.org/wiki/Halley%27s_method" target="_top">Halley</a>
 103               and Schr&#246;der) the <code class="computeroutput"><span class="identifier">tuple</span></code>
 104               should have <span class="bold"><strong>three</strong></span> elements containing
 105               the evaluation of the function and its first and second derivatives.
 106             </p>
 107 </dd>
 108 <dt><span class="term">T guess</span></dt>
 109 <dd><p>
 110               The initial starting value. A good guess is crucial to quick convergence!
 111             </p></dd>
 112 <dt><span class="term">T min</span></dt>
 113 <dd><p>
 114               The minimum possible value for the result, this is used as an initial
 115               lower bracket.
 116             </p></dd>
 117 <dt><span class="term">T max</span></dt>
 118 <dd><p>
 119               The maximum possible value for the result, this is used as an initial
 120               upper bracket.
 121             </p></dd>
 122 <dt><span class="term">int digits</span></dt>
 123 <dd><p>
 124               The desired number of binary digits precision.
 125             </p></dd>
 126 <dt><span class="term">uintmax_t&amp; max_iter</span></dt>
 127 <dd><p>
 128               An optional maximum number of iterations to perform. On exit, this
 129               is updated to the actual number of iterations performed.
 130             </p></dd>
 131 </dl>
 132 </div>
 133 <p>
 134         When using these functions you should note that:
 135       </p>
 136 <div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; ">
 137 <li class="listitem">
 138             Default <code class="computeroutput"><span class="identifier">max_iter</span> <span class="special">=</span>
 139             <span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special">&lt;</span><span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span><span class="special">&gt;::</span><span class="identifier">max</span><span class="special">)()</span></code> is effectively 'iterate for ever'.
 140           </li>
 141 <li class="listitem">
 142             They may be very sensitive to the initial guess, typically they converge
 143             very rapidly if the initial guess has two or three decimal digits correct.
 144             However convergence can be no better than <a class="link" href="roots_noderiv/bisect.html" title="Bisection">bisect</a>,
 145             or in some rare cases, even worse than <a class="link" href="roots_noderiv/bisect.html" title="Bisection">bisect</a>
 146             if the initial guess is a long way from the correct value and the derivatives
 147             are close to zero.
 148           </li>
 149 <li class="listitem">
 150             These functions include special cases to handle zero first (and second
 151             where appropriate) derivatives, and fall back to <a class="link" href="roots_noderiv/bisect.html" title="Bisection">bisect</a>
 152             in this case. However, it is helpful if functor F is defined to return
 153             an arbitrarily small value <span class="emphasis"><em>of the correct sign</em></span> rather
 154             than zero.
 155           </li>
 156 <li class="listitem">
 157             If the derivative at the current best guess for the result is infinite
 158             (or very close to being infinite) then these functions may terminate
 159             prematurely. A large first derivative leads to a very small next step,
 160             triggering the termination condition. Derivative based iteration may
 161             not be appropriate in such cases.
 162           </li>
 163 <li class="listitem">
 164             If the function is 'Really Well Behaved' (is monotonic and has only one
 165             root) the bracket bounds <span class="emphasis"><em>min</em></span> and <span class="emphasis"><em>max</em></span>
 166             may as well be set to the widest limits like zero and <code class="computeroutput"><span class="identifier">numeric_limits</span><span class="special">&lt;</span><span class="identifier">T</span><span class="special">&gt;::</span><span class="identifier">max</span><span class="special">()</span></code>.
 167           </li>
 168 <li class="listitem">
 169             But if the function more complex and may have more than one root or a
 170             pole, the choice of bounds is protection against jumping out to seek
 171             the 'wrong' root.
 172           </li>
 173 <li class="listitem">
 174             These functions fall back to <a class="link" href="roots_noderiv/bisect.html" title="Bisection">bisect</a>
 175             if the next computed step would take the next value out of bounds. The
 176             bounds are updated after each step to ensure this leads to convergence.
 177             However, a good initial guess backed up by asymptotically-tight bounds
 178             will improve performance no end - rather than relying on <a class="link" href="roots_noderiv/bisect.html" title="Bisection">bisection</a>.
 179           </li>
 180 <li class="listitem">
 181             The value of <span class="emphasis"><em>digits</em></span> is crucial to good performance
 182             of these functions, if it is set too high then at best you will get one
 183             extra (unnecessary) iteration, and at worst the last few steps will proceed
 184             by <a class="link" href="roots_noderiv/bisect.html" title="Bisection">bisection</a>.
 185             Remember that the returned value can never be more accurate than <span class="emphasis"><em>f(x)</em></span>
 186             can be evaluated, and that if <span class="emphasis"><em>f(x)</em></span> suffers from
 187             cancellation errors as it tends to zero then the computed steps will
 188             be effectively random. The value of <span class="emphasis"><em>digits</em></span> should
 189             be set so that iteration terminates before this point: remember that
 190             for second and third order methods the number of correct digits in the
 191             result is increasing quite substantially with each iteration, <span class="emphasis"><em>digits</em></span>
 192             should be set by experiment so that the final iteration just takes the
 193             next value into the zone where <span class="emphasis"><em>f(x)</em></span> becomes inaccurate.
 194             A good starting point for <span class="emphasis"><em>digits</em></span> would be 0.6*D
 195             for Newton and 0.4*D for Halley or Shr&#246;der iteration, where D is <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special">&lt;</span><span class="identifier">T</span><span class="special">&gt;::</span><span class="identifier">digits</span></code>.
 196           </li>
 197 <li class="listitem">
 198             If you need some diagnostic output to see what is going on, you can
 199             <code class="computeroutput"><span class="preprocessor">#define</span> <span class="identifier">BOOST_MATH_INSTRUMENT</span></code>
 200             before the <code class="computeroutput"><span class="preprocessor">#include</span> <span class="special">&lt;</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">tools</span><span class="special">/</span><span class="identifier">roots</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">&gt;</span></code>, and also ensure that display of
 201             all the significant digits with <code class="computeroutput"> <span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">&gt;::</span><span class="identifier">digits10</span><span class="special">)</span></code>: or even possibly significant digits
 202             with <code class="computeroutput"> <span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">&gt;::</span><span class="identifier">max_digits10</span><span class="special">)</span></code>:
 203             but be warned, this may produce copious output!
 204           </li>
 205 <li class="listitem">
 206             Finally: you may well be able to do better than these functions by hand-coding
 207             the heuristics used so that they are tailored to a specific function.
 208             You may also be able to compute the ratio of derivatives used by these
 209             methods more efficiently than computing the derivatives themselves. As
 210             ever, algebraic simplification can be a big win.
 211           </li>
 212 </ul></div>
 213 <h5>
 214 <a name="math_toolkit.roots.roots_deriv.h2"></a>
 215         <span class="phrase"><a name="math_toolkit.roots.roots_deriv.newton"></a></span><a class="link" href="roots_deriv.html#math_toolkit.roots.roots_deriv.newton">Newton
 216         Raphson Method</a>
 217       </h5>
 218 <p>
 219         Given an initial guess <span class="emphasis"><em>x0</em></span> the subsequent values are
 220         computed using:
 221       </p>
 222 <p>
 223         <span class="inlinemediaobject"><img src="../../../equations/roots1.svg"></span>
 224       </p>
 225 <p>
 226         Out of bounds steps revert to <a class="link" href="roots_noderiv/bisect.html" title="Bisection">bisection</a>
 227         of the current bounds.
 228       </p>
 229 <p>
 230         Under ideal conditions, the number of correct digits doubles with each iteration.
 231       </p>
 232 <h5>
 233 <a name="math_toolkit.roots.roots_deriv.h3"></a>
 234         <span class="phrase"><a name="math_toolkit.roots.roots_deriv.halley"></a></span><a class="link" href="roots_deriv.html#math_toolkit.roots.roots_deriv.halley">Halley's
 235         Method</a>
 236       </h5>
 237 <p>
 238         Given an initial guess <span class="emphasis"><em>x0</em></span> the subsequent values are
 239         computed using:
 240       </p>
 241 <p>
 242         <span class="inlinemediaobject"><img src="../../../equations/roots2.svg"></span>
 243       </p>
 244 <p>
 245         Over-compensation by the second derivative (one which would proceed in the
 246         wrong direction) causes the method to revert to a Newton-Raphson step.
 247       </p>
 248 <p>
 249         Out of bounds steps revert to bisection of the current bounds.
 250       </p>
 251 <p>
 252         Under ideal conditions, the number of correct digits trebles with each iteration.
 253       </p>
 254 <h5>
 255 <a name="math_toolkit.roots.roots_deriv.h4"></a>
 256         <span class="phrase"><a name="math_toolkit.roots.roots_deriv.schroder"></a></span><a class="link" href="roots_deriv.html#math_toolkit.roots.roots_deriv.schroder">Schr&#246;der's
 257         Method</a>
 258       </h5>
 259 <p>
 260         Given an initial guess x0 the subsequent values are computed using:
 261       </p>
 262 <p>
 263         <span class="inlinemediaobject"><img src="../../../equations/roots3.svg"></span>
 264       </p>
 265 <p>
 266         Over-compensation by the second derivative (one which would proceed in the
 267         wrong direction) causes the method to revert to a Newton-Raphson step. Likewise
 268         a Newton step is used whenever that Newton step would change the next value
 269         by more than 10%.
 270       </p>
 271 <p>
 272         Out of bounds steps revert to <a href="https://en.wikipedia.org/wiki/Bisection" target="_top">bisection</a>
 273         of the current bounds.
 274       </p>
 275 <p>
 276         Under ideal conditions, the number of correct digits trebles with each iteration.
 277       </p>
 278 <p>
 279         This is Schr&#246;der's general result (equation 18 from <a href="http://drum.lib.umd.edu/handle/1903/577" target="_top">Stewart,
 280         G. W. "On Infinitely Many Algorithms for Solving Equations." English
 281         translation of Schr&#246;der's original paper. College Park, MD: University of
 282         Maryland, Institute for Advanced Computer Studies, Department of Computer
 283         Science, 1993</a>.)
 284       </p>
 285 <p>
 286         This method guarantees at least quadratic convergence (the same as Newton's
 287         method), and is known to work well in the presence of multiple roots: something
 288         that neither Newton nor Halley can do.
 289       </p>
 290 <h5>
 291 <a name="math_toolkit.roots.roots_deriv.h5"></a>
 292         <span class="phrase"><a name="math_toolkit.roots.roots_deriv.examples"></a></span><a class="link" href="roots_deriv.html#math_toolkit.roots.roots_deriv.examples">Examples</a>
 293       </h5>
 294 <p>
 295         See <a class="link" href="root_finding_examples.html" title="Examples of Root-Finding (with and without derivatives)">root-finding
 296         examples</a>.
 297       </p>
 298 </div>
 299 <table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
 300 <td align="left"></td>
 301 <td align="right"><div class="copyright-footer">Copyright &#169; 2006-2010, 2012-2014 Nikhar Agrawal,
 302       Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, Hubert
 303       Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Johan R&#229;de, Gautam Sewani,
 304       Benjamin Sobotta, Thijs van den Berg, Daryle Walker and Xiaogang Zhang<p>
 305         Distributed under the Boost Software License, Version 1.0. (See accompanying
 306         file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>)
 307       </p>
 308 </div></td>
 309 </tr></table>
 310 <hr>
 311 <div class="spirit-nav">
 312 <a accesskey="p" href="roots_noderiv/implementation.html"><img src="../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../roots.html"><img src="../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="root_finding_examples.html"><img src="../../../../../../doc/src/images/next.png" alt="Next"></a>
 313 </div>
 314 </body>
 315 </html>