Commit | Line | Data |
---|---|---|
b311480e | 1 | // Copyright (c) 1997-1999 Matra Datavision |
973c2be1 | 2 | // Copyright (c) 1999-2014 OPEN CASCADE SAS |
b311480e | 3 | // |
973c2be1 | 4 | // This file is part of Open CASCADE Technology software library. |
b311480e | 5 | // |
d5f74e42 | 6 | // This library is free software; you can redistribute it and/or modify it under |
7 | // the terms of the GNU Lesser General Public License version 2.1 as published | |
973c2be1 | 8 | // by the Free Software Foundation, with special exception defined in the file |
9 | // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT | |
10 | // distribution for complete text of the license and disclaimer of any warranty. | |
b311480e | 11 | // |
973c2be1 | 12 | // Alternatively, this file may be used under the terms of Open CASCADE |
13 | // commercial license or contractual agreement. | |
b311480e | 14 | |
42cf5bc1 | 15 | |
16 | #include <math_BracketedRoot.hxx> | |
7fd59977 | 17 | #include <math_Function.hxx> |
42cf5bc1 | 18 | #include <StdFail_NotDone.hxx> |
7fd59977 | 19 | |
20 | // reference algorithme: | |
21 | // Brent method | |
22 | // numerical recipes in C p 269 | |
b311480e | 23 | math_BracketedRoot::math_BracketedRoot (math_Function& F, |
7fd59977 | 24 | const Standard_Real Bound1, |
25 | const Standard_Real Bound2, | |
26 | const Standard_Real Tolerance, | |
27 | const Standard_Integer NbIterations, | |
28 | const Standard_Real ZEPS ) { | |
29 | ||
30 | Standard_Real Fa,Fc,a,c=0,d=0,e=0; | |
31 | Standard_Real min1,min2,p,q,r,s,tol1,xm; | |
7fd59977 | 32 | |
33 | a = Bound1; | |
34 | TheRoot = Bound2; | |
96a95605 DB |
35 | F.Value(a,Fa); |
36 | F.Value(TheRoot,TheError); | |
7fd59977 | 37 | if (Fa*TheError > 0.) { Done = Standard_False;} |
38 | else { | |
39 | Fc = TheError ; | |
40 | for (NbIter = 1; NbIter <= NbIterations; NbIter++) { | |
41 | if (TheError*Fc > 0.) { | |
42 | c = a; // rename a TheRoot c and adjust bounding interval d | |
43 | Fc = Fa; | |
44 | d = TheRoot - a; | |
45 | e = d; | |
46 | } | |
47 | if ( Abs(Fc) < Abs(Fa) ) { | |
48 | a = TheRoot; | |
49 | TheRoot = c; | |
50 | c = a; | |
51 | Fa = TheError; | |
52 | TheError = Fc; | |
53 | Fc = Fa; | |
54 | } | |
55 | tol1 = 2.*ZEPS * Abs(TheRoot) + 0.5 * Tolerance; // convergence check | |
56 | xm = 0.5 * ( c - TheRoot ); | |
57 | if (Abs(xm) <= tol1 || TheError == 0. ) { | |
58 | Done = Standard_True; | |
59 | return; | |
60 | } | |
61 | if (Abs(e) >= tol1 && Abs(Fa) > Abs(TheError) ) { | |
62 | s = TheError / Fa; // attempt inverse quadratic interpolation | |
63 | if (a == c) { | |
64 | p = 2.*xm*s; | |
65 | q = 1. - s; | |
66 | } | |
67 | else { | |
68 | q = Fa / Fc; | |
69 | r = TheError / Fc; | |
70 | p = s * (2.*xm *q * (q - r) - (TheRoot - a)*(r - 1.)); | |
71 | q = (q -1.) * (r - 1.) * (s - 1.); | |
72 | } | |
73 | if ( p > 0. ) { q = -q;} // check whether in bounds | |
74 | p = Abs(p); | |
75 | min1 = 3.* xm* q - Abs(tol1 *q); | |
76 | min2 = Abs(e * q); | |
77 | if (2.* p < (min1 < min2 ? min1 : min2) ) { | |
78 | e = d ; // accept interpolation | |
79 | d = p / q; | |
80 | } | |
81 | else { | |
82 | d = xm; // interpolation failed,use bissection | |
83 | e = d; | |
84 | } | |
85 | } | |
86 | else { // bounds decreasing too slowly ,use bissection | |
87 | d = xm; | |
88 | e =d; | |
89 | } | |
90 | a = TheRoot ; // move last best guess to a | |
91 | Fa = TheError; | |
92 | if (Abs(d) > tol1) { // evaluate new trial root | |
93 | TheRoot += d; | |
94 | } | |
95 | else { | |
96 | TheRoot += (xm > 0. ? Abs(tol1) : -Abs(tol1)); | |
97 | } | |
96a95605 | 98 | F.Value(TheRoot,TheError); |
7fd59977 | 99 | } |
100 | Done = Standard_False; | |
101 | } | |
102 | } | |
103 | ||
104 | ||
105 | void math_BracketedRoot::Dump(Standard_OStream& o) const { | |
106 | ||
107 | o << "math_BracketedRoot "; | |
108 | if(Done) { | |
109 | o << " Status = Done \n"; | |
04232180 | 110 | o << " Number of iterations = " << NbIter << std::endl; |
111 | o << " The Root is: " << TheRoot << std::endl; | |
112 | o << " The value at the root is: " << TheError << std::endl; | |
7fd59977 | 113 | } |
114 | else { | |
115 | o << " Status = not Done \n"; | |
116 | } | |
117 | } |