roots_deriv.html 33 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379
  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.79.1">
  7. <link rel="home" href="../index.html" title="Math Toolkit 2.11.0">
  8. <link rel="up" href="../root_finding.html" title="Chapter&#160;10.&#160;Root Finding &amp; Minimization Algorithms">
  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="../root_finding.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><h2 class="title" style="clear: both">
  27. <a name="math_toolkit.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. </h2></div></div></div>
  30. <h5>
  31. <a name="math_toolkit.roots_deriv.h0"></a>
  32. <span class="phrase"><a name="math_toolkit.roots_deriv.synopsis"></a></span><a class="link" href="roots_deriv.html#math_toolkit.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. <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>
  42. <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>
  43. <span class="comment">// Halley</span>
  44. <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>
  45. <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>
  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> <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>
  48. <span class="comment">// Schr'''&amp;#xf6;'''der</span>
  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">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>
  51. <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>
  52. <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>
  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">Complex</span><span class="special">&gt;</span>
  54. <span class="identifier">Complex</span> <span class="identifier">complex_newton</span><span class="special">(</span><span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">Complex</span> <span class="identifier">guess</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">max_iterations</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">typename</span> <span class="identifier">Complex</span><span class="special">::</span><span class="identifier">value_type</span><span class="special">&gt;::</span><span class="identifier">digits</span><span class="special">);</span>
  55. <span class="keyword">template</span><span class="special">&lt;</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">&gt;</span>
  56. <span class="keyword">auto</span> <span class="identifier">quadratic_roots</span><span class="special">(</span><span class="identifier">T</span> <span class="keyword">const</span> <span class="special">&amp;</span> <span class="identifier">a</span><span class="special">,</span> <span class="identifier">T</span> <span class="keyword">const</span> <span class="special">&amp;</span> <span class="identifier">b</span><span class="special">,</span> <span class="identifier">T</span> <span class="keyword">const</span> <span class="special">&amp;</span> <span class="identifier">c</span><span class="special">);</span>
  57. <span class="special">}}}</span> <span class="comment">// namespaces boost::math::tools.</span>
  58. </pre>
  59. <h5>
  60. <a name="math_toolkit.roots_deriv.h1"></a>
  61. <span class="phrase"><a name="math_toolkit.roots_deriv.description"></a></span><a class="link" href="roots_deriv.html#math_toolkit.roots_deriv.description">Description</a>
  62. </h5>
  63. <p>
  64. These functions all perform iterative root-finding <span class="bold"><strong>using
  65. derivatives</strong></span>:
  66. </p>
  67. <div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; ">
  68. <li class="listitem">
  69. <code class="computeroutput"><span class="identifier">newton_raphson_iterate</span></code>
  70. performs second-order <a class="link" href="roots_deriv.html#math_toolkit.roots_deriv.newton">Newton-Raphson
  71. iteration</a>.
  72. </li>
  73. <li class="listitem">
  74. <code class="computeroutput"><span class="identifier">halley_iterate</span></code> and <code class="computeroutput"><span class="identifier">schroder_iterate</span></code> perform third-order
  75. <a class="link" href="roots_deriv.html#math_toolkit.roots_deriv.halley">Halley</a> and <a class="link" href="roots_deriv.html#math_toolkit.roots_deriv.schroder">Schr&#246;der</a> iteration.
  76. </li>
  77. <li class="listitem">
  78. <code class="computeroutput"><span class="identifier">complex_newton</span></code> performs
  79. Newton's method on complex-analytic functions.
  80. </li>
  81. <li class="listitem">
  82. <code class="computeroutput"><span class="identifier">solve_quadratic</span></code> solves
  83. quadratic equations using various tricks to keep catastrophic cancellation
  84. from occurring in computation of the discriminant.
  85. </li>
  86. </ul></div>
  87. <div class="variablelist">
  88. <p class="title"><b>Parameters of the real-valued root finding functions</b></p>
  89. <dl class="variablelist">
  90. <dt><span class="term">F f</span></dt>
  91. <dd>
  92. <p>
  93. Type F must be a callable function object (or C++ lambda) that accepts
  94. one parameter and returns a <a class="link" href="internals/tuples.html" title="Tuples">std::pair,
  95. std::tuple, boost::tuple or boost::fusion::tuple</a>:
  96. </p>
  97. <p>
  98. For second-order iterative method (<a href="http://en.wikipedia.org/wiki/Newton_Raphson" target="_top">Newton
  99. Raphson</a>) the <code class="computeroutput"><span class="identifier">tuple</span></code>
  100. should have <span class="bold"><strong>two</strong></span> elements containing
  101. the evaluation of the function and its first derivative.
  102. </p>
  103. <p>
  104. For the third-order methods (<a href="http://en.wikipedia.org/wiki/Halley%27s_method" target="_top">Halley</a>
  105. and Schr&#246;der) the <code class="computeroutput"><span class="identifier">tuple</span></code>
  106. should have <span class="bold"><strong>three</strong></span> elements containing
  107. the evaluation of the function and its first and second derivatives.
  108. </p>
  109. </dd>
  110. <dt><span class="term">T guess</span></dt>
  111. <dd><p>
  112. The initial starting value. A good guess is crucial to quick convergence!
  113. </p></dd>
  114. <dt><span class="term">T min</span></dt>
  115. <dd><p>
  116. The minimum possible value for the result, this is used as an initial
  117. lower bracket.
  118. </p></dd>
  119. <dt><span class="term">T max</span></dt>
  120. <dd><p>
  121. The maximum possible value for the result, this is used as an initial
  122. upper bracket.
  123. </p></dd>
  124. <dt><span class="term">int digits</span></dt>
  125. <dd><p>
  126. The desired number of binary digits precision.
  127. </p></dd>
  128. <dt><span class="term">uintmax_t&amp; max_iter</span></dt>
  129. <dd><p>
  130. An optional maximum number of iterations to perform. On exit, this is
  131. updated to the actual number of iterations performed.
  132. </p></dd>
  133. </dl>
  134. </div>
  135. <p>
  136. When using these functions you should note that:
  137. </p>
  138. <div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; ">
  139. <li class="listitem">
  140. Default <code class="computeroutput"><span class="identifier">max_iter</span> <span class="special">=</span>
  141. <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'.
  142. </li>
  143. <li class="listitem">
  144. They may be very sensitive to the initial guess, typically they converge
  145. very rapidly if the initial guess has two or three decimal digits correct.
  146. However convergence can be no better than <a class="link" href="roots_noderiv/bisect.html" title="Bisection">bisect</a>,
  147. or in some rare cases, even worse than <a class="link" href="roots_noderiv/bisect.html" title="Bisection">bisect</a>
  148. if the initial guess is a long way from the correct value and the derivatives
  149. are close to zero.
  150. </li>
  151. <li class="listitem">
  152. These functions include special cases to handle zero first (and second
  153. where appropriate) derivatives, and fall back to <a class="link" href="roots_noderiv/bisect.html" title="Bisection">bisect</a>
  154. in this case. However, it is helpful if functor F is defined to return
  155. an arbitrarily small value <span class="emphasis"><em>of the correct sign</em></span> rather
  156. than zero.
  157. </li>
  158. <li class="listitem">
  159. The functions will raise an <a class="link" href="error_handling.html#math_toolkit.error_handling.evaluation_error">evaluation_error</a>
  160. if arguments <code class="computeroutput"><span class="identifier">min</span></code> and <code class="computeroutput"><span class="identifier">max</span></code> are the wrong way around or if they
  161. converge to a local minima.
  162. </li>
  163. <li class="listitem">
  164. If the derivative at the current best guess for the result is infinite
  165. (or very close to being infinite) then these functions may terminate prematurely.
  166. A large first derivative leads to a very small next step, triggering the
  167. termination condition. Derivative based iteration may not be appropriate
  168. in such cases.
  169. </li>
  170. <li class="listitem">
  171. If the function is 'Really Well Behaved' (is monotonic and has only one
  172. root) the bracket bounds <span class="emphasis"><em>min</em></span> and <span class="emphasis"><em>max</em></span>
  173. 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>.
  174. </li>
  175. <li class="listitem">
  176. But if the function more complex and may have more than one root or a pole,
  177. the choice of bounds is protection against jumping out to seek the 'wrong'
  178. root.
  179. </li>
  180. <li class="listitem">
  181. These functions fall back to <a class="link" href="roots_noderiv/bisect.html" title="Bisection">bisect</a>
  182. if the next computed step would take the next value out of bounds. The
  183. bounds are updated after each step to ensure this leads to convergence.
  184. However, a good initial guess backed up by asymptotically-tight bounds
  185. will improve performance no end - rather than relying on <a class="link" href="roots_noderiv/bisect.html" title="Bisection">bisection</a>.
  186. </li>
  187. <li class="listitem">
  188. The value of <span class="emphasis"><em>digits</em></span> is crucial to good performance
  189. of these functions, if it is set too high then at best you will get one
  190. extra (unnecessary) iteration, and at worst the last few steps will proceed
  191. by <a class="link" href="roots_noderiv/bisect.html" title="Bisection">bisection</a>.
  192. Remember that the returned value can never be more accurate than <span class="emphasis"><em>f(x)</em></span>
  193. can be evaluated, and that if <span class="emphasis"><em>f(x)</em></span> suffers from cancellation
  194. errors as it tends to zero then the computed steps will be effectively
  195. random. The value of <span class="emphasis"><em>digits</em></span> should be set so that
  196. iteration terminates before this point: remember that for second and third
  197. order methods the number of correct digits in the result is increasing
  198. quite substantially with each iteration, <span class="emphasis"><em>digits</em></span> should
  199. be set by experiment so that the final iteration just takes the next value
  200. into the zone where <span class="emphasis"><em>f(x)</em></span> becomes inaccurate. A good
  201. starting point for <span class="emphasis"><em>digits</em></span> would be 0.6*D for Newton
  202. 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>.
  203. </li>
  204. <li class="listitem">
  205. If you need some diagnostic output to see what is going on, you can <code class="computeroutput"><span class="preprocessor">#define</span> <span class="identifier">BOOST_MATH_INSTRUMENT</span></code>
  206. 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 all
  207. 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 with
  208. <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>:
  209. but be warned, this may produce copious output!
  210. </li>
  211. <li class="listitem">
  212. Finally: you may well be able to do better than these functions by hand-coding
  213. the heuristics used so that they are tailored to a specific function. You
  214. may also be able to compute the ratio of derivatives used by these methods
  215. more efficiently than computing the derivatives themselves. As ever, algebraic
  216. simplification can be a big win.
  217. </li>
  218. </ul></div>
  219. <h5>
  220. <a name="math_toolkit.roots_deriv.h2"></a>
  221. <span class="phrase"><a name="math_toolkit.roots_deriv.newton"></a></span><a class="link" href="roots_deriv.html#math_toolkit.roots_deriv.newton">Newton
  222. Raphson Method</a>
  223. </h5>
  224. <p>
  225. Given an initial guess <span class="emphasis"><em>x0</em></span> the subsequent values are computed
  226. using:
  227. </p>
  228. <div class="blockquote"><blockquote class="blockquote"><p>
  229. <span class="inlinemediaobject"><img src="../../equations/roots1.svg"></span>
  230. </p></blockquote></div>
  231. <p>
  232. Out-of-bounds steps revert to <a class="link" href="roots_noderiv/bisect.html" title="Bisection">bisection</a>
  233. of the current bounds.
  234. </p>
  235. <p>
  236. Under ideal conditions, the number of correct digits doubles with each iteration.
  237. </p>
  238. <h5>
  239. <a name="math_toolkit.roots_deriv.h3"></a>
  240. <span class="phrase"><a name="math_toolkit.roots_deriv.halley"></a></span><a class="link" href="roots_deriv.html#math_toolkit.roots_deriv.halley">Halley's
  241. Method</a>
  242. </h5>
  243. <p>
  244. Given an initial guess <span class="emphasis"><em>x0</em></span> the subsequent values are computed
  245. using:
  246. </p>
  247. <div class="blockquote"><blockquote class="blockquote"><p>
  248. <span class="inlinemediaobject"><img src="../../equations/roots2.svg"></span>
  249. </p></blockquote></div>
  250. <p>
  251. Over-compensation by the second derivative (one which would proceed in the
  252. wrong direction) causes the method to revert to a Newton-Raphson step.
  253. </p>
  254. <p>
  255. Out of bounds steps revert to bisection of the current bounds.
  256. </p>
  257. <p>
  258. Under ideal conditions, the number of correct digits trebles with each iteration.
  259. </p>
  260. <h5>
  261. <a name="math_toolkit.roots_deriv.h4"></a>
  262. <span class="phrase"><a name="math_toolkit.roots_deriv.schroder"></a></span><a class="link" href="roots_deriv.html#math_toolkit.roots_deriv.schroder">Schr&#246;der's
  263. Method</a>
  264. </h5>
  265. <p>
  266. Given an initial guess x0 the subsequent values are computed using:
  267. </p>
  268. <div class="blockquote"><blockquote class="blockquote"><p>
  269. <span class="inlinemediaobject"><img src="../../equations/roots3.svg"></span>
  270. </p></blockquote></div>
  271. <p>
  272. Over-compensation by the second derivative (one which would proceed in the
  273. wrong direction) causes the method to revert to a Newton-Raphson step. Likewise
  274. a Newton step is used whenever that Newton step would change the next value
  275. by more than 10%.
  276. </p>
  277. <p>
  278. Out of bounds steps revert to <a href="https://en.wikipedia.org/wiki/Bisection" target="_top">bisection</a>
  279. of the current bounds.
  280. </p>
  281. <p>
  282. Under ideal conditions, the number of correct digits trebles with each iteration.
  283. </p>
  284. <p>
  285. This is Schr&#246;der's general result (equation 18 from <a href="http://drum.lib.umd.edu/handle/1903/577" target="_top">Stewart,
  286. G. W. "On Infinitely Many Algorithms for Solving Equations." English
  287. translation of Schr&#246;der's original paper. College Park, MD: University of Maryland,
  288. Institute for Advanced Computer Studies, Department of Computer Science, 1993</a>.)
  289. </p>
  290. <p>
  291. This method guarantees at least quadratic convergence (the same as Newton's
  292. method), and is known to work well in the presence of multiple roots: something
  293. that neither Newton nor Halley can do.
  294. </p>
  295. <p>
  296. The complex Newton method works slightly differently than the rest of the methods:
  297. Since there is no way to bracket roots in the complex plane, the <code class="computeroutput"><span class="identifier">min</span></code> and <code class="computeroutput"><span class="identifier">max</span></code>
  298. arguments are not accepted. Failure to reach a root is communicated by returning
  299. <code class="computeroutput"><span class="identifier">nan</span></code>s. Remember that if a function
  300. has many roots, then which root the complex Newton's method converges to is
  301. essentially impossible to predict a priori; see the Newton's fractal for more
  302. information.
  303. </p>
  304. <p>
  305. Finally, the derivative of <span class="emphasis"><em>f</em></span> must be continuous at the
  306. root or else non-roots can be found; see <a href="https://math.stackexchange.com/questions/3017766/constructing-newton-iteration-converging-to-non-root" target="_top">here</a>
  307. for an example.
  308. </p>
  309. <p>
  310. An example usage of <code class="computeroutput"><span class="identifier">complex_newton</span></code>
  311. is given in <code class="computeroutput"><span class="identifier">examples</span><span class="special">/</span><span class="identifier">daubechies_coefficients</span><span class="special">.</span><span class="identifier">cpp</span></code>.
  312. </p>
  313. <h5>
  314. <a name="math_toolkit.roots_deriv.h5"></a>
  315. <span class="phrase"><a name="math_toolkit.roots_deriv.quadratics"></a></span><a class="link" href="roots_deriv.html#math_toolkit.roots_deriv.quadratics">Quadratics</a>
  316. </h5>
  317. <p>
  318. To solve a quadratic <span class="emphasis"><em>ax</em></span><sup>2</sup> + <span class="emphasis"><em>bx</em></span> + <span class="emphasis"><em>c</em></span>
  319. = 0, we may use
  320. </p>
  321. <pre class="programlisting"><span class="keyword">auto</span> <span class="special">[</span><span class="identifier">x0</span><span class="special">,</span> <span class="identifier">x1</span><span class="special">]</span> <span class="special">=</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">quadratic_roots</span><span class="special">(</span><span class="identifier">a</span><span class="special">,</span> <span class="identifier">b</span><span class="special">,</span> <span class="identifier">c</span><span class="special">);</span>
  322. </pre>
  323. <p>
  324. If the roots are real, they are arranged so that <code class="computeroutput"><span class="identifier">x0</span></code>
  325. &#8804; <code class="computeroutput"><span class="identifier">x1</span></code>. If the roots are
  326. complex and the inputs are real, <code class="computeroutput"><span class="identifier">x0</span></code>
  327. and <code class="computeroutput"><span class="identifier">x1</span></code> are both <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">Real</span><span class="special">&gt;::</span><span class="identifier">quiet_NaN</span><span class="special">()</span></code>. In this case we must cast <code class="computeroutput"><span class="identifier">a</span></code>, <code class="computeroutput"><span class="identifier">b</span></code>
  328. and <code class="computeroutput"><span class="identifier">c</span></code> to a complex type to
  329. extract the complex roots. If <code class="computeroutput"><span class="identifier">a</span></code>,
  330. <code class="computeroutput"><span class="identifier">b</span></code> and <code class="computeroutput"><span class="identifier">c</span></code>
  331. are integral, then the roots are of type double. The routine is much faster
  332. if the fused-multiply-add instruction is available on your architecture. If
  333. the fma is not available, the function resorts to slow emulation. Finally,
  334. speed is improved if you compile for your particular architecture. For instance,
  335. if you compile without any architecture flags, then the <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">fma</span></code> call
  336. compiles down to <code class="computeroutput"><span class="identifier">call</span> <span class="identifier">_fma</span></code>,
  337. which dynamically chooses to emulate or execute the <code class="computeroutput"><span class="identifier">vfmadd132sd</span></code>
  338. instruction based on the capabilities of the architecture. If instead, you
  339. compile with (say) <code class="computeroutput"><span class="special">-</span><span class="identifier">march</span><span class="special">=</span><span class="identifier">native</span></code> then
  340. no dynamic choice is made: The <code class="computeroutput"><span class="identifier">vfmadd132sd</span></code>
  341. instruction is always executed if available and emulation is used if not.
  342. </p>
  343. <h5>
  344. <a name="math_toolkit.roots_deriv.h6"></a>
  345. <span class="phrase"><a name="math_toolkit.roots_deriv.examples"></a></span><a class="link" href="roots_deriv.html#math_toolkit.roots_deriv.examples">Examples</a>
  346. </h5>
  347. <p>
  348. See <a class="link" href="root_finding_examples.html" title="Examples of Root-Finding (with and without derivatives)">root-finding examples</a>.
  349. </p>
  350. </div>
  351. <table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
  352. <td align="left"></td>
  353. <td align="right"><div class="copyright-footer">Copyright &#169; 2006-2019 Nikhar
  354. Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos,
  355. Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Matthew Pulver, Johan
  356. R&#229;de, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg,
  357. Daryle Walker and Xiaogang Zhang<p>
  358. Distributed under the Boost Software License, Version 1.0. (See accompanying
  359. 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>)
  360. </p>
  361. </div></td>
  362. </tr></table>
  363. <hr>
  364. <div class="spirit-nav">
  365. <a accesskey="p" href="roots_noderiv/implementation.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../root_finding.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>
  366. </div>
  367. </body>
  368. </html>