src/math/math_BissecNewton.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_BissecNewton.hxx>
  17 #include <math_FunctionWithDerivative.hxx>
  18 #include <StdFail_NotDone.hxx>
  19
  20 //=======================================================================
  21 //function : math_BissecNewton
  22 //purpose  : Constructor
  23 //=======================================================================
  24 math_BissecNewton::math_BissecNewton(const Standard_Real theXTolerance)
  25 : TheStatus(math_NotBracketed),
  26   XTol     (theXTolerance),
  27   x        (0.0),
  28   dx       (0.0),
  29   f        (0.0),
  30   df       (0.0),
  31   Done     (Standard_False)
  32 {
  33 }
  34
  35 //=======================================================================
  36 //function : ~math_BissecNewton
  37 //purpose  : Destructor
  38 //=======================================================================
  39 math_BissecNewton::~math_BissecNewton()
  40 {
  41 }
  42
  43 //=======================================================================
  44 //function : Perform
  45 //purpose  :
  46 //=======================================================================
  47 void math_BissecNewton::Perform(math_FunctionWithDerivative& F,
  48                                 const Standard_Real    Bound1,
  49                                 const Standard_Real    Bound2,
  50                                 const Standard_Integer NbIterations)
  51 {
  52   Standard_Boolean GOOD;
  53   Standard_Integer j;
  54   Standard_Real dxold, fh, fl;
  55   Standard_Real swap, temp, xh, xl;
  56
  57   GOOD = F.Values(Bound1, fl, df);
  58   if(!GOOD) {
  59     Done = Standard_False;
  60     TheStatus = math_FunctionError;
  61     return;
  62   }
  63   GOOD = F.Values(Bound2, fh, df);
  64   if(!GOOD) {
  65     Done = Standard_False;
  66     TheStatus = math_FunctionError;
  67     return;
  68   }
  69 //  Modified by Sergey KHROMOV - Wed Jan 22 12:06:45 2003 Begin
  70   Standard_Real aFTol = RealEpsilon();
  71
  72 //   if(fl * fh >= 0.0) {
  73   if(fl * fh > aFTol*aFTol) {
  74     Done = Standard_False;
  75     TheStatus = math_NotBracketed;
  76     return;
  77   }
  78 //   if(fl < 0.0) {
  79   if(fl < -aFTol || (fl < aFTol && fh < -aFTol)) {
  80     xl = Bound1;
  81     xh = Bound2;
  82   }
  83   else {
  84     xl = Bound2;
  85     xh = Bound1;
  86     swap = fl;
  87     fl = fh;
  88     fh = swap;
  89   }
  90 //  Modified by Sergey KHROMOV - Wed Jan 22 12:06:49 2003 End
  91   x = 0.5 * (Bound1 + Bound2);
  92   dxold = fabs(Bound2 - Bound1);
  93   dx = dxold;
  94   GOOD = F.Values(x, f, df);
  95   if(!GOOD) {
  96     Done = Standard_False;
  97     TheStatus = math_FunctionError;
  98     return;
  99   }
 100   for(j = 1; j <= NbIterations; j++) {
 101     if((((x - xh) * df - f) * ((x - xl) * df - f) >= 0.0)
 102        || (fabs(2.0 * f) > fabs(dxold * df))) {
 103       dxold = dx;
 104       dx = 0.5 * (xh - xl);
 105       x = xl + dx;
 106       if(Abs(dx) < XTol) {
 107         TheStatus = math_OK;
 108         Done = Standard_True;
 109         return;
 110       }
 111     }
 112     else {
 113       dxold = dx;
 114       dx = f / df;
 115       temp = x;
 116       x -= dx;
 117       if(temp == x) {
 118         TheStatus = math_OK;
 119         Done = Standard_True;
 120         return;
 121       }
 122     }
 123     if(IsSolutionReached(F)) {
 124       TheStatus = math_OK;
 125       Done = Standard_True;
 126       return;
 127     }
 128     GOOD = F.Values(x, f, df);
 129     if(!GOOD) {
 130       Done = Standard_False;
 131       TheStatus = math_FunctionError;
 132       return;
 133     }
 134     if(f < 0.0) {
 135       xl = x;
 136       fl = f;
 137     }
 138     else if(f > 0.0) {
 139       xh = x;
 140       fh = f;
 141     }
 142     else {
 143       TheStatus = math_OK;
 144       Done = Standard_True;
 145       return;
 146     }
 147   }
 148   TheStatus = math_TooManyIterations;
 149   Done = Standard_False;
 150   return;
 151 }
 152
 153 //=======================================================================
 154 //function : Dump
 155 //purpose  :
 156 //=======================================================================
 157 void math_BissecNewton::Dump(Standard_OStream& o) const {
 158
 159   o << "math_BissecNewton ";
 160   if(Done) {
 161     o << " Status = Done \n";
 162     o << " The Root  is: " << x << std::endl;
 163     o << " The value at this Root is: " << f << std::endl;
 164   }
 165   else {
 166     o << " Status = not Done \n";
 167   }
 168 }
 169