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