486 lines
19 KiB
C++
486 lines
19 KiB
C++
// 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]
|