Line data Source code
1 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 1 : 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) { 42 : ASSERT(digits < std::numeric_limits<double>::digits10, 43 : "The desired accuracy of " << digits 44 : << " base-10 digits must be smaller than " 45 : "the machine numeric limit of " 46 : << std::numeric_limits<double>::digits10 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) { 55 : throw convergence_error(MakeString{} 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 1 : 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) { 93 : ASSERT(digits < std::numeric_limits<double>::digits10, 94 : "The desired accuracy of " << digits 95 : << " base-10 digits must be smaller than " 96 : "the machine numeric limit of " 97 : << std::numeric_limits<double>::digits10 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) { return f(x, i); }, initial_guess[i], lower_bound[i], 107 : upper_bound[i], digits_binary, max_iters); 108 : if (max_iters >= max_iterations) { 109 : throw convergence_error(MakeString{} 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