7fd59977 |
1 | #include <math_BissecNewton.ixx> |
2 | #include <math_FunctionWithDerivative.hxx> |
3 | |
4 | |
5 | void math_BissecNewton::Perform(math_FunctionWithDerivative& F, |
6 | const Standard_Real Bound1, |
7 | const Standard_Real Bound2, |
8 | const Standard_Integer NbIterations) |
9 | { |
10 | Standard_Boolean GOOD; |
11 | Standard_Integer j; |
12 | Standard_Real dxold, fh, fl; |
13 | Standard_Real swap, temp, xh, xl; |
14 | |
15 | GOOD = F.Values(Bound1, fl, df); |
16 | if(!GOOD) { |
17 | Done = Standard_False; |
18 | TheStatus = math_FunctionError; |
19 | return; |
20 | } |
21 | GOOD = F.Values(Bound2, fh, df); |
22 | if(!GOOD) { |
23 | Done = Standard_False; |
24 | TheStatus = math_FunctionError; |
25 | return; |
26 | } |
27 | // Modified by Sergey KHROMOV - Wed Jan 22 12:06:45 2003 Begin |
28 | Standard_Real aFTol = RealEpsilon(); |
29 | |
30 | // if(fl * fh >= 0.0) { |
31 | if(fl * fh > aFTol*aFTol) { |
32 | Done = Standard_False; |
33 | TheStatus = math_NotBracketed; |
34 | return; |
35 | } |
36 | // if(fl < 0.0) { |
37 | if(fl < -aFTol || (fl < aFTol && fh < -aFTol)) { |
38 | xl = Bound1; |
39 | xh = Bound2; |
40 | } |
41 | else { |
42 | xl = Bound2; |
43 | xh = Bound1; |
44 | swap = fl; |
45 | fl = fh; |
46 | fh = swap; |
47 | } |
48 | // Modified by Sergey KHROMOV - Wed Jan 22 12:06:49 2003 End |
49 | x = 0.5 * (Bound1 + Bound2); |
50 | dxold = fabs(Bound2 - Bound1); |
51 | dx = dxold; |
52 | GOOD = F.Values(x, f, df); |
53 | if(!GOOD) { |
54 | Done = Standard_False; |
55 | TheStatus = math_FunctionError; |
56 | return; |
57 | } |
58 | for(j = 1; j <= NbIterations; j++) { |
59 | if((((x - xh) * df - f) * ((x - xl) * df - f) >= 0.0) |
60 | || (fabs(2.0 * f) > fabs(dxold * df))) { |
61 | dxold = dx; |
62 | dx = 0.5 * (xh - xl); |
63 | x = xl + dx; |
64 | if(xl == x) { |
65 | TheStatus = math_OK; |
66 | Done = Standard_True; |
67 | return; |
68 | } |
69 | } |
70 | else { |
71 | dxold = dx; |
72 | dx = f / df; |
73 | temp = x; |
74 | x -= dx; |
75 | if(temp == x) { |
76 | TheStatus = math_OK; |
77 | Done = Standard_True; |
78 | return; |
79 | } |
80 | } |
81 | if(IsSolutionReached(F)) { |
82 | TheStatus = math_OK; |
83 | Done = Standard_True; |
84 | return; |
85 | } |
86 | GOOD = F.Values(x, f, df); |
87 | if(!GOOD) { |
88 | Done = Standard_False; |
89 | TheStatus = math_FunctionError; |
90 | return; |
91 | } |
92 | if(f < 0.0) { |
93 | xl = x; |
94 | fl = f; |
95 | } |
96 | else if(f > 0.0) { |
97 | xh = x; |
98 | fh = f; |
99 | } |
100 | else { |
101 | TheStatus = math_OK; |
102 | Done = Standard_True; |
103 | return; |
104 | } |
105 | } |
106 | TheStatus = math_TooManyIterations; |
107 | Done = Standard_False; |
108 | return; |
109 | } |
110 | |
111 | Standard_Boolean math_BissecNewton::IsSolutionReached |
112 | //(math_FunctionWithDerivative& F) |
113 | (math_FunctionWithDerivative& ) |
114 | { |
115 | return Abs(dx) <= XTol; |
116 | } |
117 | |
118 | |
119 | math_BissecNewton::math_BissecNewton(math_FunctionWithDerivative& F, |
120 | const Standard_Real Bound1, |
121 | const Standard_Real Bound2, |
122 | const Standard_Real TolX, |
123 | const Standard_Integer NbIterations) { |
124 | |
125 | XTol = TolX; |
126 | Perform(F, Bound1, Bound2, NbIterations); |
127 | } |
128 | |
129 | |
130 | void math_BissecNewton::Dump(Standard_OStream& o) const { |
131 | |
132 | o << "math_BissecNewton "; |
133 | if(Done) { |
134 | o << " Status = Done \n"; |
135 | o << " The Root is: " << x << endl; |
136 | o << " The value at this Root is: " << f << endl; |
137 | } |
138 | else { |
139 | o << " Status = not Done \n"; |
140 | } |
141 | } |
142 | |