]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // root_finding_example.cpp |
2 | ||
3 | // Copyright Paul A. Bristow 2010, 2015 | |
4 | ||
5 | // Use, modification and distribution are subject to the | |
6 | // Boost Software License, Version 1.0. | |
7 | // (See accompanying file LICENSE_1_0.txt | |
8 | // or copy at http://www.boost.org/LICENSE_1_0.txt) | |
9 | ||
10 | // Example of finding roots using Newton-Raphson, Halley. | |
11 | ||
12 | // Note that this file contains Quickbook mark-up as well as code | |
13 | // and comments, don't change any of the special comment mark-ups! | |
14 | ||
15 | //#define BOOST_MATH_INSTRUMENT | |
16 | ||
17 | /* | |
18 | This example demonstrates how to use the various tools for root finding | |
19 | taking the simple cube root function (`cbrt`) as an example. | |
20 | ||
21 | It shows how use of derivatives can improve the speed. | |
22 | (But is only a demonstration and does not try to make the ultimate improvements of 'real-life' | |
23 | implementation of `boost::math::cbrt`, mainly by using a better computed initial 'guess' | |
24 | at `<boost/math/special_functions/cbrt.hpp>`). | |
25 | ||
26 | Then we show how a higher root (fifth) can be computed, | |
27 | and in `root_finding_n_example.cpp` a generic method | |
28 | for the ['n]th root that constructs the derivatives at compile-time, | |
29 | ||
30 | These methods should be applicable to other functions that can be differentiated easily. | |
31 | ||
32 | First some `#includes` that will be needed. | |
33 | ||
34 | [tip For clarity, `using` statements are provided to list what functions are being used in this example: | |
35 | you can of course partly or fully qualify the names in other ways. | |
36 | (For your application, you may wish to extract some parts into header files, | |
37 | but you should never use `using` statements globally in header files).] | |
38 | */ | |
39 | ||
40 | //[root_finding_include_1 | |
41 | ||
42 | #include <boost/math/tools/roots.hpp> | |
43 | //using boost::math::policies::policy; | |
44 | //using boost::math::tools::newton_raphson_iterate; | |
45 | //using boost::math::tools::halley_iterate; // | |
46 | //using boost::math::tools::eps_tolerance; // Binary functor for specified number of bits. | |
47 | //using boost::math::tools::bracket_and_solve_root; | |
48 | //using boost::math::tools::toms748_solve; | |
49 | ||
50 | #include <boost/math/special_functions/next.hpp> // For float_distance. | |
51 | #include <tuple> // for std::tuple and std::make_tuple. | |
52 | #include <boost/math/special_functions/cbrt.hpp> // For boost::math::cbrt. | |
53 | ||
54 | //] [/root_finding_include_1] | |
55 | ||
56 | // using boost::math::tuple; | |
57 | // using boost::math::make_tuple; | |
58 | // using boost::math::tie; | |
59 | // which provide convenient aliases for various implementations, | |
60 | // including std::tr1, depending on what is available. | |
61 | ||
62 | #include <iostream> | |
63 | //using std::cout; using std::endl; | |
64 | #include <iomanip> | |
65 | //using std::setw; using std::setprecision; | |
66 | #include <limits> | |
67 | //using std::numeric_limits; | |
68 | ||
69 | /* | |
70 | ||
71 | Let's suppose we want to find the root of a number ['a], and to start, compute the cube root. | |
72 | ||
73 | So the equation we want to solve is: | |
74 | ||
75 | __spaces ['f](x) = x[cubed] - a | |
76 | ||
77 | We will first solve this without using any information | |
78 | about the slope or curvature of the cube root function. | |
79 | ||
80 | We then show how adding what we can know about this function, first just the slope, | |
81 | the 1st derivation /f'(x)/, will speed homing in on the solution. | |
82 | ||
83 | Lastly we show how adding the curvature /f''(x)/ too will speed convergence even more. | |
84 | ||
85 | */ | |
86 | ||
87 | //[root_finding_noderiv_1 | |
88 | ||
89 | template <class T> | |
90 | struct cbrt_functor_noderiv | |
91 | { | |
92 | // cube root of x using only function - no derivatives. | |
93 | cbrt_functor_noderiv(T const& to_find_root_of) : a(to_find_root_of) | |
94 | { /* Constructor just stores value a to find root of. */ } | |
95 | T operator()(T const& x) | |
96 | { | |
97 | T fx = x*x*x - a; // Difference (estimate x^3 - a). | |
98 | return fx; | |
99 | } | |
100 | private: | |
101 | T a; // to be 'cube_rooted'. | |
102 | }; | |
103 | //] [/root_finding_noderiv_1 | |
104 | ||
105 | /* | |
106 | Implementing the cube root function itself is fairly trivial now: | |
107 | the hardest part is finding a good approximation to begin with. | |
108 | In this case we'll just divide the exponent by three. | |
109 | (There are better but more complex guess algorithms used in 'real-life'.) | |
110 | ||
111 | Cube root function is 'Really Well Behaved' in that it is monotonic | |
112 | and has only one root (we leave negative values 'as an exercise for the student'). | |
113 | */ | |
114 | ||
115 | //[root_finding_noderiv_2 | |
116 | ||
117 | template <class T> | |
118 | T cbrt_noderiv(T x) | |
119 | { | |
120 | // return cube root of x using bracket_and_solve (no derivatives). | |
121 | using namespace std; // Help ADL of std functions. | |
122 | using namespace boost::math::tools; // For bracket_and_solve_root. | |
123 | ||
124 | int exponent; | |
125 | frexp(x, &exponent); // Get exponent of z (ignore mantissa). | |
126 | T guess = ldexp(1., exponent/3); // Rough guess is to divide the exponent by three. | |
127 | T factor = 2; // How big steps to take when searching. | |
128 | ||
129 | const boost::uintmax_t maxit = 20; // Limit to maximum iterations. | |
f67539c2 | 130 | boost::uintmax_t it = maxit; // Initially our chosen max iterations, but updated with actual. |
7c673cae FG |
131 | bool is_rising = true; // So if result if guess^3 is too low, then try increasing guess. |
132 | int digits = std::numeric_limits<T>::digits; // Maximum possible binary digits accuracy for type T. | |
133 | // Some fraction of digits is used to control how accurate to try to make the result. | |
134 | int get_digits = digits - 3; // We have to have a non-zero interval at each step, so | |
135 | // maximum accuracy is digits - 1. But we also have to | |
136 | // allow for inaccuracy in f(x), otherwise the last few | |
137 | // iterations just thrash around. | |
138 | eps_tolerance<T> tol(get_digits); // Set the tolerance. | |
139 | std::pair<T, T> r = bracket_and_solve_root(cbrt_functor_noderiv<T>(x), guess, factor, is_rising, tol, it); | |
140 | return r.first + (r.second - r.first)/2; // Midway between brackets is our result, if necessary we could | |
141 | // return the result as an interval here. | |
142 | } | |
143 | ||
144 | /*` | |
145 | ||
146 | [note The final parameter specifying a maximum number of iterations is optional. | |
147 | However, it defaults to `boost::uintmax_t maxit = (std::numeric_limits<boost::uintmax_t>::max)();` | |
148 | which is `18446744073709551615` and is more than anyone would wish to wait for! | |
149 | ||
150 | So it may be wise to chose some reasonable estimate of how many iterations may be needed, | |
151 | In this case the function is so well behaved that we can chose a low value of 20. | |
152 | ||
153 | Internally when Boost.Math uses these functions, it sets the maximum iterations to | |
154 | `policies::get_max_root_iterations<Policy>();`.] | |
155 | ||
156 | Should we have wished we can show how many iterations were used in `bracket_and_solve_root` | |
157 | (this information is lost outside `cbrt_noderiv`), for example with: | |
158 | ||
159 | if (it >= maxit) | |
160 | { | |
161 | std::cout << "Unable to locate solution in " << maxit << " iterations:" | |
162 | " Current best guess is between " << r.first << " and " << r.second << std::endl; | |
163 | } | |
164 | else | |
165 | { | |
166 | std::cout << "Converged after " << it << " (from maximum of " << maxit << " iterations)." << std::endl; | |
167 | } | |
168 | ||
169 | for output like | |
170 | ||
171 | Converged after 11 (from maximum of 20 iterations). | |
172 | */ | |
173 | //] [/root_finding_noderiv_2] | |
174 | ||
175 | ||
176 | // Cube root with 1st derivative (slope) | |
177 | ||
178 | /* | |
179 | We now solve the same problem, but using more information about the function, | |
180 | to show how this can speed up finding the best estimate of the root. | |
181 | ||
182 | For the root function, the 1st differential (the slope of the tangent to a curve at any point) is known. | |
183 | ||
184 | If you need some reminders then | |
185 | [@http://en.wikipedia.org/wiki/Derivative#Derivatives_of_elementary_functions Derivatives of elementary functions] | |
186 | may help. | |
187 | ||
188 | Using the rule that the derivative of ['x[super n]] for positive n (actually all nonzero n) is ['n x[super n-1]], | |
189 | allows us to get the 1st differential as ['3x[super 2]]. | |
190 | ||
191 | To see how this extra information is used to find a root, view | |
192 | [@http://en.wikipedia.org/wiki/Newton%27s_method Newton-Raphson iterations] | |
193 | and the [@http://en.wikipedia.org/wiki/Newton%27s_method#mediaviewer/File:NewtonIteration_Ani.gif animation]. | |
194 | ||
195 | We need to define a different functor `cbrt_functor_deriv` that returns | |
196 | both the evaluation of the function to solve, along with its first derivative: | |
197 | ||
198 | To \'return\' two values, we use a `std::pair` of floating-point values | |
199 | (though we could equally have used a std::tuple): | |
200 | */ | |
201 | ||
202 | //[root_finding_1_deriv_1 | |
203 | ||
204 | template <class T> | |
205 | struct cbrt_functor_deriv | |
206 | { // Functor also returning 1st derivative. | |
207 | cbrt_functor_deriv(T const& to_find_root_of) : a(to_find_root_of) | |
208 | { // Constructor stores value a to find root of, | |
209 | // for example: calling cbrt_functor_deriv<T>(a) to use to get cube root of a. | |
210 | } | |
211 | std::pair<T, T> operator()(T const& x) | |
212 | { | |
213 | // Return both f(x) and f'(x). | |
214 | T fx = x*x*x - a; // Difference (estimate x^3 - value). | |
215 | T dx = 3 * x*x; // 1st derivative = 3x^2. | |
216 | return std::make_pair(fx, dx); // 'return' both fx and dx. | |
217 | } | |
218 | private: | |
219 | T a; // Store value to be 'cube_rooted'. | |
220 | }; | |
221 | ||
222 | /*`Our cube root function is now:*/ | |
223 | ||
224 | template <class T> | |
225 | T cbrt_deriv(T x) | |
226 | { | |
227 | // return cube root of x using 1st derivative and Newton_Raphson. | |
228 | using namespace boost::math::tools; | |
229 | int exponent; | |
230 | frexp(x, &exponent); // Get exponent of z (ignore mantissa). | |
231 | T guess = ldexp(1., exponent/3); // Rough guess is to divide the exponent by three. | |
232 | T min = ldexp(0.5, exponent/3); // Minimum possible value is half our guess. | |
233 | T max = ldexp(2., exponent/3); // Maximum possible value is twice our guess. | |
234 | const int digits = std::numeric_limits<T>::digits; // Maximum possible binary digits accuracy for type T. | |
235 | int get_digits = static_cast<int>(digits * 0.6); // Accuracy doubles with each step, so stop when we have | |
236 | // just over half the digits correct. | |
237 | const boost::uintmax_t maxit = 20; | |
238 | boost::uintmax_t it = maxit; | |
239 | T result = newton_raphson_iterate(cbrt_functor_deriv<T>(x), guess, min, max, get_digits, it); | |
240 | return result; | |
241 | } | |
242 | ||
243 | //] [/root_finding_1_deriv_1] | |
244 | ||
245 | ||
246 | /* | |
247 | [h3:cbrt_2_derivatives Cube root with 1st & 2nd derivative (slope & curvature)] | |
248 | ||
249 | Finally we define yet another functor `cbrt_functor_2deriv` that returns | |
250 | both the evaluation of the function to solve, | |
251 | along with its first *and second* derivatives: | |
252 | ||
253 | __spaces[''f](x) = 6x | |
254 | ||
255 | To \'return\' three values, we use a `tuple` of three floating-point values: | |
256 | */ | |
257 | ||
258 | //[root_finding_2deriv_1 | |
259 | ||
260 | template <class T> | |
261 | struct cbrt_functor_2deriv | |
262 | { | |
263 | // Functor returning both 1st and 2nd derivatives. | |
264 | cbrt_functor_2deriv(T const& to_find_root_of) : a(to_find_root_of) | |
265 | { // Constructor stores value a to find root of, for example: | |
266 | // calling cbrt_functor_2deriv<T>(x) to get cube root of x, | |
267 | } | |
268 | std::tuple<T, T, T> operator()(T const& x) | |
269 | { | |
270 | // Return both f(x) and f'(x) and f''(x). | |
271 | T fx = x*x*x - a; // Difference (estimate x^3 - value). | |
272 | T dx = 3 * x*x; // 1st derivative = 3x^2. | |
273 | T d2x = 6 * x; // 2nd derivative = 6x. | |
274 | return std::make_tuple(fx, dx, d2x); // 'return' fx, dx and d2x. | |
275 | } | |
276 | private: | |
277 | T a; // to be 'cube_rooted'. | |
278 | }; | |
279 | ||
280 | /*`Our cube root function is now:*/ | |
281 | ||
282 | template <class T> | |
283 | T cbrt_2deriv(T x) | |
284 | { | |
285 | // return cube root of x using 1st and 2nd derivatives and Halley. | |
286 | //using namespace std; // Help ADL of std functions. | |
287 | using namespace boost::math::tools; | |
288 | int exponent; | |
289 | frexp(x, &exponent); // Get exponent of z (ignore mantissa). | |
290 | T guess = ldexp(1., exponent/3); // Rough guess is to divide the exponent by three. | |
291 | T min = ldexp(0.5, exponent/3); // Minimum possible value is half our guess. | |
292 | T max = ldexp(2., exponent/3); // Maximum possible value is twice our guess. | |
293 | const int digits = std::numeric_limits<T>::digits; // Maximum possible binary digits accuracy for type T. | |
294 | // digits used to control how accurate to try to make the result. | |
295 | int get_digits = static_cast<int>(digits * 0.4); // Accuracy triples with each step, so stop when just | |
296 | // over one third of the digits are correct. | |
297 | boost::uintmax_t maxit = 20; | |
298 | T result = halley_iterate(cbrt_functor_2deriv<T>(x), guess, min, max, get_digits, maxit); | |
299 | return result; | |
300 | } | |
301 | ||
302 | //] [/root_finding_2deriv_1] | |
303 | ||
304 | //[root_finding_2deriv_lambda | |
305 | ||
306 | template <class T> | |
307 | T cbrt_2deriv_lambda(T x) | |
308 | { | |
309 | // return cube root of x using 1st and 2nd derivatives and Halley. | |
310 | //using namespace std; // Help ADL of std functions. | |
311 | using namespace boost::math::tools; | |
312 | int exponent; | |
313 | frexp(x, &exponent); // Get exponent of z (ignore mantissa). | |
314 | T guess = ldexp(1., exponent / 3); // Rough guess is to divide the exponent by three. | |
315 | T min = ldexp(0.5, exponent / 3); // Minimum possible value is half our guess. | |
316 | T max = ldexp(2., exponent / 3); // Maximum possible value is twice our guess. | |
317 | const int digits = std::numeric_limits<T>::digits; // Maximum possible binary digits accuracy for type T. | |
318 | // digits used to control how accurate to try to make the result. | |
319 | int get_digits = static_cast<int>(digits * 0.4); // Accuracy triples with each step, so stop when just | |
320 | // over one third of the digits are correct. | |
321 | boost::uintmax_t maxit = 20; | |
322 | T result = halley_iterate( | |
323 | // lambda function: | |
324 | [x](const T& g){ return std::make_tuple(g * g * g - x, 3 * g * g, 6 * g); }, | |
325 | guess, min, max, get_digits, maxit); | |
326 | return result; | |
327 | } | |
328 | ||
329 | //] [/root_finding_2deriv_lambda] | |
330 | /* | |
331 | ||
332 | [h3 Fifth-root function] | |
333 | Let's now suppose we want to find the [*fifth root] of a number ['a]. | |
334 | ||
335 | The equation we want to solve is : | |
336 | ||
337 | __spaces['f](x) = x[super 5] - a | |
338 | ||
339 | If your differentiation is a little rusty | |
340 | (or you are faced with an equation whose complexity is daunting), | |
341 | then you can get help, for example from the invaluable | |
342 | [@http://www.wolframalpha.com/ WolframAlpha site.] | |
343 | ||
f67539c2 | 344 | For example, entering the command: `differentiate x ^ 5` |
7c673cae FG |
345 | |
346 | or the Wolfram Language command: ` D[x ^ 5, x]` | |
347 | ||
348 | gives the output: `d/dx(x ^ 5) = 5 x ^ 4` | |
349 | ||
350 | and to get the second differential, enter: `second differentiate x ^ 5` | |
351 | ||
352 | or the Wolfram Language command: `D[x ^ 5, { x, 2 }]` | |
353 | ||
354 | to get the output: `d ^ 2 / dx ^ 2(x ^ 5) = 20 x ^ 3` | |
355 | ||
356 | To get a reference value, we can enter: [^fifth root 3126] | |
357 | ||
358 | or: `N[3126 ^ (1 / 5), 50]` | |
359 | ||
360 | to get a result with a precision of 50 decimal digits: | |
361 | ||
362 | 5.0003199590478625588206333405631053401128722314376 | |
363 | ||
364 | (We could also get a reference value using Boost.Multiprecision - see below). | |
365 | ||
366 | The 1st and 2nd derivatives of x[super 5] are: | |
367 | ||
368 | __spaces['f]\'(x) = 5x[super 4] | |
369 | ||
370 | __spaces['f]\'\'(x) = 20x[super 3] | |
371 | ||
372 | */ | |
373 | ||
374 | //[root_finding_fifth_1 | |
375 | //] [/root_finding_fifth_1] | |
376 | ||
377 | ||
378 | //[root_finding_fifth_functor_2deriv | |
379 | ||
380 | /*`Using these expressions for the derivatives, the functor is: | |
381 | */ | |
382 | ||
383 | template <class T> | |
384 | struct fifth_functor_2deriv | |
385 | { | |
386 | // Functor returning both 1st and 2nd derivatives. | |
387 | fifth_functor_2deriv(T const& to_find_root_of) : a(to_find_root_of) | |
388 | { /* Constructor stores value a to find root of, for example: */ } | |
389 | ||
390 | std::tuple<T, T, T> operator()(T const& x) | |
391 | { | |
392 | // Return both f(x) and f'(x) and f''(x). | |
393 | T fx = boost::math::pow<5>(x) - a; // Difference (estimate x^3 - value). | |
394 | T dx = 5 * boost::math::pow<4>(x); // 1st derivative = 5x^4. | |
395 | T d2x = 20 * boost::math::pow<3>(x); // 2nd derivative = 20 x^3 | |
396 | return std::make_tuple(fx, dx, d2x); // 'return' fx, dx and d2x. | |
397 | } | |
398 | private: | |
399 | T a; // to be 'fifth_rooted'. | |
400 | }; // struct fifth_functor_2deriv | |
401 | ||
402 | //] [/root_finding_fifth_functor_2deriv] | |
403 | ||
404 | //[root_finding_fifth_2deriv | |
405 | ||
406 | /*`Our fifth-root function is now: | |
407 | */ | |
408 | ||
409 | template <class T> | |
410 | T fifth_2deriv(T x) | |
411 | { | |
412 | // return fifth root of x using 1st and 2nd derivatives and Halley. | |
413 | using namespace std; // Help ADL of std functions. | |
414 | using namespace boost::math::tools; // for halley_iterate. | |
415 | ||
416 | int exponent; | |
417 | frexp(x, &exponent); // Get exponent of z (ignore mantissa). | |
418 | T guess = ldexp(1., exponent / 5); // Rough guess is to divide the exponent by five. | |
419 | T min = ldexp(0.5, exponent / 5); // Minimum possible value is half our guess. | |
420 | T max = ldexp(2., exponent / 5); // Maximum possible value is twice our guess. | |
421 | // Stop when slightly more than one of the digits are correct: | |
422 | const int digits = static_cast<int>(std::numeric_limits<T>::digits * 0.4); | |
423 | const boost::uintmax_t maxit = 50; | |
424 | boost::uintmax_t it = maxit; | |
425 | T result = halley_iterate(fifth_functor_2deriv<T>(x), guess, min, max, digits, it); | |
426 | return result; | |
427 | } | |
428 | ||
429 | //] [/root_finding_fifth_2deriv] | |
430 | ||
431 | ||
432 | int main() | |
433 | { | |
434 | std::cout << "Root finding Examples." << std::endl; | |
435 | std::cout.precision(std::numeric_limits<double>::max_digits10); | |
436 | // Show all possibly significant decimal digits for double. | |
437 | // std::cout.precision(std::numeric_limits<double>::digits10); | |
438 | // Show all guaranteed significant decimal digits for double. | |
439 | ||
440 | ||
441 | //[root_finding_main_1 | |
442 | try | |
443 | { | |
444 | double threecubed = 27.; // Value that has an *exactly representable* integer cube root. | |
445 | double threecubedp1 = 28.; // Value whose cube root is *not* exactly representable. | |
446 | ||
447 | std::cout << "cbrt(28) " << boost::math::cbrt(28.) << std::endl; // boost::math:: version of cbrt. | |
448 | std::cout << "std::cbrt(28) " << std::cbrt(28.) << std::endl; // std:: version of cbrt. | |
449 | std::cout <<" cast double " << static_cast<double>(3.0365889718756625194208095785056696355814539772481111) << std::endl; | |
450 | ||
451 | // Cube root using bracketing: | |
452 | double r = cbrt_noderiv(threecubed); | |
453 | std::cout << "cbrt_noderiv(" << threecubed << ") = " << r << std::endl; | |
454 | r = cbrt_noderiv(threecubedp1); | |
455 | std::cout << "cbrt_noderiv(" << threecubedp1 << ") = " << r << std::endl; | |
456 | //] [/root_finding_main_1] | |
457 | //[root_finding_main_2 | |
458 | ||
459 | // Cube root using 1st differential Newton-Raphson: | |
460 | r = cbrt_deriv(threecubed); | |
461 | std::cout << "cbrt_deriv(" << threecubed << ") = " << r << std::endl; | |
462 | r = cbrt_deriv(threecubedp1); | |
463 | std::cout << "cbrt_deriv(" << threecubedp1 << ") = " << r << std::endl; | |
464 | ||
465 | // Cube root using Halley with 1st and 2nd differentials. | |
466 | r = cbrt_2deriv(threecubed); | |
467 | std::cout << "cbrt_2deriv(" << threecubed << ") = " << r << std::endl; | |
468 | r = cbrt_2deriv(threecubedp1); | |
469 | std::cout << "cbrt_2deriv(" << threecubedp1 << ") = " << r << std::endl; | |
470 | ||
471 | // Cube root using lambda's: | |
472 | r = cbrt_2deriv_lambda(threecubed); | |
473 | std::cout << "cbrt_2deriv(" << threecubed << ") = " << r << std::endl; | |
474 | r = cbrt_2deriv_lambda(threecubedp1); | |
475 | std::cout << "cbrt_2deriv(" << threecubedp1 << ") = " << r << std::endl; | |
476 | ||
477 | // Fifth root. | |
478 | ||
479 | double fivepowfive = 3125; // Example of a value that has an exact integer fifth root. | |
480 | // Exact value of fifth root is exactly 5. | |
481 | std::cout << "Fifth root of " << fivepowfive << " is " << 5 << std::endl; | |
482 | ||
483 | double fivepowfivep1 = fivepowfive + 1; // Example of a value whose fifth root is *not* exactly representable. | |
484 | // Value of fifth root is 5.0003199590478625588206333405631053401128722314376 (50 decimal digits precision) | |
485 | // and to std::numeric_limits<double>::max_digits10 double precision (usually 17) is | |
486 | ||
487 | double root5v2 = static_cast<double>(5.0003199590478625588206333405631053401128722314376); | |
488 | std::cout << "Fifth root of " << fivepowfivep1 << " is " << root5v2 << std::endl; | |
489 | ||
490 | // Using Halley with 1st and 2nd differentials. | |
491 | r = fifth_2deriv(fivepowfive); | |
492 | std::cout << "fifth_2deriv(" << fivepowfive << ") = " << r << std::endl; | |
493 | r = fifth_2deriv(fivepowfivep1); | |
494 | std::cout << "fifth_2deriv(" << fivepowfivep1 << ") = " << r << std::endl; | |
495 | //] [/root_finding_main_?] | |
496 | } | |
497 | catch(const std::exception& e) | |
498 | { // Always useful to include try & catch blocks because default policies | |
499 | // are to throw exceptions on arguments that cause errors like underflow, overflow. | |
500 | // Lacking try & catch blocks, the program will abort without a message below, | |
501 | // which may give some helpful clues as to the cause of the exception. | |
502 | std::cout << | |
503 | "\n""Message from thrown exception was:\n " << e.what() << std::endl; | |
504 | } | |
505 | return 0; | |
506 | } // int main() | |
507 | ||
508 | //[root_finding_example_output | |
509 | /*` | |
510 | Normal output is: | |
511 | ||
512 | [pre | |
513 | root_finding_example.cpp | |
514 | Generating code | |
515 | Finished generating code | |
516 | root_finding_example.vcxproj -> J:\Cpp\MathToolkit\test\Math_test\Release\root_finding_example.exe | |
517 | Cube Root finding (cbrt) Example. | |
518 | Iterations 10 | |
519 | cbrt_1(27) = 3 | |
520 | Iterations 10 | |
521 | Unable to locate solution in chosen iterations: Current best guess is between 3.0365889718756613 and 3.0365889718756627 | |
522 | cbrt_1(28) = 3.0365889718756618 | |
523 | cbrt_1(27) = 3 | |
524 | cbrt_2(28) = 3.0365889718756627 | |
525 | Iterations 4 | |
526 | cbrt_3(27) = 3 | |
527 | Iterations 5 | |
528 | cbrt_3(28) = 3.0365889718756627 | |
529 | ||
530 | ] [/pre] | |
531 | ||
532 | to get some (much!) diagnostic output we can add | |
533 | ||
534 | #define BOOST_MATH_INSTRUMENT | |
535 | ||
536 | [pre | |
537 | ||
538 | ] | |
539 | */ | |
540 | //] [/root_finding_example_output] | |
541 | ||
542 | /* | |
543 | ||
544 | cbrt(28) 3.0365889718756622 | |
545 | std::cbrt(28) 3.0365889718756627 | |
546 | ||
547 | */ |