corecel/math/RegulaFalsiRootFinder.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2024 UT-Battelle, LLC, and other Celeritas developers.
0003 // See the top-level COPYRIGHT file for details.
0004 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0005 //---------------------------------------------------------------------------//
0006 //! \file corecel/math/RegulaFalsiRootFinder.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include <cmath>
0011 #include <type_traits>
0012
0013 #include "corecel/Config.hh"
0014
0015 #include "corecel/Assert.hh"
0016 #include "corecel/Macros.hh"
0017 #include "corecel/Types.hh"
0018 #include "corecel/math/Algorithms.hh"
0019
0020 namespace celeritas
0021 {
0022 //---------------------------------------------------------------------------//
0023 /*!
0024  * Perform Regula Falsi iterations given a root function \em func and
0025  * tolerance \em tol .
0026  *
0027  * Using a \em left and \em right bound a Regula Falsi approximates the \em
0028  * root as: \f[ root = (left * func(right) - right * func(left)) / (func(right)
0029  * - func(left)) \f]
0030  *
0031  * Then value of \em func at the root is calculated compared to values of
0032  * \em func at the bounds. The \em root is then used update the bounds based on
0033  * the sign of \em func(root) and whether it matches the sign of \em func(left)
0034  * or \em func(right) . Performing this update of the bounds allows for the
0035  * iteration on root, using the convergence criteria based on \em func(root)
0036  * proximity to 0.
0037  */
0038 template<class F>
0039 class RegulaFalsiRootFinder
0040 {
0041   public:
0042     // Contruct with function to solve and solution tolerance
0043     inline CELER_FUNCTION RegulaFalsiRootFinder(F&& func, real_type tol);
0044
0045     // Solve for a root between two points
0046     inline CELER_FUNCTION real_type operator()(real_type left, real_type right);
0047
0048   private:
0049     F func_;
0050     real_type tol_;
0051
0052     // Maximum amount of iterations
0053     static constexpr inline int max_iters_ = 50;
0054 };
0055
0056 //---------------------------------------------------------------------------//
0057 // TEMPLATE DEDUCTION
0058 //---------------------------------------------------------------------------//
0059
0060 template<class F, class... Args>
0061 CELER_FUNCTION RegulaFalsiRootFinder(F&&, Args...) -> RegulaFalsiRootFinder<F>;
0062
0063 //---------------------------------------------------------------------------//
0064 // INLINE DEFINITIONS
0065 //---------------------------------------------------------------------------//
0066 /*!
0067  * Construct from function.
0068  */
0069 template<class F>
0070 CELER_FUNCTION
0071 RegulaFalsiRootFinder<F>::RegulaFalsiRootFinder(F&& func, real_type tol)
0072     : func_{celeritas::forward<F>(func)}, tol_{tol}
0073 {
0074     CELER_EXPECT(tol_ > 0);
0075 }
0076
0077 //---------------------------------------------------------------------------//
0078 /*!
0079  * Solve for a root between the two points.
0080  */
0081 template<class F>
0082 CELER_FUNCTION real_type RegulaFalsiRootFinder<F>::operator()(real_type left,
0083                                                               real_type right)
0084 {
0085     // Initialize Iteration parameters
0086     real_type f_left = func_(left);
0087     real_type f_right = func_(right);
0088     real_type f_root = 1;
0089     real_type root = 0;
0090     int remaining_iters = max_iters_;
0091
0092     // Iterate on root
0093     do
0094     {
0095         // Estimate root and update value
0096         root = (left * f_right - right * f_left) / (f_right - f_left);
0097         f_root = func_(root);
0098
0099         // Update the bound which produces the same sign as the root
0100         if (signum(f_left) == signum(f_root))
0101         {
0102             left = root;
0103             f_left = f_root;
0104         }
0105         else
0106         {
0107             right = root;
0108             f_right = f_root;
0109         }
0110     } while (std::fabs(f_root) > tol_ && --remaining_iters > 0);
0111
0112     CELER_ENSURE(remaining_iters > 0);
0113     return root;
0114 }
0115
0116 //---------------------------------------------------------------------------//
0117 }  // namespace celeritas