1 // Copyright (c) 1997-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
4 // This file is part of Open CASCADE Technology software library.
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
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.
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
16 #include <math_BracketedRoot.hxx>
17 #include <math_Function.hxx>
18 #include <StdFail_NotDone.hxx>
20 // reference algorithme:
22 // numerical recipes in C p 269
23 math_BracketedRoot::math_BracketedRoot (math_Function& F,
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 ) {
30 Standard_Real Fa,Fc,a,c=0,d=0,e=0;
31 Standard_Real min1,min2,p,q,r,s,tol1,xm;
36 F.Value(TheRoot,TheError);
37 if (Fa*TheError > 0.) { Done = Standard_False;}
40 for (NbIter = 1; NbIter <= NbIterations; NbIter++) {
41 if (TheError*Fc > 0.) {
42 c = a; // rename a TheRoot c and adjust bounding interval d
47 if ( Abs(Fc) < Abs(Fa) ) {
55 tol1 = 2.*ZEPS * Abs(TheRoot) + 0.5 * Tolerance; // convergence check
56 xm = 0.5 * ( c - TheRoot );
57 if (Abs(xm) <= tol1 || TheError == 0. ) {
61 if (Abs(e) >= tol1 && Abs(Fa) > Abs(TheError) ) {
62 s = TheError / Fa; // attempt inverse quadratic interpolation
70 p = s * (2.*xm *q * (q - r) - (TheRoot - a)*(r - 1.));
71 q = (q -1.) * (r - 1.) * (s - 1.);
73 if ( p > 0. ) { q = -q;} // check whether in bounds
75 min1 = 3.* xm* q - Abs(tol1 *q);
77 if (2.* p < (min1 < min2 ? min1 : min2) ) {
78 e = d ; // accept interpolation
82 d = xm; // interpolation failed,use bissection
86 else { // bounds decreasing too slowly ,use bissection
90 a = TheRoot ; // move last best guess to a
92 if (Abs(d) > tol1) { // evaluate new trial root
96 TheRoot += (xm > 0. ? Abs(tol1) : -Abs(tol1));
98 F.Value(TheRoot,TheError);
100 Done = Standard_False;
105 void math_BracketedRoot::Dump(Standard_OStream& o) const {
107 o << "math_BracketedRoot ";
109 o << " Status = Done \n";
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;
115 o << " Status = not Done \n";