]>
Commit | Line | Data |
---|---|---|
b1c2e580 | 1 | /* \r |
2 | \r | |
3 | Nikolai Amelin, Ludmila Malinina, Timur Pocheptsov (C) JINR/Dubna\r | |
4 | amelin@sunhe.jinr.ru, malinina@sunhe.jinr.ru, pocheptsov@sunhe.jinr.ru \r | |
5 | November. 2, 2006 \r | |
6 | \r | |
7 | */ \r | |
8 | //This equation solver class is taken from GEANT4 and modified!!\r | |
03896fc4 | 9 | \r |
10 | #include <iostream>\r | |
11 | #include <TMath.h>\r | |
b1c2e580 | 12 | \r |
13 | template <class Function>\r | |
03896fc4 | 14 | Bool_t EquationSolver<Function>::Brent(Function& theFunction) {\r |
b1c2e580 | 15 | const Double_t precision = 3.0e-8;\r |
16 | \r | |
17 | // Check the interval before start\r | |
03896fc4 | 18 | if(fA > fB || TMath::Abs(fA-fB) <= fTolerance) {\r |
19 | std::cerr << "EquationSolver::Brent: The interval must be properly set." << std::endl;\r | |
b1c2e580 | 20 | return false;\r |
21 | }\r | |
22 | Double_t fa = theFunction(fA);\r | |
23 | Double_t fb = theFunction(fB);\r | |
24 | if(fa*fb > 0.0) {\r | |
03896fc4 | 25 | std::cerr << "EquationSolver::Brent: The interval must include a root." << std::endl;\r |
b1c2e580 | 26 | return false;\r |
27 | }\r | |
28 | \r | |
29 | Double_t c = fB;\r | |
30 | Double_t fc = fb;\r | |
31 | Double_t d = 0.0;\r | |
32 | //Double_t fd = 0.0;\r | |
33 | Double_t e = 0.0;\r | |
34 | //Double_t fe = 0.0;\r | |
35 | \r | |
36 | for(Int_t i=0; i < fMaxIter; i++) {\r | |
37 | // Rename a,b,c and adjust bounding interval d\r | |
38 | if(fb*fc > 0.0) {\r | |
39 | c = fA;\r | |
40 | fc = fa;\r | |
41 | d = fB - fA;\r | |
42 | e = d;\r | |
43 | }\r | |
03896fc4 | 44 | if(TMath::Abs(fc) < TMath::Abs(fb)) {\r |
b1c2e580 | 45 | fA = fB;\r |
46 | fB = c;\r | |
47 | c = fA;\r | |
48 | fa = fb;\r | |
49 | fb = fc;\r | |
50 | fc = fa;\r | |
51 | }\r | |
03896fc4 | 52 | Double_t Tol1 = 2.0*precision*TMath::Abs(fB) + 0.5*fTolerance;\r |
b1c2e580 | 53 | Double_t xm = 0.5*(c-fB);\r |
03896fc4 | 54 | if(TMath::Abs(xm) <= Tol1 || fb == 0.0) {\r |
b1c2e580 | 55 | fRoot = fB;\r |
56 | return true;\r | |
57 | }\r | |
58 | // Inverse quadratic interpolation\r | |
03896fc4 | 59 | if(TMath::Abs(e) >= Tol1 && TMath::Abs(fa) > TMath::Abs(fb)) {\r |
b1c2e580 | 60 | Double_t s = fb/fa;\r |
61 | Double_t p = 0.0;\r | |
62 | Double_t q = 0.0;\r | |
63 | if(fA == c) {\r | |
64 | p = 2.0*xm*s;\r | |
65 | q = 1.0 - s;\r | |
66 | } \r | |
67 | else {\r | |
68 | q = fa/fc;\r | |
69 | Double_t r = fb/fc;\r | |
70 | p = s*(2.0*xm*q*(q-r)-(fB-fA)*(r-1.0));\r | |
71 | q = (q-1.0)*(r-1.0)*(s-1.0);\r | |
72 | }\r | |
73 | // Check bounds\r | |
74 | if(p > 0.0) q = -q;\r | |
03896fc4 | 75 | p = TMath::Abs(p);\r |
76 | Double_t min1 = 3.0*xm*q-TMath::Abs(Tol1*q);\r | |
77 | Double_t min2 = TMath::Abs(e*q);\r | |
78 | if (2.0*p < TMath::Min(min1,min2)) {\r | |
b1c2e580 | 79 | // Interpolation\r |
80 | e = d;\r | |
81 | d = p/q;\r | |
82 | } \r | |
83 | else {\r | |
84 | // Bisection\r | |
85 | d = xm;\r | |
86 | e = d;\r | |
87 | }\r | |
88 | } \r | |
89 | else {\r | |
90 | // Bounds decreasing too slowly, use bisection\r | |
91 | d = xm;\r | |
92 | e = d;\r | |
93 | }\r | |
94 | // Move last guess to a \r | |
95 | fA = fB;\r | |
96 | fa = fb;\r | |
03896fc4 | 97 | if (TMath::Abs(d) > Tol1) fB += d;\r |
b1c2e580 | 98 | else {\r |
03896fc4 | 99 | if (xm >= 0.0) fB += TMath::Abs(Tol1);\r |
100 | else fB -= TMath::Abs(Tol1);\r | |
b1c2e580 | 101 | }\r |
102 | fb = theFunction(fB);\r | |
103 | }\r | |
03896fc4 | 104 | std::cerr << "EquationSolver::Brent: Number of iterations exceeded." << std::endl;\r |
b1c2e580 | 105 | return false;\r |
106 | }\r | |
107 | \r | |
108 | template <class Function> \r | |
03896fc4 | 109 | void EquationSolver<Function>::SetIntervalLimits(const Double_t Limit1, const Double_t Limit2) {\r |
110 | if(TMath::Abs(Limit1-Limit2) <= fTolerance) {\r | |
111 | std::cerr << "EquationSolver::SetIntervalLimits: Interval must be wider than tolerance." << std::endl;\r | |
b1c2e580 | 112 | return;\r |
113 | }\r | |
114 | if(Limit1 < Limit2) {\r | |
115 | fA = Limit1;\r | |
116 | fB = Limit2;\r | |
117 | } \r | |
118 | else {\r | |
119 | fA = Limit2;\r | |
120 | fB = Limit1;\r | |
121 | }\r | |
122 | return;\r | |
123 | }\r |