rootsolver.h

Go to the documentation of this file.

/*
 *
 *    Copyright (c) 2023
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 *
 */
 
#ifndef SRC_INCLUDE_SMASH_ROOTSOLVER_H_
#define SRC_INCLUDE_SMASH_ROOTSOLVER_H_
 
#include <functional>
#include <memory>
#include <optional>
#include <string>
 
#include "gsl/gsl_errno.h"
#include "gsl/gsl_math.h"
#include "gsl/gsl_roots.h"
 
#include "logging.h"
 
namespace smash {
static constexpr int LRootSolver = LogArea::RootSolver::id;
class RootSolver1D {
 public:
  explicit RootSolver1D(std::function<double(double)> eq) {
    root_eq_ = std::make_unique<std::function<double(double)>>(eq);
  }
 
  ~RootSolver1D() { root_eq_ = nullptr; }
 
  std::optional<double> try_find_root(double initial_guess_low,
                                      double initial_guess_high,
                                      size_t itermax) {
    // check if root is in the given interval
    if ((*root_eq_)(initial_guess_low) * (*root_eq_)(initial_guess_high) > 0) {
      logg[LRootSolver].debug()
          << "Function has same sign at both ends of the interval ["
          << initial_guess_low << ", " << initial_guess_high
          << "]. Root can't be found in this interval.";
      return std::nullopt;
    }
    gsl_function function_GSL = {&(gsl_func), nullptr};
    int status = GSL_CONTINUE;
    size_t iter = 0;
    Root_finder_ = gsl_root_fsolver_alloc(Solver_name_);
    gsl_root_fsolver_set(Root_finder_, &function_GSL, initial_guess_low,
                         initial_guess_high);
    do {
      iter++;
      status = gsl_root_fsolver_iterate(Root_finder_);
      if (status != GSL_SUCCESS) {
        logg[LRootSolver].debug("GSL ERROR in root finding: " +
                                static_cast<std::string>(gsl_strerror(status)));
        break;
      }
      double xlow = gsl_root_fsolver_x_lower(Root_finder_);
      double xhigh = gsl_root_fsolver_x_upper(Root_finder_);
      status = gsl_root_test_interval(xlow, xhigh, 0, solution_precision_);
      if (status == GSL_SUCCESS) {
        double root = 0.5 * (xlow + xhigh);
        gsl_root_fsolver_free(Root_finder_);
        return root;
      }
    } while (status == GSL_CONTINUE && iter < itermax);
    gsl_root_fsolver_free(Root_finder_);
    return std::nullopt;
  }
 
 private:
  const gsl_root_fsolver_type *Solver_name_ = gsl_root_fsolver_brent;
 
  gsl_root_fsolver *Root_finder_ = nullptr;
 
  static inline std::unique_ptr<std::function<double(double)>> root_eq_ =
      nullptr;
 
  double solution_precision_ = 1e-7;
 
  static double gsl_func(const double x, void *) { return (*root_eq_)(x); }
};
 
}  // namespace smash
 
#endif  // SRC_INCLUDE_SMASH_ROOTSOLVER_H_