Skip to content

Commit

Permalink
Add 4th order solver (Householder method) with three analytic derivat…
Browse files Browse the repository at this point in the history 
…ives ; closes #1185
  • Loading branch information
@ibell
ibell committed Aug 13, 2016
1 parent 1e58037 commit e7fc3eb
Show file tree
Hide file tree
Showing 3 changed files with 91 additions and 12 deletions.
19 changes: 15 additions & 4 deletions include/Solvers.h
View file
Expand Up @@ -35,6 +35,12 @@ class FuncWrapper1DWithTwoDerivs : public FuncWrapper1DWithDeriv
public:
virtual double second_deriv(double) = 0;
};

class FuncWrapper1DWithThreeDerivs : public FuncWrapper1DWithTwoDerivs
{
public:
virtual double third_deriv(double) = 0;
};

class FuncWrapperND
{
Expand All @@ -52,9 +58,10 @@ double Brent(FuncWrapper1D* f, double a, double b, double macheps, double t, int
double Secant(FuncWrapper1D* f, double x0, double dx, double ftol, int maxiter);
double BoundedSecant(FuncWrapper1D* f, double x0, double xmin, double xmax, double dx, double ftol, int maxiter);
double Newton(FuncWrapper1DWithDeriv* f, double x0, double ftol, int maxiter);
double Halley(FuncWrapper1DWithTwoDerivs* f, double x0, double ftol, int maxiter);
double Halley(FuncWrapper1DWithTwoDerivs* f, double x0, double ftol, int maxiter, double xtol_rel = 1e-12);
double Householder4(FuncWrapper1DWithThreeDerivs* f, double x0, double ftol, int maxiter, double xtol_rel = 1e-12);

// Single-Dimensional solvers
// Single-Dimensional solvers, refere
inline double Brent(FuncWrapper1D &f, double a, double b, double macheps, double t, int maxiter){
return Brent(&f, a, b, macheps, t, maxiter);
}
Expand All @@ -67,9 +74,13 @@ inline double BoundedSecant(FuncWrapper1D &f, double x0, double xmin, double xma
inline double Newton(FuncWrapper1DWithDeriv &f, double x0, double ftol, int maxiter){
return Newton(&f, x0, ftol, maxiter);
}
inline double Halley(FuncWrapper1DWithTwoDerivs &f, double x0, double ftol, int maxiter){
return Halley(&f, x0, ftol, maxiter);
inline double Halley(FuncWrapper1DWithTwoDerivs &f, double x0, double ftol, int maxiter, double xtol_rel = 1e-12){
return Halley(&f, x0, ftol, maxiter, xtol_rel);
}
inline double Householder4(FuncWrapper1DWithThreeDerivs &f, double x0, double ftol, int maxiter, double xtol_rel = 1e-12){
return Householder4(&f, x0, ftol, maxiter, xtol_rel);
}


// Multi-Dimensional solvers
std::vector<double> NDNewtonRaphson_Jacobian(FuncWrapperND *f, std::vector<double> &x0, double tol, int maxiter);
Expand Down
16 changes: 10 additions & 6 deletions src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp
View file
Expand Up @@ -2083,7 +2083,7 @@ CoolPropDbl HelmholtzEOSMixtureBackend::solver_rho_Tp(CoolPropDbl T, CoolPropDbl
phases phase;

// Define the residual to be driven to zero
class solver_TP_resid : public FuncWrapper1DWithTwoDerivs
class solver_TP_resid : public FuncWrapper1DWithThreeDerivs
{
public:
HelmholtzEOSMixtureBackend *HEOS;
Expand All @@ -2104,7 +2104,11 @@ CoolPropDbl HelmholtzEOSMixtureBackend::solver_rho_Tp(CoolPropDbl T, CoolPropDbl
};
double second_deriv(double rhomolar){
// d2p/drho2|T / pspecified
return R_u*T/rhomolar*(2*delta*HEOS->dalphar_dDelta() + 4*pow(delta, 2)*HEOS->d2alphar_dDelta2() + pow(delta, 3)*HEOS->calc_d3alphar_dDelta3())/p;
return R_u*T/rhor*(2*HEOS->dalphar_dDelta() + 4*delta*HEOS->d2alphar_dDelta2() + pow(delta, 2)*HEOS->calc_d3alphar_dDelta3())/p;
};
double third_deriv(double rhomolar){
// d3p/drho3|T / pspecified
return R_u*T/POW2(rhor)*(6*HEOS->d2alphar_dDelta2() + 4*delta*HEOS->d3alphar_dDelta3() + POW2(delta)*HEOS->calc_d4alphar_dDelta4())/p;
};
};
solver_TP_resid resid(this,T,p);
Expand Down Expand Up @@ -2150,8 +2154,8 @@ CoolPropDbl HelmholtzEOSMixtureBackend::solver_rho_Tp(CoolPropDbl T, CoolPropDbl
}
}
else{
// Try with Halley's method starting at a very high density
rhomolar = Halley(resid, 3*rhomolar_reducing(), 1e-8, 100);
// Try with 4th order Householder method starting at a very high density
rhomolar = Householder4(&resid, 3*rhomolar_reducing(), 1e-8, 100);
}
return rhomolar;
}
Expand All @@ -2168,8 +2172,8 @@ CoolPropDbl HelmholtzEOSMixtureBackend::solver_rho_Tp(CoolPropDbl T, CoolPropDbl

try{

// First we try with Halley's method with analytic derivative
double rhomolar = Halley(resid, rhomolar_guess, 1e-8, 100);
// First we try with 4th order Householder method with analytic derivatives
double rhomolar = Householder4(resid, rhomolar_guess, 1e-8, 100);
if (!ValidNumber(rhomolar)){
throw ValueError();
}
Expand Down
68 changes: 66 additions & 2 deletions src/Solvers.cpp
View file
Expand Up @@ -142,9 +142,10 @@ x_{n+1} = x_n - \frac {2 f(x_n) f'(x_n)} {2 {[f'(x_n)]}^2 - f(x_n) f''(x_n)}
@param x0 The inital guess for the solution
@param ftol The absolute value of the tolerance accepted for the objective function
@param maxiter Maximum number of iterations
@param xtol_rel The minimum allowable (relative) step size
@returns If no errors are found, the solution, otherwise the value _HUGE, the value for infinity
*/
double Halley(FuncWrapper1DWithTwoDerivs* f, double x0, double ftol, int maxiter)
double Halley(FuncWrapper1DWithTwoDerivs* f, double x0, double ftol, int maxiter, double xtol_rel)
{
double x, dx, fval=999, dfdx, d2fdx2;
int iter=1;
Expand All @@ -171,7 +172,7 @@ double Halley(FuncWrapper1DWithTwoDerivs* f, double x0, double ftol, int maxiter

x += dx;

if (std::abs(dx/x) < 10*DBL_EPSILON){
if (std::abs(dx/x) < xtol_rel){
return x;
}

Expand All @@ -183,6 +184,69 @@ double Halley(FuncWrapper1DWithTwoDerivs* f, double x0, double ftol, int maxiter
}
return x;
}

/**
In the 4-th order Householder method, three derivatives of the input variable are needed, it yields the following method:
\f[
x_{n+1} = x_n - f(x_n)\left( \frac {[f'(x_n)]^2 - f(x_n)f''(x_n)/2 } {[f'(x_n)]^3-f(x_n)f'(x_n)f''(x_n)+f'''(x_n)*[f(x_n)]^2/6 } \right)
\f]
http://numbers.computation.free.fr/Constants/Algorithms/newton.ps
@param f A pointer to an instance of the FuncWrapper1DWithThreeDerivs class that implements the call() and three derivatives
@param x0 The inital guess for the solution
@param ftol The absolute value of the tolerance accepted for the objective function
@param maxiter Maximum number of iterations
@param xtol_rel The minimum allowable (relative) step size
@returns If no errors are found, the solution, otherwise the value _HUGE, the value for infinity
*/
double Householder4(FuncWrapper1DWithThreeDerivs* f, double x0, double ftol, int maxiter, double xtol_rel)
{
double x, dx, fval=999, dfdx, d2fdx2, d3fdx3;
int iter=1;
f->errstring.clear();
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);
d3fdx3 = f->third_deriv(x);

if (!ValidNumber(fval)){
throw ValueError("Residual function in Householder4 returned invalid number");
};
if (!ValidNumber(dfdx)){
throw ValueError("Derivative function in Householder4 returned invalid number");
};
if (!ValidNumber(d2fdx2)){
throw ValueError("Second derivative function in Householder4 returned invalid number");
};
if (!ValidNumber(d3fdx3)){
throw ValueError("Third derivative function in Householder4 returned invalid number");
};

dx = -fval*(POW2(dfdx)-fval*d2fdx2/2.0)/(POW3(dfdx)-fval*dfdx*d2fdx2+d3fdx3*POW2(fval)/6.0);

x += dx;

if (std::abs(dx/x) < xtol_rel){
return x;
}

if (iter>maxiter){
f->errstring= "reached maximum number of iterations";
throw SolutionError(format("Householder4 reached maximum number of iterations"));
}
iter=iter+1;
}
return x;
}

/**
In the secant function, a 1-D Newton-Raphson solver is implemented. An initial guess for the solution is provided.
Expand Down

0 comments on commit e7fc3eb

Please sign in to comment.