Geant4  10.02
G4Solver.icc
Go to the documentation of this file.
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: Excedded maximum number of iterations whithout convegence." << 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: Excedded maximum number of iterations whithout convegence." << 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