| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 template <class Function>
27 G4bool G4Solver<Function>::Bisection(Function & theFunction)
28 {
29 // Check the interval before start
30 if (a > b || std::abs(a-b) <= tolerance)
31 {
32 G4cerr << "G4Solver::Bisection: The interval must be properly set." << G4endl;
33 return false;
34 }
35 G4double fa = theFunction(a);
36 G4double fb = theFunction(b);
37 if (fa*fb > 0.0)
38 {
39 G4cerr << "G4Solver::Bisection: The interval must include a root." << G4endl;
40 return false;
41 }
42
43 G4double eps=tolerance*(b-a);
44
45
46 // Finding the root
47 for (G4int i = 0; i < MaxIter; i++)
48 {
49 G4double c = (a+b)/2.0;
50 if ((b-a) < eps)
51 {
52 root = c;
53 return true;
54 }
55 G4double fc = theFunction(c);
56 if (fc == 0.0)
57 {
58 root = c;
59 return true;
60 }
61 if (fa*fc < 0.0)
62 {
63 a=c;
64 fa=fc;
65 }
66 else
67 {
68 b=c;
69 fb=fc;
70 }
71 }
72 G4cerr << "G4Solver::Bisection: Exceeded maximum number of iterations without convergence." << G4endl;
73 return false;
74 }
75
76
77 template <class Function>
78 G4bool G4Solver<Function>::RegulaFalsi(Function & theFunction)
79 {
80 // Check the interval before start
81 if (a > b || std::abs(a-b) <= tolerance)
82 {
83 G4cerr << "G4Solver::RegulaFalsi: The interval must be properly set." << G4endl;
84 return false;
85 }
86 G4double fa = theFunction(a);
87 G4double fb = theFunction(b);
88 if (fa*fb > 0.0)
89 {
90 G4cerr << "G4Solver::RegulaFalsi: The interval must include a root." << G4endl;
91 return false;
92 }
93
94 G4double eps=tolerance*(b-a);
95
96
97 // Finding the root
98 for (G4int i = 0; i < MaxIter; i++)
99 {
100 G4double c = (a*fb-b*fa)/(fb-fa);
101 G4double delta = std::min(std::abs(c-a),std::abs(b-c));
102 if (delta < eps)
103 {
104 root = c;
105 return true;
106 }
107 G4double fc = theFunction(c);
108 if (fc == 0.0)
109 {
110 root = c;
111 return true;
112 }
113 if (fa*fc < 0.0)
114 {
115 b=c;
116 fb=fc;
117 }
118 else
119 {
120 a=c;
121 fa=fc;
122 }
123 }
124 G4cerr << "G4Solver::Bisection: Exceeded maximum number of iterations without convergence." << G4endl;
125 return false;
126
127 }
128
129 template <class Function>
130 G4bool G4Solver<Function>::Brent(Function & theFunction)
131 {
132
133 const G4double precision = 3.0e-8;
134
135 // Check the interval before start
136 if (a > b || std::abs(a-b) <= tolerance)
137 {
138 G4cerr << "G4Solver::Brent: The interval must be properly set." << G4endl;
139 return false;
140 }
141 G4double fa = theFunction(a);
142 G4double fb = theFunction(b);
143 if (fa*fb > 0.0)
144 {
145 G4cerr << "G4Solver::Brent: The interval must include a root." << G4endl;
146 return false;
147 }
148
149 G4double c = b;
150 G4double fc = fb;
151 G4double d = 0.0;
152 G4double e = 0.0;
153
154 for (G4int i=0; i < MaxIter; i++)
155 {
156 // Rename a,b,c and adjust bounding interval d
157 if (fb*fc > 0.0)
158 {
159 c = a;
160 fc = fa;
161 d = b - a;
162 e = d;
163 }
164 if (std::abs(fc) < std::abs(fb))
165 {
166 a = b;
167 b = c;
168 c = a;
169 fa = fb;
170 fb = fc;
171 fc = fa;
172 }
173 G4double Tol1 = 2.0*precision*std::abs(b) + 0.5*tolerance;
174 G4double xm = 0.5*(c-b);
175 if (std::abs(xm) <= Tol1 || fb == 0.0)
176 {
177 root = b;
178 return true;
179 }
180 // Inverse quadratic interpolation
181 if (std::abs(e) >= Tol1 && std::abs(fa) > std::abs(fb))
182 {
183 G4double ss = fb/fa;
184 G4double p = 0.0;
185 G4double q = 0.0;
186 if (a == c)
187 {
188 p = 2.0*xm*ss;
189 q = 1.0 - ss;
190 }
191 else
192 {
193 q = fa/fc;
194 G4double r = fb/fc;
195 p = ss*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
196 q = (q-1.0)*(r-1.0)*(ss-1.0);
197 }
198 // Check bounds
199 if (p > 0.0) q = -q;
200 p = std::abs(p);
201 G4double min1 = 3.0*xm*q-std::abs(Tol1*q);
202 G4double min2 = std::abs(e*q);
203 if (2.0*p < std::min(min1,min2))
204 {
205 // Interpolation
206 e = d;
207 d = p/q;
208 }
209 else
210 {
211 // Bisection
212 d = xm;
213 e = d;
214 }
215 }
216 else
217 {
218 // Bounds decreasing too slowly, use bisection
219 d = xm;
220 e = d;
221 }
222 // Move last guess to a
223 a = b;
224 fa = fb;
225 if (std::abs(d) > Tol1) b += d;
226 else
227 {
228 if (xm >= 0.0) b += std::abs(Tol1);
229 else b -= std::abs(Tol1);
230 }
231 fb = theFunction(b);
232 }
233 G4cerr << "G4Solver::Brent: Number of iterations exceeded." << G4endl;
234 return false;
235 }
236
237
238
239 template <class Function>
240 G4bool G4Solver<Function>::Crenshaw(Function & theFunction)
241 {
242 // Check the interval before start
243 if (a > b || std::abs(a-b) <= tolerance)
244 {
245 G4cerr << "G4Solver::Crenshaw: The interval must be properly set." << G4endl;
246 return false;
247 }
248
249 G4double fa = theFunction(a);
250 if (fa == 0.0)
251 {
252 root = a;
253 return true;
254 }
255
256 G4double Mlast = a;
257
258 G4double fb = theFunction(b);
259 if (fb == 0.0)
260 {
261 root = b;
262 return true;
263 }
264
265 if (fa*fb > 0.0)
266 {
267 G4cerr << "G4Solver::Crenshaw: The interval must include a root." << G4endl;
268 return false;
269 }
270
271
272 for (G4int i=0; i < MaxIter; i++)
273 {
274 G4double c = 0.5 * (b + a);
275 G4double fc = theFunction(c);
276 if (fc == 0.0 || std::abs(c - a) < tolerance)
277 {
278 root = c;
279 return true;
280 }
281
282 if (fc * fa > 0.0)
283 {
284 G4double tmp = a;
285 a = b;
286 b = tmp;
287 tmp = fa;
288 fa = fb;
289 fb = tmp;
290 }
291
292 G4double fc0 = fc - fa;
293 G4double fb1 = fb - fc;
294 G4double fb0 = fb - fa;
295 if (fb * fb0 < 2.0 * fc * fc0)
296 {
297 b = c;
298 fb = fc;
299 }
300 else
301 {
302 G4double B = (c - a) / fc0;
303 G4double C = (fc0 - fb1) / (fb1 * fb0);
304 G4double M = a - B * fa * (1.0 - C * fc);
305 G4double fM = theFunction(M);
306 if (fM == 0.0 || std::abs(M - Mlast) < tolerance)
307 {
308 root = M;
309 return true;
310 }
311 Mlast = M;
312 if (fM * fa < 0.0)
313 {
314 b = M;
315 fb = fM;
316 }
317 else
318 {
319 a = M;
320 fa = fM;
321 b = c;
322 fb = fc;
323 }
324 }
325 }
326 return false;
327 }
328
329
330
331
332 template <class Function>
333 void G4Solver<Function>::SetIntervalLimits(const G4double Limit1, const G4double Limit2)
334 {
335 if (std::abs(Limit1-Limit2) <= tolerance)
336 {
337 G4cerr << "G4Solver::SetIntervalLimits: Interval must be wider than tolerance." << G4endl;
338 return;
339 }
340 if (Limit1 < Limit2)
341 {
342 a = Limit1;
343 b = Limit2;
344 }
345 else
346 {
347 a = Limit2;
348 b = Limit1;
349 }
350 return;
351 }
352
353