source: trunk/source/processes/hadronic/models/de_excitation/multifragmentation/include/G4Solver.icc @ 1007

Last change on this file since 1007 was 819, checked in by garnier, 16 years ago

import all except CVS
File size: 7.6 KB
Line 
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//
26template <class Function>
27G4bool 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
77template <class Function>
78G4bool 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
129template <class Function>
130G4bool 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 s = fb/fa;
184            G4double p = 0.0;
185            G4double q = 0.0;
186            if (a == c)
187            {
188                p = 2.0*xm*s;
189                q = 1.0 - s;
190            }
191            else
192            {
193                q = fa/fc;
194                G4double r = fb/fc;
195                p = s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
196                q = (q-1.0)*(r-1.0)*(s-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
239template <class Function>
240G4bool 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
332template <class Function>       
333void 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
Note: See TracBrowser for help on using the repository browser.