123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598 |
- [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 root finding __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.
- [section:cbrt_eg Finding the Cubed Root With and 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:
- [expression ['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.
- The 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`.
- It is tempting to simply chose `std::numeric_limits<T>::digits`,
- but this may mean some inefficient and unnecessary iterations as the function thrashes around
- trying to locate the last bit. In theory, since the precision doubles with each step
- it is sufficient to stop when half the bits are correct: as the last step will have doubled
- that to full precision. Of course the function has no way to tell if that is actually the case
- unless it does one more step to be sure. In practice setting the precision to slightly more
- than `std::numeric_limits<T>::digits / 2` is a good choice.]
- Note that it is up to the caller of the function to check the iteration count
- after the call to see if iteration stoped as a result of running out of iterations
- rather than meeting the required precision.
- 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:
- most importantly it uses a much better initial guess which reduces the iteration count to
- just 1 in almost all cases.
- [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:
- [expression ['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_deriv.schroder `schroder_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!]
- Full code of this example is at
- [@../../example/root_finding_example.cpp root_finding_example.cpp],
- [endsect] [/section:cbrt_eg Finding the Cubed Root With and Without Derivatives]
- [section:lambda Using C++11 Lambda's]
- Since all the root finding functions accept a function-object, they can be made to
- work (often in a lot less code) with C++11 lambda's. Here's the much reduced code for our "toy" cube root function:
- [root_finding_2deriv_lambda]
- Full code of this example is at
- [@../../example/root_finding_example.cpp root_finding_example.cpp],
- [endsect] [/section:lambda Using C++11 Lambda's]
- [section:5th_root_eg Computing the Fifth Root]
- Let's now suppose we want to find the [*fifth root] of a number ['a].
- The equation we want to solve is :
- [expression ['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:
- [expression ['f\'(x) = 5x[super 4]]]
- [expression ['f\'\'(x) = 20x[super 3]]]
- [root_finding_fifth_functor_2deriv]
- [root_finding_fifth_2deriv]
- Full code of this example is at
- [@../../example/root_finding_example.cpp root_finding_example.cpp] and
- [@../../example/root_finding_n_example.cpp root_finding_n_example.cpp].
- [endsect] [/section:5th_root_eg Computing the Fifth Root]
- [section: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.)` compute
- `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 exponent 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. As a reminder,
- here's our toy cube root function using 2 derivatives and C++11 lambda functions to find the root:
- [root_finding_2deriv_lambda]
- 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]
- Full code of this example is at
- [@../../example/root_finding_multiprecision_example.cpp root_finding_multiprecision_example.cpp].
- [endsect] [/section:multiprecision_root Root-finding using Boost.Multiprecision]
- [section: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 `boost::math::pow<N>`
- where template parameter `N` is an integer and a compile time constant. Our functor and function now have an additional template parameter `N`,
- for the root required.
- [note Since the powers and derivatives are fixed at compile time, 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 this example is at
- [@../../example/root_finding_n_example.cpp root_finding_n_example.cpp].
- [endsect]
- [section:elliptic_eg A More complex example - Inverting the Elliptic Integrals]
- The arc length of an ellipse with radii ['a] and ['b] is given by:
- [pre L(a, b) = 4aE(k)]
- with:
- [pre k = [sqrt](1 - b[super 2]/a[super 2])]
- where ['E(k)] is the complete elliptic integral of the second kind - see __ellint_2.
- Let's suppose we know the arc length and one radii, we can then calculate the other
- radius by inverting the formula above. We'll begin by encoding the above formula
- into a functor that our root-finding algorithms can call.
- Note that while not
- completely obvious from the formula above, the function is completely symmetrical
- in the two radii - which can be interchanged at will - in this case we need to
- make sure that `a >= b` so that we don't accidentally take the square root of a negative number:
- [import ../../example/root_elliptic_finding.cpp]
- [elliptic_noderv_func]
- We'll also need a decent estimate to start searching from, the approximation:
- [pre L(a, b) [approx] 4[sqrt](a[super 2] + b[super 2])]
- Is easily inverted to give us what we need, which using derivative-free root
- finding leads to the algorithm:
- [elliptic_root_noderiv]
- This function generally finds the root within 8-10 iterations, so given that the runtime
- is completely dominated by the cost of calling the ellliptic integral it would be nice to
- reduce that count somewhat. We'll try to do that by using a derivative-based method;
- the derivatives of this function are rather hard to work out by hand, but fortunately
- [@http://www.wolframalpha.com/input/?i=d%2Fda+\[4+*+a+*+EllipticE%281+-+b^2%2Fa^2%29\]
- Wolfram Alpha] can do the grunt work for us to give:
- [pre d/da L(a, b) = 4(a[super 2]E(k) - b[super 2]K(k)) / (a[super 2] - b[super 2])]
- Note that now we have [*two] elliptic integral calls to get the derivative, so our
- functor will be at least twice as expensive to call as the derivative-free one above:
- we'll have to reduce the iteration count quite substantially to make a difference!
- Here's the revised functor:
- [elliptic_1deriv_func]
- The root-finding code is now almost the same as before, but we'll make use of
- Newton-iteration to get the result:
- [elliptic_1deriv]
- The number of iterations required for `double` precision is now usually around 4 -
- so we've slightly more than halved the number of iterations, but made the
- functor twice as expensive to call!
- Interestingly though, the second derivative requires no more expensive
- elliptic integral calls than the first does, in other words it comes
- essentially "for free", in which case we might as well make use of it
- and use Halley-iteration. This is quite a typical situation when
- inverting special-functions. Here's the revised functor:
- [elliptic_2deriv_func]
- The actual root-finding code is almost the same as before, except we can
- use Halley, rather than Newton iteration:
- [elliptic_2deriv]
- While this function uses only slightly fewer iterations (typically around 3)
- to find the root, compared to the original derivative-free method, we've moved from
- 8-10 elliptic integral calls to 6.
- Full code of this example is at
- [@../../example/root_elliptic_finding.cpp root_elliptic_finding.cpp].
- [endsect]
- [endsect] [/section:root_examples Examples of Root Finding (with and without derivatives)]
- [section:bad_guess The Effect of a Poor Initial Guess]
- It's instructive to take our "toy" example algorithms, and use deliberately bad initial guesses to see how the
- various root finding algorithms fair. We'll start with the cubed root, and using the cube root of 500 as the test case:
- [table
- [[Initial Guess=][-500% ([approx]1.323)][-100% ([approx]3.97)][-50% ([approx]3.96)][-20% ([approx]6.35)][-10% ([approx]7.14)][-5% ([approx]7.54)][5% ([approx]8.33)][10% ([approx]8.73)][20% ([approx]9.52)][50% ([approx]11.91)][100% ([approx]15.87)][500 ([approx]47.6)]]
- [[bracket_and_solve_root][12][8][8][10][11][11][11][11][11][11][7][13]]
- [[newton_iterate][12][7][7][5][5][4][4][5][5][6][7][9]]
- [[halley_iterate][7][4][4][3][3][3][3][3][3][4][4][6]]
- [[schroder_iterate][11][6][6][4][3][3][3][3][4][5][5][8]]
- ]
- As you can see `bracket_and_solve_root` is relatively insensitive to starting location - as long as you don't start many orders of magnitude away from the root it will
- take roughly the same number of steps to bracket the root and solve it. On the other hand the derivative-based methods are slow to start, but once they have some digits
- correct they increase precision exceptionally fast: they are therefore quite sensitive to the initial starting location.
- The next table shows the number of iterations required to find the second radius of an ellipse with first radius 50 and arc-length 500:
- [table
- [[Initial Guess=][-500% ([approx]20.6)][-100% ([approx]61.81)][-50% ([approx]61.81)][-20% ([approx]98.9)][-10% ([approx]111.3)][-5% ([approx]117.4)][5% ([approx]129.8)][10% ([approx]136)][20% ([approx]148.3)][50% ([approx]185.4)][100% ([approx]247.2)][500 ([approx]741.7)]]
- [[bracket_and_solve_root][11][5][5][8][8][7][7][8][9][8][6][10]]
- [[newton_iterate][4][4][4][3][3][3][3][3][3][4][4][4]]
- [[halley_iterate][4][3][3][3][3][2][2][3][3][3][3][3]]
- [[schroder_iterate][4][3][3][3][3][2][2][3][3][3][3][3]]
- ]
- Interestingly this function is much more resistant to a poor initial guess when using derivatives.
- [endsect]
- [section:bad_roots Examples Where Root Finding Goes Wrong]
- There are many reasons why root root finding can fail, here are just a few of the more common examples:
- [h3 Local Minima]
- If you start in the wrong place, such as z[sub 0] here:
- [$../roots/bad_root_1.svg]
- Then almost any root-finding algorithm will descend into a local minima rather than find the root.
- [h3 Flatlining]
- In this example, we're starting from a location (z[sub 0]) where the first derivative is essentially zero:
- [$../roots/bad_root_2.svg]
- In this situation the next iteration will shoot off to infinity (assuming we're using derivatives that is). Our
- code guards against this by insisting that the root is always bracketed, and then never stepping outside those bounds.
- In a case like this, no root finding algorithm can do better than bisecting until the root is found.
- Note that there is no scale on the graph, we have seen examples of this situation occur in practice ['even when
- several decimal places of the initial guess z[sub 0] are correct.]
- This is really a special case of a more common situation where root finding with derivatives is ['divergent]. Consider
- starting at z[sub 0] in this case:
- [$../roots/bad_root_4.svg]
- An initial Newton step would take you further from the root than you started, as will all subsequent steps.
- [h3 Micro-stepping / Non-convergence]
- Consider starting at z[sub 0] in this situation:
- [$../roots/bad_root_3.svg]
- The first derivative is essentially infinite, and the second close to zero (and so offers no correction if we use it),
- as a result we take a very small first step. In the worst case situation, the first step is so small
- - perhaps even so small that subtracting from z[sub 0] has no effect at the current working precision - that our algorithm
- will assume we are at the root already and terminate. Otherwise we will take lot's of very small steps which never converge
- on the root: our algorithms will protect against that by reverting to bisection.
- An example of this situation would be trying to find the root of e[super -1/z[super 2]] - this function has a single
- root at ['z = 0], but for ['z[sub 0] < 0] neither Newton nor Halley steps will ever converge on the root, and for ['z[sub 0] > 0]
- the steps are actually divergent.
- [endsect]
- [/
- 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).
- ]
|