src/math/math_BissecNewton.cxx

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