boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-30 20:12:09 +00:00
Code Issues Projects Releases Wiki Activity

rooting links corrected, build clean, but some links still broken.

Browse Source
This commit is contained in:
pabristow
2015-03-12 18:05:36 +00:00
parent 4257aa91ff
commit 2b57792dbc
12 changed files with 539 additions and 390 deletions
Download Patch File Download Diff File Expand all files Collapse all files

614
example/root_finding_example.cpp
View File

@@ -1,369 +1,490 @@
// root_finding_example.cpp
// Copyright Paul A. Bristow 2010, 2014.
// 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, Schroeder, TOMS748 .
// 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!
// To get (copious) diagnostic output, add make this define here or elsewhere.
//#ifdef _MSC_VER
//# pragma warning(disable: 4180) // qualifier has no effect (in Fusion).
//#endif
//#define BOOST_MATH_INSTRUMENT
//[root_finding_headers
/*
This example demonstrates how to use the various tools for root finding
taking the simple cube root function (cbrt) as an example.
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>` ).
implementation of `boost::math::cbrt`, mainly by using a better computed initial 'guess'
at `<boost/math/special_functions/cbrt.hpp>`).
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.
Then we show how a high roots (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;
//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 <tuple>
#include <utility> // pair, make_pair
#include <boost/math/special_functions/next.hpp> // For float_distance.
#include <boost/math/tools/tuple.hpp> // for tuple and make_tuple.
#include <boost/math/special_functions/cbrt.hpp> // For boost::math::cbrt.
//] [/root_finding_headers]
//] [/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;
//using std::cout; using std::endl;
#include <iomanip>
using std::setw; using std::setprecision;
//using std::setw; using std::setprecision;
#include <limits>
using std::numeric_limits;
//using std::numeric_limits;
/*
Let's suppose we want to find the cube root of a number.
Let's suppose we want to find the root of a number ['a], and to start, compute the cube root.
The equation we want to solve is:
So the equation we want to solve is:
__spaces ['f](x) = x[cubed]
__spaces ['f](x) = x[cubed] - a
We will first solve this without using any information
about the slope or curvature of the cbrt function.
about the slope or curvature of the cube root function.
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.
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.
The 1st and 2nd derivatives of x[cubed] are:
Lastly we show how adding the curvature /f''(x)/ too will speed convergence even more.
__spaces ['f]\'(x) = 2x[sup2]
__spaces ['f]\'\'(x) = 6x
*/
//[root_finding_cbrt_functor_noderiv
//[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) : value(to_find_root_of)
{ // Constructor stores value to find root of.
// For example: calling cbrt_functor<T>(x) to get cube root of x.
{ // 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)
{ //! \returns f(x) - value.
T fx = x*x*x - value; // Difference (estimate x^3 - value).
{
T fx = x*x*x - a; // Difference (estimate x^3 - a).
return fx;
}
private:
T value; // to be 'cube_rooted'.
T a; // to be 'cube_rooted'.
};
//] [/root_finding_noderiv_1
//] [/root_finding_cbrt_functor_noderiv]
//cout << ", std::numeric_limits<" << typeid(T).name() << ">::digits = " << digits
// << ", accuracy " << get_digits << " bits."<< endl;
/*`Implementing the cube root function itself is fairly trivial now:
/*
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
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_cbrt_noderiv
//[root_finding_noderiv_2
template <class T>
T cbrt_noderiv(T x)
{ //! \returns cube root of x using bracket_and_solve (no derivatives).
{ // return cube root of x using bracket_and_solve (no derivatives).
using namespace std; // Help ADL of std functions.
using namespace boost::math; // For bracket_and_solve_root.
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; // 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.
T factor = 2; // How big steps to take when searching.
//boost::uintmax_t maxit = (std::numeric_limits<boost::uintmax_t>::max)();
// (std::numeric_limits<boost::uintmax_t>::max)() = 18446744073709551615
// which is more than we might wish to wait for!!!
// so we can choose some reasonable estimate of how many iterations may be needed.
const boost::uintmax_t maxit = 10;
const boost::uintmax_t maxit = 20; // Limit to maximum iterations.
boost::uintmax_t it = maxit; // Initally our chosen max iterations, but updated with actual.
// We could also have used a maximum iterations provided by any policy:
// boost::uintmax_t max_it = policies::get_max_root_iterations<Policy>();
bool is_rising = true; // So if result if guess^3 is too low, try increasing guess.
eps_tolerance<double> tol(digits);
std::pair<T, T> r =
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) /4; // Near maximum (3/4) possible accuracy.
eps_tolerance<double> 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);
// Can show how many iterations (this information is lost outside cbrt_noderiv).
cout << "Iterations " << maxit << endl;
if(it >= maxit)
{ //
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);
std::cout << distance << " bits separate brackets." << std::endl;
for (int i = 0; i < distance; i++)
{
std::cout << float_advance(r.first, i) << std::endl;
}
return r.first + (r.second - r.first) / 2; // Midway between bracketed interval.
} // T cbrt_noderiv(T x)
//] [/root_finding_cbrt_noderiv]
// maxit = 10
// Unable to locate solution in chosen iterations: Current best guess is between 3.0365889718756613 and 3.0365889718756627
return r.first + (r.second - r.first)/2; // Midway between brackets.
}
/*`
[note Default `boost::uintmax_t maxit = (std::numeric_limits<boost::uintmax_t>::max)();`
and `(std::numeric_limits<boost::uintmax_t>::max)() = 18446744073709551615`
which 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, here only 20.]
[note We could also have used a maximum iterations provided by any policy:
`boost::uintmax_t max_it = policies::get_max_root_iterations<Policy>();`]
[tip One can show how many iterations 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 this function, the 1st differential (the slope of the tangent to a curve at any point) is known.
For the root 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.
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^n for positive n (actually all nonzero n) is nx^n-1,
allows use to get the 1st differential as 3x^2.
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 the root, view this demo:
[@http://en.wikipedia.org/wiki/Newton%27s_methodNewton Newton-Raphson iterations].
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 that returns
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 pair of floating-point values:
To \'return\' two values, we use a `std::pair` of floating-point values:
*/
//[root_finding_cbrt_functor_1stderiv
//[root_finding_1_deriv_1
template <class T>
struct cbrt_functor_1stderiv
{ // Functor returning function and 1st derivative.
cbrt_functor_1stderiv(T const& target) : value(target)
{ // Constructor stores the value to be 'cube_rooted'.
struct cbrt_functor_deriv
{ // Functor also returning 1st derviative.
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>(x) to use to get cube root of x.
}
std::pair<T, T> operator()(T const& z) // z is best estimate so far.
{ // Return both f(x) and derivative f'(x).
T fx = z*z*z - value; // Difference estimate fx = x^3 - value.
T d1x = 3 * z*z; // 1st derivative d1x = 3x^2.
return std::make_pair(fx, d1x); // 'return' both fx and d1x.
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 value; // to be 'cube_rooted'.
}; // cbrt_functor_1stderiv
T a; // to be 'cube_rooted'.
};
//] [/root_finding_cbrt_functor_1stderiv]
/*`Our cube root function using cbrt_functor_1stderiv is now:*/
//[root_finding_cbrt_1deriv
/*`Our cube root function is now:*/
template <class T>
T cbrt_1deriv(T x)
{ //! \return cube 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/3); // Rough guess is to divide the exponent by three.
// Set an initial bracket interval.
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.
// 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.
boost::uintmax_t maxit = 20; // Optionally limit the number of iterations.
//cout << "Max Iterations " << maxit << endl; //
T result = newton_raphson_iterate(cbrt_functor_1stderiv<T>(x), guess, min, max, digits, maxit);
// Can check and show how many iterations (updated by newton_raphson_iterate).
// cout << "Iterations " << maxit << endl;
return result;
} // cbrt_1deriv
//] [/root_finding_cbrt_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_cbrt_functor_2deriv
template <class T>
struct cbrt_functor_2deriv
{ // Functor returning both 1st and 2nd derivatives.
cbrt_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 - value; // 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 value; // to be 'cube_rooted'.
}; // struct cbrt_functor_2deriv
//] [/root_finding_cbrt_functor_2deriv]
/*`Our cube function is now:*/
//[root_finding_cbrt_2deriv
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; // halley_iterate
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.
int digits = std::numeric_limits<T>::digits /2 ; // Half maximum possible binary digits accuracy for type T.
boost::uintmax_t maxit = 10;
T result = halley_iterate(cbrt_functor_2deriv<T>(x), guess, min, max, digits, maxit);
// Can show how many iterations (updated by halley_iterate).
// cout << "Iterations " << maxit << endl;
const int digits = std::numeric_limits<T>::digits; // Maximum possible binary digits accuracy for type T.
int get_digits = (digits * 3) /4; // Near maximum (3/4) possible accuracy.
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;
} // cbrt_2deriv(x)
}
//] [/root_finding_cbrt_2deriv]
//] [/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,
}
boost::math::tuple<T, T, T> operator()(T const& x)
{ // Return both f(x) and f'(x) and f''(x).
using boost::math::make_tuple;
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 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 = digits/2; // Near maximum (1/2) possible accuracy.
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]
/*
[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 = x*x*x*x*x - a; // 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 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.
const 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);
return result;
}
//] [/root_finding_fifth_2deriv]
// Default is
//boost::uintmax_t maxit = (std::numeric_limits<boost::uintmax_t>::max)();
// the default is (std::numeric_limits<boost::uintmax_t>::max)() = 18446744073709551615
// which is more than we might wish to wait for!!! so we can reduce it.
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_example_1
cout << "Cube Root finding (cbrt) 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);
//[root_finding_main_1
try
{ // Always use try'n'catch blocks with Boost.Math to get any error messages.
{
double threecubed = 27.; // value that has an exactly representable integer cube root.
double threecubedp1 = 28.; // whose cube root is *not* exactly representable.
double v27 = 27; // Example of a value that has an exact integer cube root.
double v28 = 28; // Example of a 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;
// Using bracketing:
double r = cbrt_noderiv(v27);
cout << "cbrt_noderiv(" << v27 << ") = " << r << endl;
r = cbrt_noderiv(v28);
cout << "cbrt_noderiv(" << v28 << ") = " << r << 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
// Using 1st differential Newton-Raphson:
r = cbrt_1deriv(v27);
cout << "cbrt_1deriv(" << v27 << ") = " << r << endl;
r = cbrt_1deriv(v28);
cout << "cbrt_1deriv(" << v28 << ") = " << r << endl;
// 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;
// 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 = cbrt_2deriv(v27);
cout << "cbrt_2deriv(" << v27 << ") = " << r << endl;
r = cbrt_2deriv(v28);
cout << "cbrt_2deriv(" << v28 << ") = " << r << endl;
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.
{ // 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:
Normal output is:
[pre
Description: Autorun "J:\Cpp\MathToolkit\test\Math_test\Release\root_finding_example.exe"1> Cube Root finding (cbrt) Example.
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_noderiv(27) = 3
cbrt_1(27) = 3
Iterations 10
Unable to locate solution in chosen iterations: Current best guess is between 3.0365889718756613 and 3.0365889718756627
cbrt_noderiv(28) = 3.0365889718756618
cbrt_1deriv(27) = 3
cbrt_1deriv(28) = 3.0365889718756627
cbrt_2deriv(27) = 3
cbrt_2deriv(28) = 3.0365889718756627
]
[/pre]
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
@@ -374,3 +495,10 @@ to get some (much!) diagnostic output we can add
]
*/
//] [/root_finding_example_output]
/*
cbrt(28) 3.0365889718756622
std::cbrt(28) 3.0365889718756627
*/