math/doc/roots/root_finding_examples.qbk

[section:root_finding_examples Examples of Root-Finding (with and without derivatives)]

[import ../../example/root_finding_example.cpp]
[import ../../example/root_finding_n_example.cpp]
[import ../../example/root_finding_multiprecision_example.cpp]

The examples demonstrate how to use the various tools for
[@http://en.wikipedia.org/wiki/Root-finding_algorithm root finding].

We start with the simple cube root function `cbrt` ( C++ standard function name
[@http://en.cppreference.com/w/cpp/numeric/math/cbrt cbrt])
showing __cbrt_no_derivatives.

We then show how use of derivatives can improve the speed of convergence.

(But these examples are only a demonstration and do not try to make
the ultimate improvements of an 'industrial-strength'
implementation, for example, of `boost::math::cbrt`, mainly by using a better computed initial 'guess'
at [@boost:/libs/math/include/boost/math/special_functions/cbrt.hpp  cbrt.hpp]).

Then we show how a higher root (__fifth_root) [super 5][radic] can be computed,
and in
[@../../example/root_finding_n_example.cpp root_finding_n_example.cpp]
a generic method for the __nth_root that constructs the derivatives at compile-time.

These methods should be applicable to other functions that can be differentiated easily.

[h3:no_derivatives Finding a root without derivatives]

First some `#includes` that will be needed.

[root_finding_include_1]

[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).]

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.

Fortunately, the 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').

We then show how adding what we can know about this function, first just the slope
or 1st derivative ['f'(x)], will speed homing in on the solution.

Lastly, we show how adding the curvature ['f''(x)] too will speed convergence even more.

[h3:cbrt_no_derivatives Cube root function without derivatives]

First we define a function object (functor):

[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'.)

[root_finding_noderiv_2]

This snippet from `main()` in [@../../example/root_finding_example.cpp root_finding_example.cpp]
shows how it can be used.

[root_finding_main_1]

[pre
  cbrt_noderiv(27) = 3
  cbrt_noderiv(28) = 3.0365889718756618
]

The result of `bracket_and_solve_root` is a [@http://www.cplusplus.com/reference/utility/pair/ pair]
of values that could be displayed.

Tthe number of bits separating them can be found using `float_distance(r.first, r.second)`.
The distance is zero (closest representable) for 3[super 3] = 27
but `float_distance(r.first, r.second) = 3` for cube root of 28 with this function.
The result (avoiding overflow) is midway between these two values.

[h3:cbrt_1st_derivative Cube root function 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.

This algorithm is similar to this [@http://en.wikipedia.org/wiki/Nth_root_algorithm nth root algorithm].

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 define a better 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 [@http://en.cppreference.com/w/cpp/utility/pair std::pair]
of floating-point values.

[root_finding_1_deriv_1]

The result of [@boost:/libs/math/include/boost/math/tools/roots.hpp `newton_raphson_iterate`]
function is a single value.

[tip There is a compromise between accuracy and speed when chosing the value of `digits`.
The simplest approach is to use
 `std::numeric_limits<T>::digits`,  the maximum possible,
 but this may mean some inefficient iterations.
 A good compromise may be to use `3 * std::numeric_limits<T>::digits/4`.   See below.]

It is possible to tell if the call met the accuracy set in `get_digits`
by comparing the number of iterations, updated in `it`, with the maximum allowed in `max_it`.

Using the test data in [@../../test/test_cbrt.cpp  /test/test_cbrt.cpp] this found the cube root
exact to the last digit in every case, and in no more than 6 iterations at double
precision.  However, you will note that a high precision was used in this
example, exactly what was warned against earlier on in these docs!  In this
particular case it is possible to compute ['f(x)] exactly and without undue
cancellation error, so a high limit is not too much of an issue.

However, reducing the limit to `std::numeric_limits<T>::digits * 2 / 3` gave full
precision in all but one of the test cases (and that one was out by just one bit).
The maximum number of iterations remained 6, but in most cases was reduced by one.

Note also that the above code omits a probable optimization by computing z[sup2]
and reusing it, omits error handling, and does not handle
negative values of z correctly. (These are left as the customary exercise for the reader!)

The `boost::math::cbrt` function also includes these and other improvements.

[h3:cbrt_2_derivatives Cube root with 1st & 2nd derivative (slope & curvature)]

Next we define yet another even better functor `cbrt_functor_2deriv` that returns
both the evaluation of the function to solve,
along with its first [*and second] derivative:

__spaces['f''(x) = 6x]

using information about both slope and curvature to speed convergence.

To [''return'] three values, we use a `tuple` of three floating-point values:
[root_finding_2deriv_1]

The function `halley_iterate` also returns a single value,
and the number of iterations will reveal if it met the convergence criterion set by `get_digits`.

The no-derivative method gives a result of

  cbrt_noderiv(28) = 3.0365889718756618

with a 3 bits distance between the bracketed values, whereas the derivative methods both converge to a single value

  cbrt_2deriv(28) = 3.0365889718756627

which we can compare with the [@boost:/libs/math/doc/html/math_toolkit/powers/cbrt.html boost::math::cbrt]

  cbrt(28)              = 3.0365889718756627

Note that the iterations are set to stop at just one-half of full precision,
and yet, even so, not one of the test cases had a single bit wrong.
What's more, the maximum number of iterations was now just 4.

Just to complete the picture, we could have called
[link math_toolkit.roots.roots_deriv.schroeder `schroeder_iterate`] in the last
example: and in fact it makes no difference to the accuracy or number of iterations
in this particular case.  However, the relative performance of these two methods
may vary depending upon the nature of ['f(x)], and the accuracy to which the initial
guess can be computed.  There appear to be no generalisations that can be made
except "try them and see".

Finally, had we called `cbrt` with [@http://shoup.net/ntl/doc/RR.txt NTL::RR]
set to 1000 bit precision (about 300 decimal digits),
then full precision can be obtained with just 7 iterations.
To put that in perspective,
an increase in precision by a factor of 20, has less than doubled the number of
iterations.  That just goes to emphasise that most of the iterations are used
up getting the first few digits correct: after that these methods can churn out
further digits with remarkable efficiency.

Or to put it another way: ['nothing beats a really good initial guess!]

[h3:fifth_root 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 function whose complexity makes differentiation 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 __multiprecision_root).

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_functor_2deriv]
[root_finding_fifth_2deriv]

[h3:multiprecision_root  Root-finding using Boost.Multiprecision]

The apocryphally astute reader might, by now, be asking "How do we know if this computes the 'right' answer?".

For most values, there is, sadly, no 'right' answer.
This is because values can only rarely be ['exactly represented] by C++ floating-point types.
What we do want is the 'best' representation - one that is the nearest __representable value.
(For more about how numbers are represented see __floating_point).

Of course, we might start with finding an external reference source like
__WolframAlpha, as above, but this is not always possible.

Another way to reassure is to compute 'reference' values at higher precision
with which to compare the results of our iterative computations using built-in like `double`.
They should agree within the tolerance that was set.

The result  of `static_cast`ing to `double` from a higher-precision type like `cpp_bin_float_50` is guaranteed
to be the [*nearest representable] `double` value.

For example, the cube root functions in our example for cbrt(28.) return

`std::cbrt<double>(28.)  = 3.0365889718756627`

WolframAlpha says  `3.036588971875662519420809578505669635581453977248111123242141...`

`static_cast<double>(3.03658897187566251942080957850) = 3.0365889718756627`

This example `cbrt(28.)   =  3.0365889718756627`

[tip To ensure that all potentially significant decimal digits are displayed use `std::numeric_limits<T>::max_digits10`
(or if not available on older platforms or compilers use `2+std::numeric_limits<double>::digits*3010/10000`).[br]

Ideally, values should agree to `std::numeric-limits<T>::digits10` decimal digits.

This also means that a 'reference' value to be [*input] or `static_cast` should have
at least `max_digits10` decimal digits (17 for 64-bit `double`).
]

If we wish to compute [*higher-precision values] then, on some platforms, we may be able to use `long double`
with a higher precision than `double` to compare with the very common `double`
and/or a more efficient built-in quad floating-point type like `__float128`.

Almost all platforms can easily use __multiprecision,
for example, __cpp_dec_float or a binary type __cpp_bin_float types,
to compute values at very much higher precision.

[note With multiprecision types, it is debatable whether to use the type `T` for computing the initial guesses.
Type `double` is like to be accurate enough for the method used in these examples.
This would limit the range of possible values to that of `double`.
There is also the cost of conversion to and from type `T` to consider.
In these examples, `double` is used via `typedef double guess_type`.]

Since the functors and functions used above are templated on the value type,
we can very simply use them with any of the __multiprecision types.

Some examples below are 50 decimal digit decimal and binary types
(and on some platforms a much faster `float128` or `quad_float` type )
that we can use with these includes:

[root_finding_multiprecision_include_1]

Some using statements simplify their use:

[root_finding_multiprecision_example_1]

They can be used thus:

[root_finding_multiprecision_example_2]

A reference value computed by __WolframAlpha is

   N[2^(1/3), 50]  1.2599210498948731647672106072782283505702514647015

which agrees exactly.

To [*show] values to their full precision, it is necessary to adjust the `std::ostream` `precision` to suit the type, for example:

[root_finding_multiprecision_show_1]

[root_finding_multiprecision_example_3]

which outputs:

[pre
cbrt(2) = 1.2599210498948731647672106072782283505702514647015

value = 2, cube root =1.25992104989487
value = 2, cube root =1.25992104989487
value = 2, cube root =1.2599210498948731647672106072782283505702514647015
]

[tip Be [*very careful] about the floating-point type `T` that is passed to the root-finding function.
Carelessly passing a integer by writing
`cpp_dec_float_50  r = cbrt_2deriv(2);` or `show_cube_root(2);`
will provoke many warnings and compile errors.

Even `show_cube_root(2.F);` will produce warnings because `typedef double guess_type` defines the type
used to compute the guess and bracket values as `double`.

Even more treacherous is passing a `double` as in `cpp_dec_float_50  r = cbrt_2deriv(2.);`
which silently gives the 'wrong' result, computing a `double` result and [*then] converting to `cpp_dec_float_50`!
All digits beyond `max_digits10` will be incorrect.
Making the `cbrt` type explicit with `cbrt_2deriv<cpp_dec_float_50>(2.);` will give you the desired 50 decimal digit precision result.
] [/tip]

[h3:nth_root Generalizing to Compute the nth root]

If desired, we can now further generalize to compute the ['n]th root by computing the derivatives [*at compile-time]
using the rules for differentiation and the compile-time `static` function `pow<N>`
where template parameter `N` is an integer.  Our functor and function now have an additional template parameter `N`,
for the root required.

[note Since this is a compile-time static process, the resulting code is as efficient as as if hand-coded as the cube and fifth-root examples above.
A good compiler should also optimise any repeated multiplications.]

Our ['n]th root functor is

[root_finding_nth_functor_2deriv]

and our ['n]th root  function is

[root_finding_nth_function_2deriv]

[root_finding_n_example_2]

produces an output similar to this

[root_finding_example_output_1]

[tip Take care with the type passed to the function. It is best to pass a `double` or greater-precision floating-point type.

Passing an integer value, for example, `nth_2deriv<5>(2)` will be rejected, while `nth_2deriv<5, double>(2)` converts the integer to `double`.

Avoid passing a `float` value that will provoke warnings (actually spurious) from the compiler about potential loss of data,
as noted above.]

[warning Asking for unreasonable roots, for example, `show_nth_root<1000000>(2.);` may lead to
[@http://en.wikipedia.org/wiki/Loss_of_significance Loss of significance] like
`Type double value = 2, 1000000th root = 1.00000069314783`.
Use of the the `pow` function is more sensible for this unusual need.
]

Full code of these examples is at
[@../../example/root_finding_example.cpp  root_finding_example.cpp],
[@../../example/root_finding_n_example.cpp  root_finding_n_example.cpp] and
[@../../example/root_finding_multiprecision_example.cpp  root_finding_multiprecision_example.cpp].

[endsect] [/section:root_examples Examples of Root Finding (with and without derivatives)]

[/
  Copyright 2015 John Maddock and Paul A. Bristow.
  Distributed under 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).
]