root_finding_example.cpp 20 KB

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