src/math/math_BracketedRoot.cxx

   1 // Copyright (c) 1997-1999 Matra Datavision
   2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   3 //
   4 // This file is part of Open CASCADE Technology software library.
   5 //
   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.
  11 //
  12 // Alternatively, this file may be used under the terms of Open CASCADE
  13 // commercial license or contractual agreement.
  14
  15
  16 #include <math_BracketedRoot.hxx>
  17 #include <math_Function.hxx>
  18 #include <StdFail_NotDone.hxx>
  19
  20 // reference algorithme:
  21 //                   Brent method
  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 ) {
  29
  30     Standard_Real Fa,Fc,a,c=0,d=0,e=0;
  31     Standard_Real min1,min2,p,q,r,s,tol1,xm;
  32
  33     a = Bound1;
  34     TheRoot = Bound2;
  35     F.Value(a,Fa);
  36     F.Value(TheRoot,TheError);
  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           }
  98           F.Value(TheRoot,TheError);
  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";
 110          o << " Number of iterations = " << NbIter << endl;
 111          o << " The Root is: " << TheRoot << endl;
 112          o << " The value at the root is: " << TheError << endl;
 113        }
 114        else {
 115          o << " Status = not Done \n";
 116        }
 117 }