NewtonRaphson.hpp
Go to the documentation of this file.
1 // Distributed under the MIT License.
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 `DataVector`s.
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:42
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