Secant behaves properly now when y1 approx. y2; closes #1015

CoolProp · Apr 7, 2016 · c3f2867 · c3f2867
1 parent 2abcb60
commit c3f2867
Showing 1 changed file with 12 additions and 5 deletions.
diff --git a/src/Solvers.cpp b/src/Solvers.cpp
@@ -152,13 +152,14 @@ double Halley(FuncWrapper1DWithTwoDerivs* f, double x0, double ftol, int maxiter
     x = x0;
     while (iter < 2 || std::abs(fval) > ftol)
     {
+        if (f->input_not_in_range(x)){
+            throw ValueError(format("Input [%g] is out of range",x));
+        }
+
         fval = f->call(x);
         dfdx = f->deriv(x);
         d2fdx2 = f->second_deriv(x);
 
-        if (f->input_not_in_range(x)){
-            throw ValueError(format("Input [%g] is out of range",x));
-        }
         if (!ValidNumber(fval)){
             throw ValueError("Residual function in Halley returned invalid number");
         };
@@ -168,8 +169,6 @@ double Halley(FuncWrapper1DWithTwoDerivs* f, double x0, double ftol, int maxiter
 
         dx = -(2*fval*dfdx)/(2*POW2(dfdx)-fval*d2fdx2);
 
-
-
         x += dx;
 
         if (std::abs(dx/x) < 10*DBL_EPSILON){
@@ -212,6 +211,10 @@ double Secant(FuncWrapper1D* f, double x0, double dx, double tol, int maxiter)
         if (iter==1){x1=x0; x=x1;}
         if (iter==2){x2=x0+dx; x=x2;}
         if (iter>2) {x=x2;}
+
+            if (f->input_not_in_range(x)){
+                throw ValueError(format("Input [%g] is out of range",x));
+            }
 
             fval = f->call(x);
 
@@ -231,6 +234,10 @@ double Secant(FuncWrapper1D* f, double x0, double dx, double tol, int maxiter)
                 return x;
             }
             y2=fval;
+            double deltay = y2-y1;
+            if (iter > 2 && std::abs(deltay)<1e-14){
+                return x;
+            }
             x3=x2-y2/(y2-y1)*(x2-x1);
             y1=y2; x1=x2; x2=x3;