G4Solver.icc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 template <class Function>
00027 G4bool G4Solver<Function>::Bisection(Function & theFunction)
00028 {
00029     // Check the interval before start
00030     if (a > b || std::abs(a-b) <= tolerance) 
00031     {
00032         G4cerr << "G4Solver::Bisection: The interval must be properly set." << G4endl;
00033         return false;
00034     }
00035     G4double fa = theFunction(a);
00036     G4double fb = theFunction(b);
00037     if (fa*fb > 0.0) 
00038     {
00039         G4cerr << "G4Solver::Bisection: The interval must include a root." << G4endl;
00040         return false;
00041     }
00042     
00043     G4double eps=tolerance*(b-a);
00044     
00045     
00046     // Finding the root
00047     for (G4int i = 0; i < MaxIter; i++) 
00048     {
00049         G4double c = (a+b)/2.0;
00050         if ((b-a) < eps) 
00051         {
00052             root = c;
00053             return true;
00054         }
00055         G4double fc = theFunction(c);
00056         if (fc == 0.0) 
00057         {
00058             root = c;
00059             return true;
00060         }
00061         if (fa*fc < 0.0) 
00062         {
00063             a=c;
00064             fa=fc;
00065         } 
00066         else 
00067         {
00068             b=c;
00069             fb=fc;
00070         }
00071     }
00072     G4cerr << "G4Solver::Bisection: Excedded maximum number of iterations whithout convegence." << G4endl;
00073     return false;
00074 }
00075 
00076 
00077 template <class Function>
00078 G4bool G4Solver<Function>::RegulaFalsi(Function & theFunction)
00079 {
00080     // Check the interval before start
00081     if (a > b || std::abs(a-b) <= tolerance) 
00082     {
00083         G4cerr << "G4Solver::RegulaFalsi: The interval must be properly set." << G4endl;
00084         return false;
00085     }
00086     G4double fa = theFunction(a);
00087     G4double fb = theFunction(b);
00088     if (fa*fb > 0.0) 
00089     {
00090         G4cerr << "G4Solver::RegulaFalsi: The interval must include a root." << G4endl;
00091         return false;
00092     }
00093     
00094     G4double eps=tolerance*(b-a);
00095         
00096         
00097     // Finding the root
00098     for (G4int i = 0; i < MaxIter; i++) 
00099     {
00100         G4double c = (a*fb-b*fa)/(fb-fa);
00101         G4double delta = std::min(std::abs(c-a),std::abs(b-c));
00102         if (delta < eps) 
00103         {
00104             root = c;
00105             return true;
00106         }
00107         G4double fc = theFunction(c);
00108         if (fc == 0.0) 
00109         {
00110             root = c;
00111             return true;
00112         }
00113         if (fa*fc < 0.0) 
00114         {
00115             b=c;
00116             fb=fc;
00117         } 
00118         else 
00119         {
00120             a=c;
00121             fa=fc;
00122         }
00123     }
00124     G4cerr << "G4Solver::Bisection: Excedded maximum number of iterations whithout convegence." << G4endl;
00125     return false;
00126 
00127 }       
00128 
00129 template <class Function>
00130 G4bool G4Solver<Function>::Brent(Function & theFunction)
00131 {
00132     
00133     const G4double precision = 3.0e-8;
00134     
00135     // Check the interval before start
00136     if (a > b || std::abs(a-b) <= tolerance) 
00137     {
00138         G4cerr << "G4Solver::Brent: The interval must be properly set." << G4endl;
00139         return false;
00140     }
00141     G4double fa = theFunction(a);
00142     G4double fb = theFunction(b);
00143     if (fa*fb > 0.0) 
00144     {
00145         G4cerr << "G4Solver::Brent: The interval must include a root." << G4endl;
00146         return false;
00147     }
00148     
00149     G4double c = b;
00150     G4double fc = fb;
00151     G4double d = 0.0;
00152     G4double e = 0.0;
00153     
00154     for (G4int i=0; i < MaxIter; i++) 
00155     {
00156         // Rename a,b,c and adjust bounding interval d
00157         if (fb*fc > 0.0) 
00158         {
00159             c = a;
00160             fc = fa;
00161             d = b - a;
00162             e = d;
00163         }
00164         if (std::abs(fc) < std::abs(fb)) 
00165         {
00166             a = b;
00167             b = c;
00168             c = a;
00169             fa = fb;
00170             fb = fc;
00171             fc = fa;
00172         }
00173         G4double Tol1 = 2.0*precision*std::abs(b) + 0.5*tolerance;
00174         G4double xm = 0.5*(c-b);
00175         if (std::abs(xm) <= Tol1 || fb == 0.0) 
00176         {
00177             root = b;
00178             return true;
00179         }
00180         // Inverse quadratic interpolation
00181         if (std::abs(e) >= Tol1 && std::abs(fa) > std::abs(fb)) 
00182         {
00183             G4double ss = fb/fa;
00184             G4double p = 0.0;
00185             G4double q = 0.0;
00186             if (a == c) 
00187             {
00188                 p = 2.0*xm*ss;
00189                 q = 1.0 - ss;
00190             } 
00191             else 
00192             {
00193                 q = fa/fc;
00194                 G4double r = fb/fc;
00195                 p = ss*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
00196                 q = (q-1.0)*(r-1.0)*(ss-1.0);
00197             }
00198             // Check bounds
00199             if (p > 0.0) q = -q;
00200             p = std::abs(p);
00201             G4double min1 = 3.0*xm*q-std::abs(Tol1*q);
00202             G4double min2 = std::abs(e*q);
00203             if (2.0*p < std::min(min1,min2)) 
00204             {
00205                 // Interpolation
00206                 e = d;
00207                 d = p/q;
00208             } 
00209             else 
00210             {
00211                 // Bisection
00212                 d = xm;
00213                 e = d;
00214             }
00215         } 
00216         else 
00217         {
00218             // Bounds decreasing too slowly, use bisection
00219             d = xm;
00220             e = d;
00221         }
00222         // Move last guess to a 
00223         a = b;
00224         fa = fb;
00225         if (std::abs(d) > Tol1) b += d;
00226         else 
00227         {
00228             if (xm >= 0.0) b += std::abs(Tol1);
00229             else b -= std::abs(Tol1);
00230         }
00231         fb = theFunction(b);
00232     }
00233     G4cerr << "G4Solver::Brent: Number of iterations exceeded." << G4endl;
00234     return false;
00235 }
00236 
00237 
00238 
00239 template <class Function>
00240 G4bool G4Solver<Function>::Crenshaw(Function & theFunction)
00241 {
00242     // Check the interval before start
00243     if (a > b || std::abs(a-b) <= tolerance) 
00244     {
00245         G4cerr << "G4Solver::Crenshaw: The interval must be properly set." << G4endl;
00246         return false;
00247     }
00248 
00249     G4double fa = theFunction(a);
00250     if (fa == 0.0) 
00251     {
00252         root = a;
00253         return true;
00254     }
00255 
00256     G4double Mlast = a;
00257 
00258     G4double fb = theFunction(b);
00259     if (fb == 0.0)
00260     {
00261         root = b;
00262         return true;
00263     }
00264 
00265     if (fa*fb > 0.0) 
00266     {
00267         G4cerr << "G4Solver::Crenshaw: The interval must include a root." << G4endl;
00268         return false;
00269     }
00270 
00271     
00272     for (G4int i=0; i < MaxIter; i++) 
00273       {
00274         G4double c = 0.5 * (b + a);
00275         G4double fc = theFunction(c);
00276         if (fc == 0.0 || std::abs(c - a) < tolerance) 
00277           {
00278             root = c;
00279             return true;
00280           }
00281 
00282         if (fc * fa  > 0.0)
00283           {
00284             G4double tmp = a;
00285             a = b;
00286             b = tmp;
00287             tmp = fa;
00288             fa = fb;
00289             fb = tmp;
00290           }
00291         
00292         G4double fc0 = fc - fa;
00293         G4double fb1 = fb - fc;
00294         G4double fb0 = fb - fa;
00295         if (fb * fb0 < 2.0 * fc * fc0)
00296           {
00297             b = c;
00298             fb = fc;
00299           }
00300         else 
00301           {
00302             G4double B = (c - a) / fc0;
00303             G4double C = (fc0 - fb1) / (fb1 * fb0);
00304             G4double M = a - B * fa * (1.0 - C * fc);
00305             G4double fM = theFunction(M);
00306             if (fM == 0.0 || std::abs(M - Mlast) < tolerance)
00307               {
00308                 root = M;
00309                 return true;
00310               }
00311             Mlast = M;
00312             if (fM * fa < 0.0)
00313               {
00314                 b = M;
00315                 fb = fM;
00316               }
00317             else 
00318             {
00319                 a = M;
00320                 fa = fM;
00321                 b = c;
00322                 fb = fc;
00323             }
00324         }
00325     }
00326     return false;
00327 }
00328 
00329 
00330 
00331 
00332 template <class Function>       
00333 void G4Solver<Function>::SetIntervalLimits(const G4double Limit1, const G4double Limit2)
00334 {
00335     if (std::abs(Limit1-Limit2) <= tolerance) 
00336     {
00337         G4cerr << "G4Solver::SetIntervalLimits: Interval must be wider than tolerance." << G4endl;
00338         return;
00339     }
00340     if (Limit1 < Limit2) 
00341     {
00342         a = Limit1;
00343         b = Limit2;
00344     } 
00345     else 
00346     {
00347         a = Limit2;
00348         b = Limit1;
00349     }
00350     return;
00351 }
00352 

Generated on Mon May 27 17:49:52 2013 for Geant4 by  doxygen 1.4.7