123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547 |
- // root_finding_example.cpp
- // Copyright Paul A. Bristow 2010, 2015
- // Use, modification and distribution are subject to the
- // Boost Software License, Version 1.0.
- // (See accompanying file LICENSE_1_0.txt
- // or copy at http://www.boost.org/LICENSE_1_0.txt)
- // Example of finding roots using Newton-Raphson, Halley.
- // Note that this file contains Quickbook mark-up as well as code
- // and comments, don't change any of the special comment mark-ups!
- //#define BOOST_MATH_INSTRUMENT
- /*
- This example demonstrates how to use the various tools for root finding
- taking the simple cube root function (`cbrt`) as an example.
- It shows how use of derivatives can improve the speed.
- (But is only a demonstration and does not try to make the ultimate improvements of 'real-life'
- implementation of `boost::math::cbrt`, mainly by using a better computed initial 'guess'
- at `<boost/math/special_functions/cbrt.hpp>`).
- Then we show how a higher root (fifth) can be computed,
- and in `root_finding_n_example.cpp` a generic method
- for the ['n]th root that constructs the derivatives at compile-time,
- These methods should be applicable to other functions that can be differentiated easily.
- First some `#includes` that will be needed.
- [tip For clarity, `using` statements are provided to list what functions are being used in this example:
- you can of course partly or fully qualify the names in other ways.
- (For your application, you may wish to extract some parts into header files,
- but you should never use `using` statements globally in header files).]
- */
- //[root_finding_include_1
- #include <boost/math/tools/roots.hpp>
- //using boost::math::policies::policy;
- //using boost::math::tools::newton_raphson_iterate;
- //using boost::math::tools::halley_iterate; //
- //using boost::math::tools::eps_tolerance; // Binary functor for specified number of bits.
- //using boost::math::tools::bracket_and_solve_root;
- //using boost::math::tools::toms748_solve;
- #include <boost/math/special_functions/next.hpp> // For float_distance.
- #include <tuple> // for std::tuple and std::make_tuple.
- #include <boost/math/special_functions/cbrt.hpp> // For boost::math::cbrt.
- //] [/root_finding_include_1]
- // using boost::math::tuple;
- // using boost::math::make_tuple;
- // using boost::math::tie;
- // which provide convenient aliases for various implementations,
- // including std::tr1, depending on what is available.
- #include <iostream>
- //using std::cout; using std::endl;
- #include <iomanip>
- //using std::setw; using std::setprecision;
- #include <limits>
- //using std::numeric_limits;
- /*
- Let's suppose we want to find the root of a number ['a], and to start, compute the cube root.
- So the equation we want to solve is:
- __spaces ['f](x) = x[cubed] - a
- We will first solve this without using any information
- about the slope or curvature of the cube root function.
- We then show how adding what we can know about this function, first just the slope,
- the 1st derivation /f'(x)/, will speed homing in on the solution.
- Lastly we show how adding the curvature /f''(x)/ too will speed convergence even more.
- */
- //[root_finding_noderiv_1
- template <class T>
- struct cbrt_functor_noderiv
- {
- // cube root of x using only function - no derivatives.
- cbrt_functor_noderiv(T const& to_find_root_of) : a(to_find_root_of)
- { /* Constructor just stores value a to find root of. */ }
- T operator()(T const& x)
- {
- T fx = x*x*x - a; // Difference (estimate x^3 - a).
- return fx;
- }
- private:
- T a; // to be 'cube_rooted'.
- };
- //] [/root_finding_noderiv_1
- /*
- Implementing the cube root function itself is fairly trivial now:
- the hardest part is finding a good approximation to begin with.
- In this case we'll just divide the exponent by three.
- (There are better but more complex guess algorithms used in 'real-life'.)
- Cube root function is 'Really Well Behaved' in that it is monotonic
- and has only one root (we leave negative values 'as an exercise for the student').
- */
- //[root_finding_noderiv_2
- template <class T>
- T cbrt_noderiv(T x)
- {
- // return cube root of x using bracket_and_solve (no derivatives).
- using namespace std; // Help ADL of std functions.
- using namespace boost::math::tools; // For bracket_and_solve_root.
- int exponent;
- frexp(x, &exponent); // Get exponent of z (ignore mantissa).
- T guess = ldexp(1., exponent/3); // Rough guess is to divide the exponent by three.
- T factor = 2; // How big steps to take when searching.
- const boost::uintmax_t maxit = 20; // Limit to maximum iterations.
- boost::uintmax_t it = maxit; // Initally our chosen max iterations, but updated with actual.
- bool is_rising = true; // So if result if guess^3 is too low, then try increasing guess.
- int digits = std::numeric_limits<T>::digits; // Maximum possible binary digits accuracy for type T.
- // Some fraction of digits is used to control how accurate to try to make the result.
- int get_digits = digits - 3; // We have to have a non-zero interval at each step, so
- // maximum accuracy is digits - 1. But we also have to
- // allow for inaccuracy in f(x), otherwise the last few
- // iterations just thrash around.
- eps_tolerance<T> tol(get_digits); // Set the tolerance.
- std::pair<T, T> r = bracket_and_solve_root(cbrt_functor_noderiv<T>(x), guess, factor, is_rising, tol, it);
- return r.first + (r.second - r.first)/2; // Midway between brackets is our result, if necessary we could
- // return the result as an interval here.
- }
- /*`
- [note The final parameter specifying a maximum number of iterations is optional.
- However, it defaults to `boost::uintmax_t maxit = (std::numeric_limits<boost::uintmax_t>::max)();`
- which is `18446744073709551615` and is more than anyone would wish to wait for!
- So it may be wise to chose some reasonable estimate of how many iterations may be needed,
- In this case the function is so well behaved that we can chose a low value of 20.
- Internally when Boost.Math uses these functions, it sets the maximum iterations to
- `policies::get_max_root_iterations<Policy>();`.]
- Should we have wished we can show how many iterations were used in `bracket_and_solve_root`
- (this information is lost outside `cbrt_noderiv`), for example with:
- if (it >= maxit)
- {
- std::cout << "Unable to locate solution in " << maxit << " iterations:"
- " Current best guess is between " << r.first << " and " << r.second << std::endl;
- }
- else
- {
- std::cout << "Converged after " << it << " (from maximum of " << maxit << " iterations)." << std::endl;
- }
- for output like
- Converged after 11 (from maximum of 20 iterations).
- */
- //] [/root_finding_noderiv_2]
- // Cube root with 1st derivative (slope)
- /*
- We now solve the same problem, but using more information about the function,
- to show how this can speed up finding the best estimate of the root.
- For the root function, the 1st differential (the slope of the tangent to a curve at any point) is known.
- If you need some reminders then
- [@http://en.wikipedia.org/wiki/Derivative#Derivatives_of_elementary_functions Derivatives of elementary functions]
- may help.
- Using the rule that the derivative of ['x[super n]] for positive n (actually all nonzero n) is ['n x[super n-1]],
- allows us to get the 1st differential as ['3x[super 2]].
- To see how this extra information is used to find a root, view
- [@http://en.wikipedia.org/wiki/Newton%27s_method Newton-Raphson iterations]
- and the [@http://en.wikipedia.org/wiki/Newton%27s_method#mediaviewer/File:NewtonIteration_Ani.gif animation].
- We need to define a different functor `cbrt_functor_deriv` that returns
- both the evaluation of the function to solve, along with its first derivative:
- To \'return\' two values, we use a `std::pair` of floating-point values
- (though we could equally have used a std::tuple):
- */
- //[root_finding_1_deriv_1
- template <class T>
- struct cbrt_functor_deriv
- { // Functor also returning 1st derivative.
- cbrt_functor_deriv(T const& to_find_root_of) : a(to_find_root_of)
- { // Constructor stores value a to find root of,
- // for example: calling cbrt_functor_deriv<T>(a) to use to get cube root of a.
- }
- std::pair<T, T> operator()(T const& x)
- {
- // Return both f(x) and f'(x).
- T fx = x*x*x - a; // Difference (estimate x^3 - value).
- T dx = 3 * x*x; // 1st derivative = 3x^2.
- return std::make_pair(fx, dx); // 'return' both fx and dx.
- }
- private:
- T a; // Store value to be 'cube_rooted'.
- };
- /*`Our cube root function is now:*/
- template <class T>
- T cbrt_deriv(T x)
- {
- // return cube root of x using 1st derivative and Newton_Raphson.
- using namespace boost::math::tools;
- int exponent;
- frexp(x, &exponent); // Get exponent of z (ignore mantissa).
- T guess = ldexp(1., exponent/3); // Rough guess is to divide the exponent by three.
- T min = ldexp(0.5, exponent/3); // Minimum possible value is half our guess.
- T max = ldexp(2., exponent/3); // Maximum possible value is twice our guess.
- const int digits = std::numeric_limits<T>::digits; // Maximum possible binary digits accuracy for type T.
- int get_digits = static_cast<int>(digits * 0.6); // Accuracy doubles with each step, so stop when we have
- // just over half the digits correct.
- const boost::uintmax_t maxit = 20;
- boost::uintmax_t it = maxit;
- T result = newton_raphson_iterate(cbrt_functor_deriv<T>(x), guess, min, max, get_digits, it);
- return result;
- }
- //] [/root_finding_1_deriv_1]
- /*
- [h3:cbrt_2_derivatives Cube root with 1st & 2nd derivative (slope & curvature)]
- Finally we define yet another functor `cbrt_functor_2deriv` that returns
- both the evaluation of the function to solve,
- along with its first *and second* derivatives:
- __spaces[''f](x) = 6x
- To \'return\' three values, we use a `tuple` of three floating-point values:
- */
- //[root_finding_2deriv_1
- template <class T>
- struct cbrt_functor_2deriv
- {
- // Functor returning both 1st and 2nd derivatives.
- cbrt_functor_2deriv(T const& to_find_root_of) : a(to_find_root_of)
- { // Constructor stores value a to find root of, for example:
- // calling cbrt_functor_2deriv<T>(x) to get cube root of x,
- }
- std::tuple<T, T, T> operator()(T const& x)
- {
- // Return both f(x) and f'(x) and f''(x).
- T fx = x*x*x - a; // Difference (estimate x^3 - value).
- T dx = 3 * x*x; // 1st derivative = 3x^2.
- T d2x = 6 * x; // 2nd derivative = 6x.
- return std::make_tuple(fx, dx, d2x); // 'return' fx, dx and d2x.
- }
- private:
- T a; // to be 'cube_rooted'.
- };
- /*`Our cube root function is now:*/
- template <class T>
- T cbrt_2deriv(T x)
- {
- // return cube root of x using 1st and 2nd derivatives and Halley.
- //using namespace std; // Help ADL of std functions.
- using namespace boost::math::tools;
- int exponent;
- frexp(x, &exponent); // Get exponent of z (ignore mantissa).
- T guess = ldexp(1., exponent/3); // Rough guess is to divide the exponent by three.
- T min = ldexp(0.5, exponent/3); // Minimum possible value is half our guess.
- T max = ldexp(2., exponent/3); // Maximum possible value is twice our guess.
- const int digits = std::numeric_limits<T>::digits; // Maximum possible binary digits accuracy for type T.
- // digits used to control how accurate to try to make the result.
- int get_digits = static_cast<int>(digits * 0.4); // Accuracy triples with each step, so stop when just
- // over one third of the digits are correct.
- boost::uintmax_t maxit = 20;
- T result = halley_iterate(cbrt_functor_2deriv<T>(x), guess, min, max, get_digits, maxit);
- return result;
- }
- //] [/root_finding_2deriv_1]
- //[root_finding_2deriv_lambda
- template <class T>
- T cbrt_2deriv_lambda(T x)
- {
- // return cube root of x using 1st and 2nd derivatives and Halley.
- //using namespace std; // Help ADL of std functions.
- using namespace boost::math::tools;
- int exponent;
- frexp(x, &exponent); // Get exponent of z (ignore mantissa).
- T guess = ldexp(1., exponent / 3); // Rough guess is to divide the exponent by three.
- T min = ldexp(0.5, exponent / 3); // Minimum possible value is half our guess.
- T max = ldexp(2., exponent / 3); // Maximum possible value is twice our guess.
- const int digits = std::numeric_limits<T>::digits; // Maximum possible binary digits accuracy for type T.
- // digits used to control how accurate to try to make the result.
- int get_digits = static_cast<int>(digits * 0.4); // Accuracy triples with each step, so stop when just
- // over one third of the digits are correct.
- boost::uintmax_t maxit = 20;
- T result = halley_iterate(
- // lambda function:
- [x](const T& g){ return std::make_tuple(g * g * g - x, 3 * g * g, 6 * g); },
- guess, min, max, get_digits, maxit);
- return result;
- }
- //] [/root_finding_2deriv_lambda]
- /*
- [h3 Fifth-root function]
- Let's now suppose we want to find the [*fifth root] of a number ['a].
- The equation we want to solve is :
- __spaces['f](x) = x[super 5] - a
- If your differentiation is a little rusty
- (or you are faced with an equation whose complexity is daunting),
- then you can get help, for example from the invaluable
- [@http://www.wolframalpha.com/ WolframAlpha site.]
- For example, entering the commmand: `differentiate x ^ 5`
- or the Wolfram Language command: ` D[x ^ 5, x]`
- gives the output: `d/dx(x ^ 5) = 5 x ^ 4`
- and to get the second differential, enter: `second differentiate x ^ 5`
- or the Wolfram Language command: `D[x ^ 5, { x, 2 }]`
- to get the output: `d ^ 2 / dx ^ 2(x ^ 5) = 20 x ^ 3`
- To get a reference value, we can enter: [^fifth root 3126]
- or: `N[3126 ^ (1 / 5), 50]`
- to get a result with a precision of 50 decimal digits:
- 5.0003199590478625588206333405631053401128722314376
- (We could also get a reference value using Boost.Multiprecision - see below).
- The 1st and 2nd derivatives of x[super 5] are:
- __spaces['f]\'(x) = 5x[super 4]
- __spaces['f]\'\'(x) = 20x[super 3]
- */
- //[root_finding_fifth_1
- //] [/root_finding_fifth_1]
- //[root_finding_fifth_functor_2deriv
- /*`Using these expressions for the derivatives, the functor is:
- */
- template <class T>
- struct fifth_functor_2deriv
- {
- // Functor returning both 1st and 2nd derivatives.
- fifth_functor_2deriv(T const& to_find_root_of) : a(to_find_root_of)
- { /* Constructor stores value a to find root of, for example: */ }
- std::tuple<T, T, T> operator()(T const& x)
- {
- // Return both f(x) and f'(x) and f''(x).
- T fx = boost::math::pow<5>(x) - a; // Difference (estimate x^3 - value).
- T dx = 5 * boost::math::pow<4>(x); // 1st derivative = 5x^4.
- T d2x = 20 * boost::math::pow<3>(x); // 2nd derivative = 20 x^3
- return std::make_tuple(fx, dx, d2x); // 'return' fx, dx and d2x.
- }
- private:
- T a; // to be 'fifth_rooted'.
- }; // struct fifth_functor_2deriv
- //] [/root_finding_fifth_functor_2deriv]
- //[root_finding_fifth_2deriv
- /*`Our fifth-root function is now:
- */
- template <class T>
- T fifth_2deriv(T x)
- {
- // return fifth root of x using 1st and 2nd derivatives and Halley.
- using namespace std; // Help ADL of std functions.
- using namespace boost::math::tools; // for halley_iterate.
- int exponent;
- frexp(x, &exponent); // Get exponent of z (ignore mantissa).
- T guess = ldexp(1., exponent / 5); // Rough guess is to divide the exponent by five.
- T min = ldexp(0.5, exponent / 5); // Minimum possible value is half our guess.
- T max = ldexp(2., exponent / 5); // Maximum possible value is twice our guess.
- // Stop when slightly more than one of the digits are correct:
- const int digits = static_cast<int>(std::numeric_limits<T>::digits * 0.4);
- const boost::uintmax_t maxit = 50;
- boost::uintmax_t it = maxit;
- T result = halley_iterate(fifth_functor_2deriv<T>(x), guess, min, max, digits, it);
- return result;
- }
- //] [/root_finding_fifth_2deriv]
- int main()
- {
- std::cout << "Root finding Examples." << std::endl;
- std::cout.precision(std::numeric_limits<double>::max_digits10);
- // Show all possibly significant decimal digits for double.
- // std::cout.precision(std::numeric_limits<double>::digits10);
- // Show all guaranteed significant decimal digits for double.
- //[root_finding_main_1
- try
- {
- double threecubed = 27.; // Value that has an *exactly representable* integer cube root.
- double threecubedp1 = 28.; // Value whose cube root is *not* exactly representable.
- std::cout << "cbrt(28) " << boost::math::cbrt(28.) << std::endl; // boost::math:: version of cbrt.
- std::cout << "std::cbrt(28) " << std::cbrt(28.) << std::endl; // std:: version of cbrt.
- std::cout <<" cast double " << static_cast<double>(3.0365889718756625194208095785056696355814539772481111) << std::endl;
- // Cube root using bracketing:
- double r = cbrt_noderiv(threecubed);
- std::cout << "cbrt_noderiv(" << threecubed << ") = " << r << std::endl;
- r = cbrt_noderiv(threecubedp1);
- std::cout << "cbrt_noderiv(" << threecubedp1 << ") = " << r << std::endl;
- //] [/root_finding_main_1]
- //[root_finding_main_2
- // Cube root using 1st differential Newton-Raphson:
- r = cbrt_deriv(threecubed);
- std::cout << "cbrt_deriv(" << threecubed << ") = " << r << std::endl;
- r = cbrt_deriv(threecubedp1);
- std::cout << "cbrt_deriv(" << threecubedp1 << ") = " << r << std::endl;
- // Cube root using Halley with 1st and 2nd differentials.
- r = cbrt_2deriv(threecubed);
- std::cout << "cbrt_2deriv(" << threecubed << ") = " << r << std::endl;
- r = cbrt_2deriv(threecubedp1);
- std::cout << "cbrt_2deriv(" << threecubedp1 << ") = " << r << std::endl;
- // Cube root using lambda's:
- r = cbrt_2deriv_lambda(threecubed);
- std::cout << "cbrt_2deriv(" << threecubed << ") = " << r << std::endl;
- r = cbrt_2deriv_lambda(threecubedp1);
- std::cout << "cbrt_2deriv(" << threecubedp1 << ") = " << r << std::endl;
- // Fifth root.
- double fivepowfive = 3125; // Example of a value that has an exact integer fifth root.
- // Exact value of fifth root is exactly 5.
- std::cout << "Fifth root of " << fivepowfive << " is " << 5 << std::endl;
- double fivepowfivep1 = fivepowfive + 1; // Example of a value whose fifth root is *not* exactly representable.
- // Value of fifth root is 5.0003199590478625588206333405631053401128722314376 (50 decimal digits precision)
- // and to std::numeric_limits<double>::max_digits10 double precision (usually 17) is
- double root5v2 = static_cast<double>(5.0003199590478625588206333405631053401128722314376);
- std::cout << "Fifth root of " << fivepowfivep1 << " is " << root5v2 << std::endl;
- // Using Halley with 1st and 2nd differentials.
- r = fifth_2deriv(fivepowfive);
- std::cout << "fifth_2deriv(" << fivepowfive << ") = " << r << std::endl;
- r = fifth_2deriv(fivepowfivep1);
- std::cout << "fifth_2deriv(" << fivepowfivep1 << ") = " << r << std::endl;
- //] [/root_finding_main_?]
- }
- catch(const std::exception& e)
- { // Always useful to include try & catch blocks because default policies
- // are to throw exceptions on arguments that cause errors like underflow, overflow.
- // Lacking try & catch blocks, the program will abort without a message below,
- // which may give some helpful clues as to the cause of the exception.
- std::cout <<
- "\n""Message from thrown exception was:\n " << e.what() << std::endl;
- }
- return 0;
- } // int main()
- //[root_finding_example_output
- /*`
- Normal output is:
- [pre
- root_finding_example.cpp
- Generating code
- Finished generating code
- root_finding_example.vcxproj -> J:\Cpp\MathToolkit\test\Math_test\Release\root_finding_example.exe
- Cube Root finding (cbrt) Example.
- Iterations 10
- cbrt_1(27) = 3
- Iterations 10
- Unable to locate solution in chosen iterations: Current best guess is between 3.0365889718756613 and 3.0365889718756627
- cbrt_1(28) = 3.0365889718756618
- cbrt_1(27) = 3
- cbrt_2(28) = 3.0365889718756627
- Iterations 4
- cbrt_3(27) = 3
- Iterations 5
- cbrt_3(28) = 3.0365889718756627
- ] [/pre]
- to get some (much!) diagnostic output we can add
- #define BOOST_MATH_INSTRUMENT
- [pre
- ]
- */
- //] [/root_finding_example_output]
- /*
- cbrt(28) 3.0365889718756622
- std::cbrt(28) 3.0365889718756627
- */
|