include/Geant4/G4Solver.icc

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 template <class Function>
0027 G4bool G4Solver<Function>::Bisection(Function & theFunction)
0028 {
0029     // Check the interval before start
0030     if (a > b || std::abs(a-b) <= tolerance)
0031     {
0032     G4cerr << "G4Solver::Bisection: The interval must be properly set." << G4endl;
0033     return false;
0034     }
0035     G4double fa = theFunction(a);
0036     G4double fb = theFunction(b);
0037     if (fa*fb > 0.0)
0038     {
0039     G4cerr << "G4Solver::Bisection: The interval must include a root." << G4endl;
0040     return false;
0041     }
0042
0043     G4double eps=tolerance*(b-a);
0044
0045
0046     // Finding the root
0047     for (G4int i = 0; i < MaxIter; i++)
0048     {
0049     G4double c = (a+b)/2.0;
0050     if ((b-a) < eps)
0051     {
0052         root = c;
0053         return true;
0054     }
0055     G4double fc = theFunction(c);
0056     if (fc == 0.0)
0057     {
0058         root = c;
0059         return true;
0060     }
0061     if (fa*fc < 0.0)
0062     {
0063         a=c;
0064         fa=fc;
0065     }
0066     else
0067     {
0068         b=c;
0069         fb=fc;
0070     }
0071     }
0072     G4cerr << "G4Solver::Bisection: Exceeded maximum number of iterations without convergence." << G4endl;
0073     return false;
0074 }
0075
0076
0077 template <class Function>
0078 G4bool G4Solver<Function>::RegulaFalsi(Function & theFunction)
0079 {
0080     // Check the interval before start
0081     if (a > b || std::abs(a-b) <= tolerance)
0082     {
0083     G4cerr << "G4Solver::RegulaFalsi: The interval must be properly set." << G4endl;
0084     return false;
0085     }
0086     G4double fa = theFunction(a);
0087     G4double fb = theFunction(b);
0088     if (fa*fb > 0.0)
0089     {
0090     G4cerr << "G4Solver::RegulaFalsi: The interval must include a root." << G4endl;
0091     return false;
0092     }
0093
0094     G4double eps=tolerance*(b-a);
0095
0096
0097     // Finding the root
0098     for (G4int i = 0; i < MaxIter; i++)
0099     {
0100     G4double c = (a*fb-b*fa)/(fb-fa);
0101     G4double delta = std::min(std::abs(c-a),std::abs(b-c));
0102     if (delta < eps)
0103     {
0104         root = c;
0105         return true;
0106     }
0107     G4double fc = theFunction(c);
0108     if (fc == 0.0)
0109     {
0110         root = c;
0111         return true;
0112     }
0113     if (fa*fc < 0.0)
0114     {
0115         b=c;
0116         fb=fc;
0117     }
0118     else
0119     {
0120         a=c;
0121         fa=fc;
0122     }
0123     }
0124     G4cerr << "G4Solver::Bisection: Exceeded maximum number of iterations without convergence." << G4endl;
0125     return false;
0126
0127 }
0128
0129 template <class Function>
0130 G4bool G4Solver<Function>::Brent(Function & theFunction)
0131 {
0132
0133     const G4double precision = 3.0e-8;
0134
0135     // Check the interval before start
0136     if (a > b || std::abs(a-b) <= tolerance)
0137     {
0138     G4cerr << "G4Solver::Brent: The interval must be properly set." << G4endl;
0139     return false;
0140     }
0141     G4double fa = theFunction(a);
0142     G4double fb = theFunction(b);
0143     if (fa*fb > 0.0)
0144     {
0145     G4cerr << "G4Solver::Brent: The interval must include a root." << G4endl;
0146     return false;
0147     }
0148
0149     G4double c = b;
0150     G4double fc = fb;
0151     G4double d = 0.0;
0152     G4double e = 0.0;
0153
0154     for (G4int i=0; i < MaxIter; i++)
0155     {
0156     // Rename a,b,c and adjust bounding interval d
0157     if (fb*fc > 0.0)
0158     {
0159         c = a;
0160         fc = fa;
0161         d = b - a;
0162         e = d;
0163     }
0164     if (std::abs(fc) < std::abs(fb))
0165     {
0166         a = b;
0167         b = c;
0168         c = a;
0169         fa = fb;
0170         fb = fc;
0171         fc = fa;
0172     }
0173     G4double Tol1 = 2.0*precision*std::abs(b) + 0.5*tolerance;
0174     G4double xm = 0.5*(c-b);
0175     if (std::abs(xm) <= Tol1 || fb == 0.0)
0176     {
0177         root = b;
0178         return true;
0179     }
0180     // Inverse quadratic interpolation
0181     if (std::abs(e) >= Tol1 && std::abs(fa) > std::abs(fb))
0182     {
0183         G4double ss = fb/fa;
0184         G4double p = 0.0;
0185         G4double q = 0.0;
0186         if (a == c)
0187         {
0188         p = 2.0*xm*ss;
0189         q = 1.0 - ss;
0190         }
0191         else
0192         {
0193         q = fa/fc;
0194         G4double r = fb/fc;
0195         p = ss*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
0196         q = (q-1.0)*(r-1.0)*(ss-1.0);
0197         }
0198         // Check bounds
0199         if (p > 0.0) q = -q;
0200         p = std::abs(p);
0201         G4double min1 = 3.0*xm*q-std::abs(Tol1*q);
0202         G4double min2 = std::abs(e*q);
0203         if (2.0*p < std::min(min1,min2))
0204         {
0205         // Interpolation
0206         e = d;
0207         d = p/q;
0208         }
0209         else
0210         {
0211         // Bisection
0212         d = xm;
0213         e = d;
0214         }
0215     }
0216     else
0217     {
0218         // Bounds decreasing too slowly, use bisection
0219         d = xm;
0220         e = d;
0221     }
0222     // Move last guess to a
0223     a = b;
0224     fa = fb;
0225     if (std::abs(d) > Tol1) b += d;
0226     else
0227     {
0228         if (xm >= 0.0) b += std::abs(Tol1);
0229         else b -= std::abs(Tol1);
0230     }
0231     fb = theFunction(b);
0232     }
0233     G4cerr << "G4Solver::Brent: Number of iterations exceeded." << G4endl;
0234     return false;
0235 }
0236
0237
0238
0239 template <class Function>
0240 G4bool G4Solver<Function>::Crenshaw(Function & theFunction)
0241 {
0242     // Check the interval before start
0243     if (a > b || std::abs(a-b) <= tolerance)
0244     {
0245     G4cerr << "G4Solver::Crenshaw: The interval must be properly set." << G4endl;
0246     return false;
0247     }
0248
0249     G4double fa = theFunction(a);
0250     if (fa == 0.0)
0251     {
0252     root = a;
0253     return true;
0254     }
0255
0256     G4double Mlast = a;
0257
0258     G4double fb = theFunction(b);
0259     if (fb == 0.0)
0260     {
0261     root = b;
0262     return true;
0263     }
0264
0265     if (fa*fb > 0.0)
0266     {
0267     G4cerr << "G4Solver::Crenshaw: The interval must include a root." << G4endl;
0268     return false;
0269     }
0270
0271
0272     for (G4int i=0; i < MaxIter; i++)
0273       {
0274     G4double c = 0.5 * (b + a);
0275     G4double fc = theFunction(c);
0276     if (fc == 0.0 || std::abs(c - a) < tolerance)
0277       {
0278         root = c;
0279         return true;
0280       }
0281
0282     if (fc * fa  > 0.0)
0283       {
0284         G4double tmp = a;
0285         a = b;
0286         b = tmp;
0287         tmp = fa;
0288         fa = fb;
0289         fb = tmp;
0290       }
0291
0292     G4double fc0 = fc - fa;
0293     G4double fb1 = fb - fc;
0294     G4double fb0 = fb - fa;
0295     if (fb * fb0 < 2.0 * fc * fc0)
0296       {
0297         b = c;
0298         fb = fc;
0299       }
0300     else
0301       {
0302         G4double B = (c - a) / fc0;
0303         G4double C = (fc0 - fb1) / (fb1 * fb0);
0304         G4double M = a - B * fa * (1.0 - C * fc);
0305         G4double fM = theFunction(M);
0306         if (fM == 0.0 || std::abs(M - Mlast) < tolerance)
0307           {
0308         root = M;
0309         return true;
0310           }
0311         Mlast = M;
0312         if (fM * fa < 0.0)
0313           {
0314         b = M;
0315         fb = fM;
0316           }
0317         else
0318         {
0319         a = M;
0320         fa = fM;
0321         b = c;
0322         fb = fc;
0323         }
0324     }
0325     }
0326     return false;
0327 }
0328
0329
0330
0331
0332 template <class Function>
0333 void G4Solver<Function>::SetIntervalLimits(const G4double Limit1, const G4double Limit2)
0334 {
0335     if (std::abs(Limit1-Limit2) <= tolerance)
0336     {
0337     G4cerr << "G4Solver::SetIntervalLimits: Interval must be wider than tolerance." << G4endl;
0338     return;
0339     }
0340     if (Limit1 < Limit2)
0341     {
0342     a = Limit1;
0343     b = Limit2;
0344     }
0345     else
0346     {
0347     a = Limit2;
0348     b = Limit1;
0349     }
0350     return;
0351 }
0352