include/Geant4/G4FunctionSolver.hh

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // G4FunctionSolver
0027 //
0028 // Class description:
0029 //
0030 // Templated utility class to solve equation F(x) = 0 on interval [a, b]
0031 // Parameters of the class: tolerance - relative accuracy of the method,
0032 // maxIter - max number of iterations. User should provide a class with
0033 // only one method: T_Function::Function(G4double x). The main method
0034 // of the class is 'G4bool G4FunctionSolver::FindRoot(G4double& x)', where
0035 // user defines initial value x, which is modifided to be a final solution
0036 // of the equation.
0037 //
0038 // If the equation cannot be resolved with required accuracy FindRoot()
0039 // method returns "false".
0040 //
0041 // Parameters of the solver may be changed before the new call via Set
0042 // methods.
0043 //
0044 // Created 04.10.2025 V.Ivanchenko
0045 //
0046
0047 #ifndef G4FunctionSolver_h
0048 #define G4FunctionSolver_h 1
0049
0050 #include "globals.hh"
0051
0052 template <class T_Function> class G4FunctionSolver
0053 {
0054 public:
0055
0056   G4FunctionSolver(T_Function* ff, const G4int iterations, const G4double tol)
0057     : maxIter(iterations), tolerance(tol), tF(ff) {};
0058
0059   // copy constructor
0060   G4FunctionSolver(const G4FunctionSolver& right) = delete;
0061
0062   // destructor
0063   ~G4FunctionSolver() = default;
0064
0065   // operators
0066   G4FunctionSolver& operator=(const G4FunctionSolver& right) = delete;
0067   G4bool operator==(const G4FunctionSolver& right) const = delete;
0068   G4bool operator!=(const G4FunctionSolver& right) const = delete;
0069
0070   inline void SetMaxIterations(const G4int iterations)
0071   { maxIter = iterations; }
0072
0073   inline void SetTolerance(const G4double epsilon)
0074   { tolerance = epsilon; }
0075
0076   inline void SetIntervalLimits(const G4double Limit1, const G4double Limit2)
0077   {
0078     aa = std::min(Limit1, Limit2);
0079     bb = std::max(Limit1, Limit2);
0080   }
0081
0082   // Calculates the root of the equation Function(x)=0
0083   inline G4bool FindRoot(G4double& x)
0084   {
0085     G4double a = aa;
0086     G4double b = bb;
0087
0088     // check the interval before the start
0089     x = std::min(std::max(x, a), b);
0090
0091     // check initial function
0092     G4double fc = tF->Function(x);
0093     if (0.0 == fc) { return true; }
0094
0095     // define accuracy in X
0096     G4double epsX = tolerance*x;
0097     // define accuracy in Y
0098     G4double epsY = tolerance*fc;
0099
0100     // the interval is too small
0101     if (std::abs(a - b) <= epsX)
0102     {
0103       x = 0.5*(a + b);
0104       return true;
0105     }
0106
0107     // check edges
0108     G4double fa = tF->Function(a);
0109     G4double fb = tF->Function(b);
0110
0111     // root should be inside interval
0112     if (fa*fb >= 0.0)
0113     {
0114       x = (std::abs(fa) <= std::abs(fb)) ? a : b;
0115       return (std::min(std::abs(fa), std::abs(fb)) < epsY);
0116     }
0117
0118     // fa*fb < 0.0 - finding the root by iterative procedure,
0119     // the loop is completed if function is below epsY
0120     // or if x become close to edges with accuracy epsX
0121     for (G4int i = 0; i < maxIter; ++i)
0122     {
0123       x = (a*fb - b*fa)/(fb - fa);
0124       fc = tF->Function(x);
0125       if (std::abs(fc) < epsY) { return true; }
0126
0127       G4double delta = std::min((x - a), (b - x));
0128       if (delta < epsX) { return true; }
0129       else if (fa*fc < 0.0) { b = x; fb = fc; }
0130       else { a = x; fa = fc; }
0131     }
0132     // number of iterations exceed the limit
0133     return false;
0134   }
0135
0136 private:
0137
0138   // maximum number of iterations
0139   G4int maxIter;
0140   // relative accuracy in X and Y
0141   G4double tolerance;
0142
0143   // interval limits [a,b]
0144   G4double aa{0.0};
0145   G4double bb{0.0};
0146
0147   T_Function* tF;
0148 };
0149
0150 #endif