NewtonRaphson.hpp
Go to the documentation of this file.
2 // See LICENSE.txt for details.
3
4 /// \file
5 /// Declares function RootFinder::newton_raphson
6
7 #pragma once
8
9 #include <boost/math/tools/roots.hpp>
10 #include <functional>
11 #include <limits>
12
13 #include "DataStructures/DataVector.hpp"
14 #include "ErrorHandling/Exceptions.hpp"
15
16 namespace RootFinder {
17 /*!
18  * \ingroup NumericalAlgorithmsGroup
19  * \brief Finds the root of the function f with the Newton-Raphson method.
20  *
21  * f is a unary invokable that takes a double which is the current value at
22  * which to evaluate f. An example is below.
23  *
24  * \snippet Test_NewtonRaphson.cpp double_newton_raphson_root_find
25  *
26  * See the [Boost](http://www.boost.org/) documentation for more details.
27  *
28  * \requires Function f is invokable with a double
29  * \note The parameter digits specifies the precision of the result in its
30  * desired number of base-10 digits.
31  *
32  * \throws convergence_error if the requested precision is not met after
33  * max_iterations iterations.
34  */
35 template <typename Function>
36 double newton_raphson(const Function& f, const double initial_guess,
37  const double lower_bound, const double upper_bound,
38  const size_t digits, const size_t max_iterations = 50) {
40  "The desired accuracy of " << digits
41  << " base-10 digits must be smaller than "
42  "the machine numeric limit of "
44  << " base-10 digits.");
45
46  boost::uintmax_t max_iters = max_iterations;
47  // clang-tidy: internal boost warning, can't fix it.
48  const auto result = boost::math::tools::newton_raphson_iterate( // NOLINT
49  f, initial_guess, lower_bound, upper_bound,
50  std::round(std::log2(std::pow(10, digits))), max_iters);
51  if (max_iters >= max_iterations) {
52  throw convergence_error(
53  "newton_raphson reached max iterations without converging");
54  }
55  return result;
56 }
57
58 /*!
59  * \ingroup NumericalAlgorithmsGroup
60  * \brief Finds the root of the function f with the Newton-Raphson method on
61  * each element in a DataVector.
62  *
63  * f is a binary invokable that takes a double as its first argument and a
64  * size_t as its second. The double is the current value at which to
65  * evaluate f, and the size_t is the current index into the DataVectors.
66  * Below is an example of how to root find different functions by indexing into
67  * a lambda-captured DataVector using the size_t passed to f.
68  *
69  * \snippet Test_NewtonRaphson.cpp datavector_newton_raphson_root_find
70  *
71  * See the [Boost](http://www.boost.org/) documentation for more details.
72  *
73  * \requires Function f be callable with a double and a size_t
74  * \note The parameter digits specifies the precision of the result in its
75  * desired number of base-10 digits.
76  *
77  * \throws convergence_error if, for any index, the requested precision is not
78  * met after max_iterations iterations.
79  */
80 template <typename Function>
81 DataVector newton_raphson(const Function& f, const DataVector& initial_guess,
82  const DataVector& lower_bound,
83  const DataVector& upper_bound, const size_t digits,
84  const size_t max_iterations = 50) {
86  "The desired accuracy of " << digits
87  << " base-10 digits must be smaller than "
88  "the machine numeric limit of "
90  << " base-10 digits.");
91  const auto digits_binary = std::round(std::log2(std::pow(10, digits)));
92
93  DataVector result_vector{lower_bound.size()};
94  for (size_t i = 0; i < result_vector.size(); ++i) {
95  boost::uintmax_t max_iters = max_iterations;
96  // clang-tidy: internal boost warning, can't fix it.
97  result_vector[i] = boost::math::tools::newton_raphson_iterate( // NOLINT
98  [&f, i ](double x) noexcept { return f(x, i); }, initial_guess[i],
99  lower_bound[i], upper_bound[i], digits_binary, max_iters);
100  if (max_iters >= max_iterations) {
101  throw convergence_error(
102  "newton_raphson reached max iterations without converging");
103  }
104  }
105  return result_vector;
106 }
107
108 } // namespace RootFinder
Definition: GslMultiRoot.cpp:10
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
Exception indicating convergence failure.
Definition: Exceptions.hpp:10
Stores a collection of function values.
Definition: DataVector.hpp:46
DataVector newton_raphson(const Function &f, const DataVector &initial_guess, const DataVector &lower_bound, const DataVector &upper_bound, const size_t digits, const size_t max_iterations=50)
Finds the root of the function f with the Newton-Raphson method on each element in a DataVector...
Definition: NewtonRaphson.hpp:81
decltype(auto) constexpr pow(const T &t) noexcept
Compute t^N where N is an integer (positive or negative)
Definition: ConstantExpressions.hpp:157