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