TUHKMgen/UHKM/EquationSolver.icc

   1 /*                                                                            \r
   2                                                                             \r
   3         Nikolai Amelin, Ludmila Malinina, Timur Pocheptsov (C) JINR/Dubna\r
   4       amelin@sunhe.jinr.ru, malinina@sunhe.jinr.ru, pocheptsov@sunhe.jinr.ru \r
   5                            November. 2, 2006                                \r
   6 \r
   7 */ \r
   8 //This equation solver class is taken from GEANT4 and modified!!\r
   9  \r
  10 \r
  11 template <class Function>\r
  12 Bool_t NAEquationSolver<Function>::Brent(Function& theFunction) {\r
  13   const Double_t precision = 3.0e-8;\r
  14 \r
  15   // Check the interval before start\r
  16   if(fA > fB || Abs(fA-fB) <= fTolerance) {\r
  17     std::cerr << "NAEquationSolver::Brent: The interval must be properly set." << std::endl;\r
  18     return false;\r
  19   }\r
  20   Double_t fa = theFunction(fA);\r
  21   Double_t fb = theFunction(fB);\r
  22   if(fa*fb > 0.0) {\r
  23     std::cerr << "NAEquationSolver::Brent: The interval must include a root." << std::endl;\r
  24     return false;\r
  25   }\r
  26     \r
  27   Double_t c = fB;\r
  28   Double_t fc = fb;\r
  29   Double_t d = 0.0;\r
  30   //Double_t fd = 0.0;\r
  31   Double_t e = 0.0;\r
  32   //Double_t fe = 0.0;\r
  33     \r
  34   for(Int_t i=0; i < fMaxIter; i++) {\r
  35     // Rename a,b,c and adjust bounding interval d\r
  36     if(fb*fc > 0.0) {\r
  37       c = fA;\r
  38       fc = fa;\r
  39       d = fB - fA;\r
  40       e = d;\r
  41     }\r
  42     if(Abs(fc) < Abs(fb)) {\r
  43       fA = fB;\r
  44       fB = c;\r
  45       c = fA;\r
  46       fa = fb;\r
  47       fb = fc;\r
  48       fc = fa;\r
  49     }\r
  50     Double_t Tol1 = 2.0*precision*Abs(fB) + 0.5*fTolerance;\r
  51     Double_t xm = 0.5*(c-fB);\r
  52     if(Abs(xm) <= Tol1 || fb == 0.0) {\r
  53       fRoot = fB;\r
  54       return true;\r
  55     }\r
  56     // Inverse quadratic interpolation\r
  57     if(Abs(e) >= Tol1 && Abs(fa) > Abs(fb)) {\r
  58       Double_t s = fb/fa;\r
  59       Double_t p = 0.0;\r
  60       Double_t q = 0.0;\r
  61       if(fA == c) {\r
  62         p = 2.0*xm*s;\r
  63         q = 1.0 - s;\r
  64       } \r
  65       else {\r
  66         q = fa/fc;\r
  67         Double_t r = fb/fc;\r
  68         p = s*(2.0*xm*q*(q-r)-(fB-fA)*(r-1.0));\r
  69         q = (q-1.0)*(r-1.0)*(s-1.0);\r
  70       }\r
  71       // Check bounds\r
  72       if(p > 0.0) q = -q;\r
  73       p = Abs(p);\r
  74       Double_t min1 = 3.0*xm*q-Abs(Tol1*q);\r
  75       Double_t min2 = Abs(e*q);\r
  76       if (2.0*p < Min(min1,min2)) {\r
  77         // Interpolation\r
  78         e = d;\r
  79         d = p/q;\r
  80       } \r
  81       else {\r
  82         // Bisection\r
  83         d = xm;\r
  84         e = d;\r
  85       }\r
  86     } \r
  87     else {\r
  88       // Bounds decreasing too slowly, use bisection\r
  89       d = xm;\r
  90       e = d;\r
  91     }\r
  92     // Move last guess to a \r
  93     fA = fB;\r
  94     fa = fb;\r
  95     if (Abs(d) > Tol1) fB += d;\r
  96     else {\r
  97       if (xm >= 0.0) fB += Abs(Tol1);\r
  98       else fB -= Abs(Tol1);\r
  99     }\r
 100     fb = theFunction(fB);\r
 101   }\r
 102   std::cerr << "NAEquationSolver::Brent: Number of iterations exceeded." << std::endl;\r
 103   return false;\r
 104 }\r
 105 \r
 106 template <class Function>   \r
 107 void NAEquationSolver<Function>::SetIntervalLimits(const Double_t Limit1, const Double_t Limit2) {\r
 108   if(Abs(Limit1-Limit2) <= fTolerance) {\r
 109     std::cerr << "NAEquationSolver::SetIntervalLimits: Interval must be wider than tolerance." << std::endl;\r
 110     return;\r
 111   }\r
 112   if(Limit1 < Limit2) {\r
 113     fA = Limit1;\r
 114     fB = Limit2;\r
 115   } \r
 116   else {\r
 117     fA = Limit2;\r
 118     fB = Limit1;\r
 119   }\r
 120   return;\r
 121 }\r