123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485 |
- // root_finding_fith.cpp
- // Copyright Paul A. Bristow 2014.
- // 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 fifth root using Newton-Raphson, Halley, Schroder, TOMS748 .
- // Note that this file contains Quickbook mark-up as well as code
- // and comments, don't change any of the special comment mark-ups!
- // To get (copious!) diagnostic output, add make this define here or elsewhere.
- //#define BOOST_MATH_INSTRUMENT
- //[root_fifth_headers
- /*
- This example demonstrates how to use the Boost.Math tools for root finding,
- taking the fifth root function (fifth_root) as an example.
- It shows how use of derivatives can improve the speed.
- First some includes that will be needed.
- Using statements are provided to list what functions are being used in this example:
- you can of course qualify the names in other ways.
- */
- #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>
- #include <tuple>
- #include <utility> // pair, make_pair
- //] [/root_finding_headers]
- #include <iostream>
- using std::cout; using std::endl;
- #include <iomanip>
- using std::setw; using std::setprecision;
- #include <limits>
- using std::numeric_limits;
- /*
- //[root_finding_fifth_1
- Let's suppose we want to find the fifth root of a number.
- The equation we want to solve is:
- __spaces ['f](x) = x[fifth]
- We will first solve this without using any information
- about the slope or curvature of the fifth function.
- 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/ site
- 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
- D[x^5, {x, 2}]
- to get the output
- d^2/dx^2(x^5) = 20 x^3
- or
- 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).
- We then show how adding what we can know, for this function, about the slope,
- the 1st derivation /f'(x)/, will speed homing in on the solution,
- and then finally how adding the curvature /f''(x)/ as well will improve even more.
- The 1st and 2nd derivatives of x[fifth] are:
- __spaces ['f]\'(x) = 2x[sup2]
- __spaces ['f]\'\'(x) = 6x
- */
- //] [/root_finding_fifth_1]
- //[root_finding_fifth_functor_noderiv
- template <class T>
- struct fifth_functor_noderiv
- { // fifth root of x using only function - no derivatives.
- fifth_functor_noderiv(T const& to_find_root_of) : value(to_find_root_of)
- { // Constructor stores value to find root of.
- // For example: calling fifth_functor<T>(x) to get fifth root of x.
- }
- T operator()(T const& x)
- { //! \returns f(x) - value.
- T fx = x*x*x*x*x - value; // Difference (estimate x^5 - value).
- return fx;
- }
- private:
- T value; // to be 'fifth_rooted'.
- };
- //] [/root_finding_fifth_functor_noderiv]
- //cout << ", std::numeric_limits<" << typeid(T).name() << ">::digits = " << digits
- // << ", accuracy " << get_digits << " bits."<< endl;
- /*`Implementing the fifth 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 five.
- (There are better but more complex guess algorithms used in 'real-life'.)
- fifth 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_fifth_noderiv
- template <class T>
- T fifth_noderiv(T x)
- { //! \returns fifth 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 / 5); // Rough guess is to divide the exponent by five.
- T factor = 2; // To multiply and divide guess to bracket.
- // digits used to control how accurate to try to make the result.
- // int digits = 3 * std::numeric_limits<T>::digits / 4; // 3/4 maximum possible binary digits accuracy for type T.
- int digits = std::numeric_limits<T>::digits; // Maximum possible binary digits accuracy for type T.
- //boost::uintmax_t maxit = (std::numeric_limits<boost::uintmax_t>::max)();
- // (std::numeric_limits<boost::uintmax_t>::max)() = 18446744073709551615
- // which is more than anyone might wish to wait for!!!
- // so better to choose some reasonable estimate of how many iterations may be needed.
- const boost::uintmax_t maxit = 50; // Chosen max iterations,
- // but updated on exit with actual iteration count.
- // We could also have used a maximum iterations provided by any policy:
- // boost::uintmax_t max_it = policies::get_max_root_iterations<Policy>();
- boost::uintmax_t it = maxit; // Initally our chosen max iterations,
- bool is_rising = true; // So if result if guess^5 is too low, try increasing guess.
- eps_tolerance<double> tol(digits);
- std::pair<T, T> r =
- bracket_and_solve_root(fifth_functor_noderiv<T>(x), guess, factor, is_rising, tol, it);
- // because the iteration count is updating,
- // you can't call with a literal maximum iterations value thus:
- //bracket_and_solve_root(fifth_functor_noderiv<T>(x), guess, factor, is_rising, tol, 20);
- // Can show how many iterations (this information is lost outside fifth_noderiv).
- cout << "Iterations " << it << endl;
- if (it >= maxit)
- { // Failed to converge (or is jumping between bracket values).
- cout << "Unable to locate solution in chosen iterations:"
- " Current best guess is between " << r.first << " and " << r.second << endl;
- }
- T distance = float_distance(r.first, r.second);
- if (distance > 0)
- { //
- std::cout << distance << " bits separate the bracketing values." << std::endl;
- for (int i = 0; i < distance; i++)
- { // Show all the values within the bracketing values.
- std::cout << float_advance(r.first, i) << std::endl;
- }
- }
- else
- { // distance == 0 and r.second == r.first
- std::cout << "Converged to a single value " << r.first << std::endl;
- }
- return r.first + (r.second - r.first) / 2; // return midway between bracketed interval.
- } // T fifth_noderiv(T x)
- //] [/root_finding_fifth_noderiv]
- // maxit = 10
- // Unable to locate solution in chosen iterations: Current best guess is between 3.0365889718756613 and 3.0365889718756627
- /*`
- 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 this function, the 1st differential (the slope of the tangent to a curve at any point) is known.
- [@http://en.wikipedia.org/wiki/Derivative#Derivatives_of_elementary_functions Derivatives]
- gives some reminders.
- Using the rule that the derivative of x^n for positive n (actually all nonzero n) is nx^n-1,
- allows use to get the 1st differential as 3x^2.
- To see how this extra information is used to find the root, view this demo:
- [@http://en.wikipedia.org/wiki/Newton%27s_methodNewton Newton-Raphson iterations].
- We need to define a different functor that returns
- both the evaluation of the function to solve, along with its first derivative:
- To \'return\' two values, we use a pair of floating-point values:
- */
- //[root_finding_fifth_functor_1stderiv
- template <class T>
- struct fifth_functor_1stderiv
- { // Functor returning function and 1st derivative.
- fifth_functor_1stderiv(T const& target) : value(target)
- { // Constructor stores the value to be 'fifth_rooted'.
- }
- std::pair<T, T> operator()(T const& z) // z is best estimate so far.
- { // Return both f(x) and first derivative f'(x).
- T fx = z*z*z*z*z - value; // Difference estimate fx = x^5 - value.
- T d1x = 5 * z*z*z*z; // 1st derivative d1x = 5x^4.
- return std::make_pair(fx, d1x); // 'return' both fx and d1x.
- }
- private:
- T value; // to be 'fifth_rooted'.
- }; // fifth_functor_1stderiv
- //] [/root_finding_fifth_functor_1stderiv]
- /*`Our fifth root function using fifth_functor_1stderiv is now:*/
- //[root_finding_fifth_1deriv
- template <class T>
- T fifth_1deriv(T x)
- { //! \return fifth root of x using 1st derivative and Newton_Raphson.
- using namespace std; // For frexp, ldexp, numeric_limits.
- using namespace boost::math::tools; // For newton_raphson_iterate.
- int exponent;
- frexp(x, &exponent); // Get exponent of x (ignore mantissa).
- T guess = ldexp(1., exponent / 5); // Rough guess is to divide the exponent by three.
- // Set an initial bracket interval.
- 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.
- // digits used to control how accurate to try to make the result.
- int digits = std::numeric_limits<T>::digits; // Maximum possible binary digits accuracy for type T.
- const boost::uintmax_t maxit = 20; // Optionally limit the number of iterations.
- boost::uintmax_t it = maxit; // limit the number of iterations.
- //cout << "Max Iterations " << maxit << endl; //
- T result = newton_raphson_iterate(fifth_functor_1stderiv<T>(x), guess, min, max, digits, it);
- // Can check and show how many iterations (updated by newton_raphson_iterate).
- cout << it << " iterations (from max of " << maxit << ")" << endl;
- return result;
- } // fifth_1deriv
- //] [/root_finding_fifth_1deriv]
- // int get_digits = (digits * 2) /3; // Two thirds of maximum possible accuracy.
- //boost::uintmax_t maxit = (std::numeric_limits<boost::uintmax_t>::max)();
- // the default (std::numeric_limits<boost::uintmax_t>::max)() = 18446744073709551615
- // which is more than we might wish to wait for!!! so we can reduce it
- /*`
- Finally need to define yet another functor that returns
- both the evaluation of the function to solve,
- along with its first and second derivatives:
- f''(x) = 3 * 3x
- To \'return\' three values, we use a tuple of three floating-point values:
- */
- //[root_finding_fifth_functor_2deriv
- template <class T>
- struct fifth_functor_2deriv
- { // Functor returning both 1st and 2nd derivatives.
- fifth_functor_2deriv(T const& to_find_root_of) : value(to_find_root_of)
- { // Constructor stores value to find root of, for example:
- }
- // using boost::math::tuple; // to return three values.
- std::tuple<T, T, T> operator()(T const& x)
- { // Return both f(x) and f'(x) and f''(x).
- T fx = x*x*x*x*x - value; // Difference (estimate x^3 - value).
- T dx = 5 * x*x*x*x; // 1st derivative = 5x^4.
- T d2x = 20 * x*x*x; // 2nd derivative = 20 x^3
- return std::make_tuple(fx, dx, d2x); // 'return' fx, dx and d2x.
- }
- private:
- T value; // to be 'fifth_rooted'.
- }; // struct fifth_functor_2deriv
- //] [/root_finding_fifth_functor_2deriv]
- /*`Our fifth function is now:*/
- //[root_finding_fifth_2deriv
- 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; // 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 three.
- 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.
- int digits = std::numeric_limits<T>::digits / 2; // Half maximum possible binary digits accuracy for type T.
- 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);
- // Can show how many iterations (updated by halley_iterate).
- cout << it << " iterations (from max of " << maxit << ")" << endl;
- return result;
- } // fifth_2deriv(x)
- //] [/root_finding_fifth_2deriv]
- int main()
- {
- //[root_finding_example_1
- cout << "fifth Root finding (fifth) Example." << endl;
- // Show all possibly significant decimal digits.
- cout.precision(std::numeric_limits<double>::max_digits10);
- // or use cout.precision(max_digits10 = 2 + std::numeric_limits<double>::digits * 3010/10000);
- try
- { // Always use try'n'catch blocks with Boost.Math to get any error messages.
- double v27 = 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 " << v27 << " is " << 5 << std::endl;
- double v28 = v27+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 " << v28 << " is " << root5v2 << std::endl;
- // Using bracketing:
- double r = fifth_noderiv(v27);
- cout << "fifth_noderiv(" << v27 << ") = " << r << endl;
- r = fifth_noderiv(v28);
- cout << "fifth_noderiv(" << v28 << ") = " << r << endl;
- // Using 1st differential Newton-Raphson:
- r = fifth_1deriv(v27);
- cout << "fifth_1deriv(" << v27 << ") = " << r << endl;
- r = fifth_1deriv(v28);
- cout << "fifth_1deriv(" << v28 << ") = " << r << endl;
- // Using Halley with 1st and 2nd differentials.
- r = fifth_2deriv(v27);
- cout << "fifth_2deriv(" << v27 << ") = " << r << endl;
- r = fifth_2deriv(v28);
- cout << "fifth_2deriv(" << v28 << ") = " << r << endl;
- }
- 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;
- }
- //] [/root_finding_example_1
- return 0;
- } // int main()
- //[root_finding_example_output
- /*`
- Normal output is:
- [pre
- 1> Description: Autorun "J:\Cpp\MathToolkit\test\Math_test\Release\root_finding_fifth.exe"
- 1> fifth Root finding (fifth) Example.
- 1> Fifth root of 3125 is 5
- 1> Fifth root of 3126 is 5.0003199590478626
- 1> Iterations 10
- 1> Converged to a single value 5
- 1> fifth_noderiv(3125) = 5
- 1> Iterations 11
- 1> 2 bits separate the bracketing values.
- 1> 5.0003199590478609
- 1> 5.0003199590478618
- 1> fifth_noderiv(3126) = 5.0003199590478618
- 1> 6 iterations (from max of 20)
- 1> fifth_1deriv(3125) = 5
- 1> 7 iterations (from max of 20)
- 1> fifth_1deriv(3126) = 5.0003199590478626
- 1> 4 iterations (from max of 50)
- 1> fifth_2deriv(3125) = 5
- 1> 4 iterations (from max of 50)
- 1> fifth_2deriv(3126) = 5.0003199590478626
- [/pre]
- to get some (much!) diagnostic output we can add
- #define BOOST_MATH_INSTRUMENT
- [pre
- 1> fifth Root finding (fifth) Example.
- 1> I:\modular-boost\boost/math/tools/toms748_solve.hpp:537 a = 4 b = 8 fa = -2101 fb = 29643 count = 18
- 1> I:\modular-boost\boost/math/tools/toms748_solve.hpp:340 a = 4.264742943548387 b = 8
- 1> I:\modular-boost\boost/math/tools/toms748_solve.hpp:352 a = 4.264742943548387 b = 5.1409225585147951
- 1> I:\modular-boost\boost/math/tools/toms748_solve.hpp:259 a = 4.264742943548387 b = 5.1409225585147951 d = 8 e = 4 fa = -1714.2037505671719 fb = 465.91652114644285 fd = 29643 fe = -2101
- 1> I:\modular-boost\boost/math/tools/toms748_solve.hpp:267 q11 = -3.735257056451613 q21 = -0.045655399937094755 q31 = 0.68893005658139972 d21 = -2.9047328414222999 d31 = -0.18724955838500826
- 1> I:\modular-boost\boost/math/tools/toms748_solve.hpp:275 q22 = -0.15074699539567221 q32 = 0.007740525571111408 d32 = -0.13385363287680208 q33 = 0.074868009790687237 c = 5.0362815354915851
- 1> I:\modular-boost\boost/math/tools/toms748_solve.hpp:388 a = 4.264742943548387 b = 5.0362815354915851
- 1> I:\modular-boost\boost/math/tools/toms748_solve.hpp:259 a = 4.264742943548387 b = 5.0362815354915851 d = 5.1409225585147951 e = 8 fa = -1714.2037505671719 fb = 115.03721886368339 fd = 465.91652114644285 fe = 29643
- 1> I:\modular-boost\boost/math/tools/toms748_solve.hpp:267 q11 = -0.045655399937094755 q21 = -0.034306988726112195 q31 = 0.7230181097615842 d21 = -0.1389480117493222 d31 = -0.048520482181613811
- 1> I:\modular-boost\boost/math/tools/toms748_solve.hpp:275 q22 = -0.00036345624935100459 q32 = 0.011175908093791367 d32 = -0.0030375853617102483 q33 = 0.00014618657296010219 c = 4.999083147976723
- 1> I:\modular-boost\boost/math/tools/toms748_solve.hpp:408 a = 4.999083147976723 b = 5.0362815354915851
- 1> I:\modular-boost\boost/math/tools/toms748_solve.hpp:433 a = 4.999083147976723 b = 5.0008904277935091
- 1> I:\modular-boost\boost/math/tools/toms748_solve.hpp:434 tol = -0.00036152225583956088
- 1> I:\modular-boost\boost/math/tools/toms748_solve.hpp:259 a = 4.999083147976723 b = 5.0008904277935091 d = 5.0362815354915851 e = 4.264742943548387 fa = -2.8641119933622576 fb = 2.7835781082976609 fd = 115.03721886368339 fe = -1714.2037505671719
- 1> I:\modular-boost\boost/math/tools/toms748_solve.hpp:267 q11 = -0.048520482181613811 q21 = -0.00087760104664616457 q31 = 0.00091652546535745522 d21 = -0.036268708744722128 d31 = -0.00089075435142862297
- 1> I:\modular-boost\boost/math/tools/toms748_solve.hpp:275 q22 = -1.9862562616034592e-005 q32 = 3.1952597740788757e-007 d32 = -1.2833778805050512e-005 q33 = 1.1763429980834706e-008 c = 5.0000000047314881
- 1> I:\modular-boost\boost/math/tools/toms748_solve.hpp:388 a = 4.999083147976723 b = 5.0000000047314881
- 1> I:\modular-boost\boost/math/tools/toms748_solve.hpp:259 a = 4.999083147976723 b = 5.0000000047314881 d = 5.0008904277935091 e = 5.0362815354915851 fa = -2.8641119933622576 fb = 1.4785900475544622e-005 fd = 2.7835781082976609 fe = 115.03721886368339
- 1> I:\modular-boost\boost/math/tools/toms748_solve.hpp:267 q11 = -0.00087760104664616457 q21 = -4.7298032238887272e-009 q31 = 0.00091685202154135855 d21 = -0.00089042779182425238 d31 = -4.7332236912279757e-009
- 1> I:\modular-boost\boost/math/tools/toms748_solve.hpp:275 q22 = -1.6486403607318402e-012 q32 = 1.7346209428817704e-012 d32 = -1.6858463963666777e-012 q33 = 9.0382569995250912e-016 c = 5
- 1> I:\modular-boost\boost/math/tools/toms748_solve.hpp:592 max_iter = 10 count = 7
- 1> Iterations 20
- 1> 0 bits separate brackets.
- 1> fifth_noderiv(3125) = 5
- ]
- */
- //] [/root_finding_example_output]
|