CommitLineData
7fd59977 1#include <math_BracketedRoot.ixx>
2#include <math_Function.hxx>
3
4// reference algorithme:
5// Brent method
6// numerical recipes in C p 269
7
8 math_BracketedRoot::math_BracketedRoot (math_Function& F,
9 const Standard_Real Bound1,
10 const Standard_Real Bound2,
11 const Standard_Real Tolerance,
12 const Standard_Integer NbIterations,
13 const Standard_Real ZEPS ) {
14
15 Standard_Real Fa,Fc,a,c=0,d=0,e=0;
16 Standard_Real min1,min2,p,q,r,s,tol1,xm;
17 Standard_Boolean Ok;
18
19 a = Bound1;
20 TheRoot = Bound2;
21 Ok = F.Value(a,Fa);
22 Ok = F.Value(TheRoot,TheError);
23 if (Fa*TheError > 0.) { Done = Standard_False;}
24 else {
25 Fc = TheError ;
26 for (NbIter = 1; NbIter <= NbIterations; NbIter++) {
27 if (TheError*Fc > 0.) {
28 c = a; // rename a TheRoot c and adjust bounding interval d
29 Fc = Fa;
30 d = TheRoot - a;
31 e = d;
32 }
33 if ( Abs(Fc) < Abs(Fa) ) {
34 a = TheRoot;
35 TheRoot = c;
36 c = a;
37 Fa = TheError;
38 TheError = Fc;
39 Fc = Fa;
40 }
41 tol1 = 2.*ZEPS * Abs(TheRoot) + 0.5 * Tolerance; // convergence check
42 xm = 0.5 * ( c - TheRoot );
43 if (Abs(xm) <= tol1 || TheError == 0. ) {
44 Done = Standard_True;
45 return;
46 }
47 if (Abs(e) >= tol1 && Abs(Fa) > Abs(TheError) ) {
48 s = TheError / Fa; // attempt inverse quadratic interpolation
49 if (a == c) {
50 p = 2.*xm*s;
51 q = 1. - s;
52 }
53 else {
54 q = Fa / Fc;
55 r = TheError / Fc;
56 p = s * (2.*xm *q * (q - r) - (TheRoot - a)*(r - 1.));
57 q = (q -1.) * (r - 1.) * (s - 1.);
58 }
59 if ( p > 0. ) { q = -q;} // check whether in bounds
60 p = Abs(p);
61 min1 = 3.* xm* q - Abs(tol1 *q);
62 min2 = Abs(e * q);
63 if (2.* p < (min1 < min2 ? min1 : min2) ) {
64 e = d ; // accept interpolation
65 d = p / q;
66 }
67 else {
68 d = xm; // interpolation failed,use bissection
69 e = d;
70 }
71 }
72 else { // bounds decreasing too slowly ,use bissection
73 d = xm;
74 e =d;
75 }
76 a = TheRoot ; // move last best guess to a
77 Fa = TheError;
78 if (Abs(d) > tol1) { // evaluate new trial root
79 TheRoot += d;
80 }
81 else {
82 TheRoot += (xm > 0. ? Abs(tol1) : -Abs(tol1));
83 }
84 Ok = F.Value(TheRoot,TheError);
85 }
86 Done = Standard_False;
87 }
88 }
89
90
91 void math_BracketedRoot::Dump(Standard_OStream& o) const {
92
93 o << "math_BracketedRoot ";
94 if(Done) {
95 o << " Status = Done \n";
96 o << " Number of iterations = " << NbIter << endl;
97 o << " The Root is: " << TheRoot << endl;
98 o << " The value at the root is: " << TheError << endl;
99 }
100 else {
101 o << " Status = not Done \n";
102 }
103}