diff options
Diffstat (limited to 'src/boost/libs/math/example/root_finding_fifth.cpp')
| -rw-r--r-- | src/boost/libs/math/example/root_finding_fifth.cpp | 485 |
1 files changed, 485 insertions, 0 deletions
diff --git a/src/boost/libs/math/example/root_finding_fifth.cpp b/src/boost/libs/math/example/root_finding_fifth.cpp new file mode 100644 index 00000000..a5f69938 --- /dev/null +++ b/src/boost/libs/math/example/root_finding_fifth.cpp @@ -0,0 +1,485 @@ +// 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] |