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
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
|
//static const char* sccsid = "@(#)math_BissecNewton.cxx 3.2 95/01/10"; // Do not delete this line. Used by sccs.
#include <math_BissecNewton.ixx>
#include <math_FunctionWithDerivative.hxx>
void math_BissecNewton::Perform(math_FunctionWithDerivative& F,
const Standard_Real Bound1,
const Standard_Real Bound2,
const Standard_Integer NbIterations)
{
Standard_Boolean GOOD;
Standard_Integer j;
Standard_Real dxold, fh, fl;
Standard_Real swap, temp, xh, xl;
GOOD = F.Values(Bound1, fl, df);
if(!GOOD) {
Done = Standard_False;
TheStatus = math_FunctionError;
return;
}
GOOD = F.Values(Bound2, fh, df);
if(!GOOD) {
Done = Standard_False;
TheStatus = math_FunctionError;
return;
}
// Modified by Sergey KHROMOV - Wed Jan 22 12:06:45 2003 Begin
Standard_Real aFTol = RealEpsilon();
// if(fl * fh >= 0.0) {
if(fl * fh > aFTol*aFTol) {
Done = Standard_False;
TheStatus = math_NotBracketed;
return;
}
// if(fl < 0.0) {
if(fl < -aFTol || (fl < aFTol && fh < -aFTol)) {
xl = Bound1;
xh = Bound2;
}
else {
xl = Bound2;
xh = Bound1;
swap = fl;
fl = fh;
fh = swap;
}
// Modified by Sergey KHROMOV - Wed Jan 22 12:06:49 2003 End
x = 0.5 * (Bound1 + Bound2);
dxold = fabs(Bound2 - Bound1);
dx = dxold;
GOOD = F.Values(x, f, df);
if(!GOOD) {
Done = Standard_False;
TheStatus = math_FunctionError;
return;
}
for(j = 1; j <= NbIterations; j++) {
if((((x - xh) * df - f) * ((x - xl) * df - f) >= 0.0)
|| (fabs(2.0 * f) > fabs(dxold * df))) {
dxold = dx;
dx = 0.5 * (xh - xl);
x = xl + dx;
if(xl == x) {
TheStatus = math_OK;
Done = Standard_True;
return;
}
}
else {
dxold = dx;
dx = f / df;
temp = x;
x -= dx;
if(temp == x) {
TheStatus = math_OK;
Done = Standard_True;
return;
}
}
if(IsSolutionReached(F)) {
TheStatus = math_OK;
Done = Standard_True;
return;
}
GOOD = F.Values(x, f, df);
if(!GOOD) {
Done = Standard_False;
TheStatus = math_FunctionError;
return;
}
if(f < 0.0) {
xl = x;
fl = f;
}
else if(f > 0.0) {
xh = x;
fh = f;
}
else {
TheStatus = math_OK;
Done = Standard_True;
return;
}
}
TheStatus = math_TooManyIterations;
Done = Standard_False;
return;
}
Standard_Boolean math_BissecNewton::IsSolutionReached
//(math_FunctionWithDerivative& F)
(math_FunctionWithDerivative& )
{
return Abs(dx) <= XTol;
}
math_BissecNewton::math_BissecNewton(math_FunctionWithDerivative& F,
const Standard_Real Bound1,
const Standard_Real Bound2,
const Standard_Real TolX,
const Standard_Integer NbIterations) {
XTol = TolX;
Perform(F, Bound1, Bound2, NbIterations);
}
void math_BissecNewton::Dump(Standard_OStream& o) const {
o << "math_BissecNewton ";
if(Done) {
o << " Status = Done \n";
o << " The Root is: " << x << endl;
o << " The value at this Root is: " << f << endl;
}
else {
o << " Status = not Done \n";
}
}
|
