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