]>
Commit | Line | Data |
---|---|---|
1 | [section:roots_deriv Root Finding With Derivatives: Newton-Raphson, Halley & Schr'''ö'''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'''ö'''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 | }}} // namespaces boost::math::tools. | |
33 | ||
34 | [h4 Description] | |
35 | ||
36 | These functions all perform iterative root-finding [*using derivatives]: | |
37 | ||
38 | * `newton_raphson_iterate` performs second-order __newton. | |
39 | ||
40 | * `halley_iterate` and `schroder_iterate` perform third-order | |
41 | __halley and __schroder iteration. | |
42 | ||
43 | The functions all take the same parameters: | |
44 | ||
45 | [variablelist Parameters of the root finding functions | |
46 | [[F f] [Type F must be a callable function object that accepts one parameter and | |
47 | returns a __tuple_type: | |
48 | ||
49 | For second-order iterative method ([@http://en.wikipedia.org/wiki/Newton_Raphson Newton Raphson]) | |
50 | the `tuple` should have [*two] elements containing the evaluation | |
51 | of the function and its first derivative. | |
52 | ||
53 | For the third-order methods | |
54 | ([@http://en.wikipedia.org/wiki/Halley%27s_method Halley] and | |
55 | Schr'''ö'''der) | |
56 | the `tuple` should have [*three] elements containing the evaluation of | |
57 | the function and its first and second derivatives.]] | |
58 | [[T guess] [The initial starting value. A good guess is crucial to quick convergence!]] | |
59 | [[T min] [The minimum possible value for the result, this is used as an initial lower bracket.]] | |
60 | [[T max] [The maximum possible value for the result, this is used as an initial upper bracket.]] | |
61 | [[int digits] [The desired number of binary digits precision.]] | |
62 | [[uintmax_t& max_iter] [An optional maximum number of iterations to perform. On exit, this is updated to the actual number of iterations performed.]] | |
63 | ] | |
64 | ||
65 | When using these functions you should note that: | |
66 | ||
67 | * Default `max_iter = (std::numeric_limits<boost::uintmax_t>::max)()` is effectively 'iterate for ever'. | |
68 | * They may be very sensitive to the initial guess, typically they converge very rapidly | |
69 | if the initial guess has two or three decimal digits correct. However convergence | |
70 | can be no better than __bisect, or in some rare cases, even worse than __bisect if the | |
71 | initial guess is a long way from the correct value and the derivatives are close to zero. | |
72 | * These functions include special cases to handle zero first (and second where appropriate) | |
73 | derivatives, and fall back to __bisect in this case. However, it is helpful | |
74 | if functor F is defined to return an arbitrarily small value ['of the correct sign] rather | |
75 | than zero. | |
76 | * If the derivative at the current best guess for the result is infinite (or | |
77 | very close to being infinite) then these functions may terminate prematurely. | |
78 | A large first derivative leads to a very small next step, triggering the termination | |
79 | condition. Derivative based iteration may not be appropriate in such cases. | |
80 | * If the function is 'Really Well Behaved' (is monotonic and has only one root) | |
81 | the bracket bounds ['min] and ['max] may as well be set to the widest limits | |
82 | like zero and `numeric_limits<T>::max()`. | |
83 | *But if the function more complex and may have more than one root or a pole, | |
84 | the choice of bounds is protection against jumping out to seek the 'wrong' root. | |
85 | * These functions fall back to __bisect if the next computed step would take the | |
86 | next value out of bounds. The bounds are updated after each step to ensure this leads | |
87 | to convergence. However, a good initial guess backed up by asymptotically-tight | |
88 | bounds will improve performance no end - rather than relying on __bisection. | |
89 | * The value of ['digits] is crucial to good performance of these functions, | |
90 | if it is set too high then at best you will get one extra (unnecessary) | |
91 | iteration, and at worst the last few steps will proceed by __bisection. | |
92 | Remember that the returned value can never be more accurate than ['f(x)] can be | |
93 | evaluated, and that if ['f(x)] suffers from cancellation errors as it | |
94 | tends to zero then the computed steps will be effectively random. The | |
95 | value of ['digits] should be set so that iteration terminates before this point: | |
96 | remember that for second and third order methods the number of correct | |
97 | digits in the result is increasing quite | |
98 | substantially with each iteration, ['digits] should be set by experiment so that the final | |
99 | iteration just takes the next value into the zone where ['f(x)] becomes inaccurate. | |
100 | A good starting point for ['digits] would be 0.6*D for Newton and 0.4*D for Halley or Shr'''ö'''der | |
101 | iteration, where D is `std::numeric_limits<T>::digits`. | |
102 | * If you need some diagnostic output to see what is going on, you can | |
103 | `#define BOOST_MATH_INSTRUMENT` before the `#include <boost/math/tools/roots.hpp>`, | |
104 | and also ensure that display of all the significant digits with | |
105 | ` cout.precision(std::numeric_limits<double>::digits10)`: | |
106 | or even possibly significant digits with | |
107 | ` cout.precision(std::numeric_limits<double>::max_digits10)`: | |
108 | but be warned, this may produce copious output! | |
109 | * Finally: you may well be able to do better than these functions by hand-coding | |
110 | the heuristics used so that they are tailored to a specific function. You may also | |
111 | be able to compute the ratio of derivatives used by these methods more efficiently | |
112 | than computing the derivatives themselves. As ever, algebraic simplification can | |
113 | be a big win. | |
114 | ||
115 | [h4:newton Newton Raphson Method] | |
116 | ||
117 | Given an initial guess ['x0] the subsequent values are computed using: | |
118 | ||
119 | [equation roots1] | |
120 | ||
121 | Out of bounds steps revert to __bisection of the current bounds. | |
122 | ||
123 | Under ideal conditions, the number of correct digits doubles with each iteration. | |
124 | ||
125 | [h4:halley Halley's Method] | |
126 | ||
127 | Given an initial guess ['x0] the subsequent values are computed using: | |
128 | ||
129 | [equation roots2] | |
130 | ||
131 | Over-compensation by the second derivative (one which would proceed | |
132 | in the wrong direction) causes the method to | |
133 | revert to a Newton-Raphson step. | |
134 | ||
135 | Out of bounds steps revert to bisection of the current bounds. | |
136 | ||
137 | Under ideal conditions, the number of correct digits trebles with each iteration. | |
138 | ||
139 | [h4:schroder Schr'''ö'''der's Method] | |
140 | ||
141 | Given an initial guess x0 the subsequent values are computed using: | |
142 | ||
143 | [equation roots3] | |
144 | ||
145 | Over-compensation by the second derivative (one which would proceed | |
146 | in the wrong direction) causes the method to | |
147 | revert to a Newton-Raphson step. Likewise a Newton step is used | |
148 | whenever that Newton step would change the next value by more than 10%. | |
149 | ||
150 | Out of bounds steps revert to __bisection_wikipedia of the current bounds. | |
151 | ||
152 | Under ideal conditions, the number of correct digits trebles with each iteration. | |
153 | ||
154 | This is Schr'''ö'''der's general result (equation 18 from [@http://drum.lib.umd.edu/handle/1903/577 Stewart, G. W. | |
155 | "On Infinitely Many Algorithms for Solving Equations." English translation of Schr'''ö'''der's original paper. | |
156 | College Park, MD: University of Maryland, Institute for Advanced Computer Studies, Department of Computer Science, 1993].) | |
157 | ||
158 | 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: | |
159 | something that neither Newton nor Halley can do. | |
160 | ||
161 | [h4 Examples] | |
162 | ||
163 | See __root_finding_examples. | |
164 | ||
165 | [endsect] [/section:roots_deriv Root Finding With Derivatives] | |
166 | ||
167 | [/ | |
168 | Copyright 2006, 2010, 2012 John Maddock and Paul A. Bristow. | |
169 | Distributed under the Boost Software License, Version 1.0. | |
170 | (See accompanying file LICENSE_1_0.txt or copy at | |
171 | http://www.boost.org/LICENSE_1_0.txt). | |
172 | ] |