From 1ccceeb7a595f1e1710996e0cbe64d24cb47a055 Mon Sep 17 00:00:00 2001 From: jowr Date: Wed, 11 Jun 2014 19:27:07 +0200 Subject: [PATCH] Added a bounded Brent solver to Polynomial2D --- include/PolyMath.h | 53 ++++++++++- src/PolyMath.cpp | 226 ++++++++++++++++++++++++++++++--------------- 2 files changed, 202 insertions(+), 77 deletions(-) diff --git a/include/PolyMath.h b/include/PolyMath.h index 99047676..64567ec3 100644 --- a/include/PolyMath.h +++ b/include/PolyMath.h @@ -222,16 +222,43 @@ public: /// @param axis unsigned integer value that represents the axis to solve for (0=x, 1=y) double solve(const double &in, const double &z_in, int axis); - /// @param y_in double value that represents the current input in x (1st dimension) or y (2nd dimension) + /// @param y_in double value that represents the current input in y (2nd dimension) /// @param z_in double value that represents the current output in the 3rd dimension double solve_x(const double &y_in, const double &z_in){return solve(y_in, z_in, 0);} - /// @param x_in double value that represents the current input in x (1st dimension) or y (2nd dimension) + /// @param x_in double value that represents the current input in x (1st dimension) /// @param z_in double value that represents the current output in the 3rd dimension - /// @param axis unsigned integer value that represents the axis to solve for (0=x, 1=y) double solve_y(const double &x_in, const double &z_in){return solve(x_in, z_in, 1);} + /// @param in double value that represents the current input in x (1st dimension) or y (2nd dimension) + /// @param z_in double value that represents the current output in the 3rd dimension + /// @param axis unsigned integer value that represents the axis to solve for (0=x, 1=y) + double solve_limits(const double &in, const double &z_in, const double &min, const double &max, const int &axis); + + /// @param y_in double value that represents the current input in y (2nd dimension) + /// @param z_in double value that represents the current output in the 3rd dimension + double solve_limits_x(const double &y_in, const double &z_in, const double &x_min, const double &x_max){return solve_limits(y_in, z_in, x_min, x_max, 0);} + + /// @param x_in double value that represents the current input in x (1st dimension) + /// @param z_in double value that represents the current output in the 3rd dimension + double solve_limits_y(const double &x_in, const double &z_in, const double &y_min, const double &y_max){return solve_limits(x_in, z_in, y_min, y_max, 1);} + + + protected: +// /// @param x_in double value that represents the current input in x (1st dimension) +// /// @param y_in double value that represents the current input in y (2nd dimension) +// /// @param z_in double value that represents the current output in the 3rd dimension +// double residual(const double &x_in, const double &y_in, const double &z_in){return this->evaluate(x_in,y_in)-z_in;} +// +// /// @param x_in double value that represents the current input in x (1st dimension) +// /// @param z_in double value that represents the current output in the 3rd dimension +// double residual_x(const double &x_in, const double &z_in){return residual(x_in, this->y_std, z_in);} +// +// /// @param y_in double value that represents the current input in y (2nd dimension) +// /// @param z_in double value that represents the current output in the 3rd dimension +// double residual_y(const double &y_in, const double &z_in){return residual(this->x_std, y_in, z_in);} + /// Simple polynomial function generator. <- Deprecated due to poor performance, use Horner-scheme instead /** Base function to produce n-th order polynomials * based on the length of the coefficient vector. @@ -252,7 +279,27 @@ protected: }; +class Poly2DResidual : public FuncWrapper1D { +protected: + enum dims {iX, iY}; + int targetDim; + /// the fixed input != targetDim + double fixed; + /// Object that evaluates the equation + Polynomial2D poly; + /// Current output value + double output; +private: + Poly2DResidual(); + +public: + Poly2DResidual(const Polynomial2D &poly, const double &fixed, const int &targetDim, const double &output); + virtual ~Poly2DResidual(){}; + + double call(double in); + double deriv(double in); +}; diff --git a/src/PolyMath.cpp b/src/PolyMath.cpp index 00afc844..370cb87f 100644 --- a/src/PolyMath.cpp +++ b/src/PolyMath.cpp @@ -18,17 +18,6 @@ namespace CoolProp{ -/// Set the coefficient matrix. -/// @param coefficients matrix containing the ordered coefficients -void Polynomial2D::setCoefficients(const Eigen::MatrixXd &coefficients){ - this->coefficients = coefficients; - this->coefficientsDerX = this->deriveCoeffs(coefficients,0); - this->coefficientsDerY = this->deriveCoeffs(coefficients,1); -} -void Polynomial2D::setCoefficients(const std::vector > &coefficients){ - this->setCoefficients(vec_to_eigen(coefficients)); -} - /// Basic checks for coefficient vectors. /** Starts with only the first coefficient dimension * and checks the matrix size against the parameters rows and columns. */ @@ -87,27 +76,27 @@ Eigen::MatrixXd Polynomial2D::integrateCoeffs(const Eigen::MatrixXd &coefficient std::size_t r = coefficients.rows(), c = coefficients.cols(); Eigen::MatrixXd newCoefficients; switch (axis) { - case 0: - newCoefficients = Eigen::MatrixXd::Zero(r+1,c); - newCoefficients.block(1,0,r,c) = coefficients.block(0,0,r,c); - for (size_t i = 0; i < r; ++i) { - for (size_t j = 0; j < c; ++j) { - newCoefficients(i+1,j) /= (i+1.); - } + case 0: + newCoefficients = Eigen::MatrixXd::Zero(r+1,c); + newCoefficients.block(1,0,r,c) = coefficients.block(0,0,r,c); + for (size_t i = 0; i < r; ++i) { + for (size_t j = 0; j < c; ++j) { + newCoefficients(i+1,j) /= (i+1.); } - break; - case 1: - newCoefficients = Eigen::MatrixXd::Zero(r,c+1); - newCoefficients.block(0,1,r,c) = coefficients.block(0,0,r,c); - for (size_t i = 0; i < r; ++i) { - for (size_t j = 0; j < c; ++j) { - newCoefficients(i,j+1) /= (j+1.); - } + } + break; + case 1: + newCoefficients = Eigen::MatrixXd::Zero(r,c+1); + newCoefficients.block(0,1,r,c) = coefficients.block(0,0,r,c); + for (size_t i = 0; i < r; ++i) { + for (size_t j = 0; j < c; ++j) { + newCoefficients(i,j+1) /= (j+1.); } - break; - default: - throw ValueError(format("You have to provide a dimension, 0 or 1, for integration, %d is not valid. ",axis)); - break; + } + break; + default: + throw ValueError(format("You have to provide a dimension, 0 or 1, for integration, %d is not valid. ",axis)); + break; } return newCoefficients; } @@ -123,27 +112,25 @@ Eigen::MatrixXd Polynomial2D::deriveCoeffs(const Eigen::MatrixXd &coefficients, std::size_t r = coefficients.rows(), c = coefficients.cols(); Eigen::MatrixXd newCoefficients(coefficients); switch (axis) { - case 0: - //newCoefficients.resize(r-1,c); - for (size_t i = 1; i < r; ++i) { - for (size_t j = 0; j < c; ++j) { - newCoefficients(i,j) *= i; - } + case 0: + for (size_t i = 1; i < r; ++i) { + for (size_t j = 0; j < c; ++j) { + newCoefficients(i,j) *= i; } - removeRow(newCoefficients,0); - break; - case 1: - //newCoefficients.resize(r,c-1); - for (size_t i = 0; i < r; ++i) { - for (size_t j = 1; j < c; ++j) { - newCoefficients(i,j) *= j; - } + } + removeRow(newCoefficients,0); + break; + case 1: + for (size_t i = 0; i < r; ++i) { + for (size_t j = 1; j < c; ++j) { + newCoefficients(i,j) *= j; } - removeColumn(newCoefficients,0); - break; - default: - throw ValueError(format("You have to provide a dimension, 0 or 1, for derivation, %d is not valid. ",axis)); - break; + } + removeColumn(newCoefficients,0); + break; + default: + throw ValueError(format("You have to provide a dimension, 0 or 1, for derivation, %d is not valid. ",axis)); + break; } return newCoefficients; } @@ -187,15 +174,15 @@ double Polynomial2D::evaluate(const Eigen::MatrixXd &coefficients, const double double Polynomial2D::derivative(const double &x_in, const double &y_in, int axis = -1){ double result = 0; switch (axis) { - case 0: - result = this->evaluate(this->coefficientsDerX, x_in,y_in); - break; - case 1: - result = this->evaluate(this->coefficientsDerY, x_in,y_in); - break; - default: - throw ValueError(format("You have to provide a dimension, 0 or 1, for derivation, %d is not valid. ",axis)); - break; + case 0: + result = this->evaluate(this->coefficientsDerX, x_in,y_in); + break; + case 1: + result = this->evaluate(this->coefficientsDerY, x_in,y_in); + break; + default: + throw ValueError(format("You have to provide a dimension, 0 or 1, for derivation, %d is not valid. ",axis)); + break; } return result; } @@ -206,21 +193,21 @@ double Polynomial2D::solve(const double &in, const double &z_in, int axis = -1){ std::size_t r = coefficients.rows(), c = coefficients.cols(); Eigen::VectorXd tmpCoefficients; switch (axis) { - case 0: - tmpCoefficients = Eigen::VectorXd::Zero(r); - for(size_t i=0; ido_debug()) std::cout << "Coefficients: " << mat_to_string(Eigen::MatrixXd(tmpCoefficients)) << std::endl; @@ -231,6 +218,45 @@ double Polynomial2D::solve(const double &in, const double &z_in, int axis = -1){ return rootsVec[0]; // TODO: implement root selection algorithm } +/// @param in double value that represents the current input in x (1st dimension) or y (2nd dimension) +/// @param z_in double value that represents the current output in the 3rd dimension +/// @param axis unsigned integer value that represents the axis to solve for (0=x, 1=y) +double Polynomial2D::solve_limits(const double &in, const double &z_in, const double &min, const double &max, const int &axis){ + + Poly2DResidual res = Poly2DResidual(*this, in, axis, z_in); + std::string errstring; + + double macheps = DBL_EPSILON; + double tol = DBL_EPSILON*1e3; + int maxiter = 50; + + return Brent(res, min, max, macheps, tol, maxiter, errstring); +} + + +// std::string errstring; +// double result = -1.0; +// switch (this->uses) { +// case iNewton: ///< Newton solver with derivative and guess value +// if (res.is2D()) { +// throw CoolProp::NotImplementedError("The Newton solver is not suitable for 2D polynomials, yet."); +// } +// result = Newton(res, this->guess, this->tol, this->maxiter, errstring); +// break; +// +// case iBrent: ///< Brent solver with bounds +// result = Brent(res, this->min, this->max, this->macheps, this->tol, this->maxiter, errstring); +// break; +// +// default: +// throw CoolProp::NotImplementedError("This solver has not been implemented or you forgot to select a solver..."); +// } +// return result; +//} + + + + /// Simple polynomial function generator. <- Deprecated due to poor performance, use Horner-scheme instead /** Base function to produce n-th order polynomials @@ -279,6 +305,49 @@ double Polynomial2D::baseHorner(std::vector< std::vector > const& coeffi return result; } + + +Poly2DResidual::Poly2DResidual(const Polynomial2D &poly, const double &fixed, const int &targetDim, const double &output){ + this->poly = poly; + this->fixed = fixed; + + switch (targetDim) { + case iX: + case iY: + this->targetDim = targetDim; + break; + default: + throw ValueError(format("You have to provide a dimension to the solver, %d is not valid. ",targetDim)); + break; + } + + this->output = output; +} + +double Poly2DResidual::call(double in){ + if (targetDim==iX) return poly.evaluate(in,fixed)-output; + if (targetDim==iY) return poly.evaluate(fixed,in)-output; + return -_HUGE; +} + +double Poly2DResidual::deriv(double in){ + if (targetDim==iX) return poly.derivative(in,fixed,0)-output; + if (targetDim==iY) return poly.derivative(fixed,in,1)-output; + return -_HUGE; +} + + +/// Set the coefficient matrix. +/// @param coefficients matrix containing the ordered coefficients +void Polynomial2D::setCoefficients(const Eigen::MatrixXd &coefficients){ + this->coefficients = coefficients; + this->coefficientsDerX = this->deriveCoeffs(coefficients,0); + this->coefficientsDerY = this->deriveCoeffs(coefficients,1); +} +void Polynomial2D::setCoefficients(const std::vector > &coefficients){ + this->setCoefficients(vec_to_eigen(coefficients)); +} + } @@ -402,7 +471,7 @@ TEST_CASE("Internal consistency checks and example use cases for PolyMath.cpp"," CHECK( check_abs(c,d,acc) ); } - //c = 2.0; + c = 2.0; c = poly.solve_x(0.0, d); { CAPTURE(T); @@ -411,6 +480,15 @@ TEST_CASE("Internal consistency checks and example use cases for PolyMath.cpp"," CHECK( check_abs(c,T,acc) ); } + c = 2.0; + c = poly.solve_limits_x(0.0, d, -50, 750); + { + CAPTURE(T); + CAPTURE(c); + CAPTURE(d); + CHECK( check_abs(c,T,acc) ); + } + // T = 0.0; // solve.setGuess(75+273.15); // T = solve.polyval(cHeat,c);