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