00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 template <class Function>
00027 G4bool G4Solver<Function>::Bisection(Function & theFunction)
00028 {
00029
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
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
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
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
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
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
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
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
00206 e = d;
00207 d = p/q;
00208 }
00209 else
00210 {
00211
00212 d = xm;
00213 e = d;
00214 }
00215 }
00216 else
00217 {
00218
00219 d = xm;
00220 e = d;
00221 }
00222
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
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