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 "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 `DataVector`s.
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