From 4e2aa1c4041a4d6bd9d0bba77a141cb4ab0ba5dc Mon Sep 17 00:00:00 2001 From: jowr Date: Fri, 13 Jun 2014 13:56:08 +0200 Subject: [PATCH] Started to work on flexible polynomial solvers for fractional exponents --- include/CoolPropTools.h | 4 +- include/MatrixMath.h | 16 +- include/PolyMath.h | 360 +++++++++++++++----------------- src/PolyMath.cpp | 447 +++++++++++++++++++++++++--------------- 4 files changed, 456 insertions(+), 371 deletions(-) diff --git a/include/CoolPropTools.h b/include/CoolPropTools.h index 39b5e57d..fda4631e 100644 --- a/include/CoolPropTools.h +++ b/include/CoolPropTools.h @@ -243,8 +243,8 @@ max = min; min = fabs(A); } - if (max>D) return ( ( 1.0-min/max*1e0 ) < D ); - else throw CoolProp::ValueError(format("Too small numbers: %f cannot be tested with an accepted error of %f. ",max,D)); + if (max>DBL_EPSILON*1e3) return ( ( 1.0-min/max*1e0 ) < D ); + else throw CoolProp::ValueError(format("Too small numbers: %f cannot be tested with an accepted error of %f for a machine precision of %f. ",max,D,DBL_EPSILON)); }; bool inline check_abs(double A, double B){ return check_abs(A,B,1e5*DBL_EPSILON); diff --git a/include/MatrixMath.h b/include/MatrixMath.h index cf2910ac..f0d46c58 100644 --- a/include/MatrixMath.h +++ b/include/MatrixMath.h @@ -35,16 +35,16 @@ namespace CoolProp{ * it is still needed. */ /// @param coefficients matrix containing the ordered coefficients -template std::vector eigen_to_vec1D(const Eigen::Matrix &coefficients, int axis = 1){ +template std::vector eigen_to_vec1D(const Eigen::Matrix &coefficients, int axis = 0){ std::vector result; size_t r = coefficients.rows(), c = coefficients.cols(); - if (axis==1) { + if (axis==0) { if (c!=1) throw ValueError(format("Your matrix has the wrong dimensions: %d,%d",r,c)); result.resize(r); for (size_t i = 0; i < r; ++i) { result[i] = coefficients(i,0); } - } else if (axis==2) { + } else if (axis==1) { if (r!=1) throw ValueError(format("Your matrix has the wrong dimensions: %d,%d",r,c)); result.resize(c); for (size_t i = 0; i < c; ++i) { @@ -84,15 +84,15 @@ template Eigen::Matrix vec_to_eigen(c return result; } /// @param coefficients matrix containing the ordered coefficients -template Eigen::Matrix vec_to_eigen(const std::vector &coefficients, int axis = 1){ +template Eigen::Matrix vec_to_eigen(const std::vector &coefficients, int axis = 0){ size_t nRows = num_rows(coefficients); Eigen::Matrix result; - if (axis==1) result.resize(nRows,1); - else if (axis==2) result.resize(1,nRows); + if (axis==0) result.resize(nRows,1); + else if (axis==1) result.resize(1,nRows); else throw ValueError(format("You have to provide axis information: %d is not valid. ",axis)); for (size_t i = 0; i < nRows; ++i) { - if (axis==1) result(i,0) = coefficients[i]; - if (axis==2) result(0,i) = coefficients[i]; + if (axis==0) result(i,0) = coefficients[i]; + if (axis==1) result(0,i) = coefficients[i]; } return result; } diff --git a/include/PolyMath.h b/include/PolyMath.h index 833b4736..b2822321 100644 --- a/include/PolyMath.h +++ b/include/PolyMath.h @@ -11,135 +11,32 @@ //#include // inner_product //#include //#include "float.h" - +#include "MatrixMath.h" #include namespace CoolProp{ - -///// The base class for all Polynomials -//class Polynomial1D{ -// -//private: -// Eigen::VectorXd coefficients; -// -//public: -// /// Constructors -// Polynomial1D(); -// Polynomial1D(const Eigen::VectorXd &coefficients){ -// this->setCoefficients(coefficients); -// } -// Polynomial1D(const std::vector &coefficients){ -// this->setCoefficients(coefficients); -// } -// -// /// Destructor. No implementation -// virtual ~Polynomial1D(){}; -// -//public: -// /// Set the coefficient vector. -// /// @param coefficients vector containing the ordered coefficients -// bool setCoefficients(const Eigen::VectorXd &coefficients); -// /// @param coefficients vector containing the ordered coefficients -// bool setCoefficients(const std::vector &coefficients); -// -// /// Basic checks for coefficient vectors. -// /** Starts with only the first coefficient dimension -// * and checks the vector length against parameter n. */ -// /// @param n unsigned integer value that represents the desired degree of the polynomial -// bool checkCoefficients(const unsigned int n); -// /// @param coefficients vector containing the ordered coefficients -// /// @param n unsigned integer value that represents the desired degree of the polynomial -// bool checkCoefficients(const Eigen::VectorXd &coefficients, const unsigned int n); -// /// @param coefficients vector containing the ordered coefficients -// /// @param n unsigned integer value that represents the desired degree of the polynomial -// bool checkCoefficients(const std::vector &coefficients, const unsigned int n); -// -//protected: -// /// Integration functions -// /** Integrating coefficients for polynomials is done by dividing the -// * original coefficients by (i+1) and elevating the order by 1 -// * through adding a zero as first coefficient. -// * Some reslicing needs to be applied to integrate along the x-axis. -// * In the brine/solution equations, reordering of the parameters -// * avoids this expensive operation. However, it is included for the -// * sake of completeness. -// */ -// Eigen::VectorXd integrateCoeffs(const Eigen::VectorXd &coefficients); -// std::vector integrateCoeffs(const std::vector &coefficients); -// -// /// Derivative coefficients calculation -// /** Deriving coefficients for polynomials is done by multiplying the -// * original coefficients with i and lowering the order by 1. -// * -// * It is not really deprecated, but untested and therefore a warning -// * is issued. Please check this method before you use it. -// */ -// Eigen::VectorXd deriveCoeffs(const Eigen::VectorXd &coefficients); -// std::vector deriveCoeffs(const std::vector &coefficients); -// -//public: -// /// The core functions to evaluate the polynomial -// /** It is here we implement the different special -// * functions that allow us to specify certain -// * types of polynomials. -// * The derivative might bee needed during the -// * solution process of the solver. It could also -// * be a protected function... -// */ -// double evaluate(const double &x_in); -// double derivative(const double &x_in); -// double solve(const double &y_in); -//}; - - /// The base class for all Polynomials class Polynomial2D { -private: - Eigen::MatrixXd coefficients; - Eigen::MatrixXd coefficientsDerX,coefficientsDerY; - double x_std; - double y_std; - public: /// Constructors - Polynomial2D(){x_std = 0.0; y_std = 0.0;}; - Polynomial2D(const Eigen::MatrixXd &coefficients){ - x_std = 0.0; y_std = 0.0; - this->setCoefficients(coefficients); - } - Polynomial2D(const std::vector > &coefficients){ - x_std = 0.0; y_std = 0.0; - this->setCoefficients(coefficients); - } + Polynomial2D(){}; /// Destructor. No implementation virtual ~Polynomial2D(){}; public: - /// Set the coefficient matrix. + /// Convert the coefficient vector. + /// @param coefficients vector containing the ordered coefficients + Eigen::MatrixXd convertCoefficients(const std::vector &coefficients){return vec_to_eigen(coefficients);} + /// Convert the coefficient matrix. /// @param coefficients matrix containing the ordered coefficients - void setCoefficients(const Eigen::MatrixXd &coefficients); - void setCoefficients(const std::vector > &coefficients); - - /// Set the standard inputs. - /// @param x_std fixed input for the first dimension - void set_x(const double &x_std){this->x_std = x_std;} - /// @param y_std fixed input for the second dimension - void set_y(const double &y_std){this->y_std = y_std;} + Eigen::MatrixXd convertCoefficients(const std::vector > &coefficients){return 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. */ - /// @param rows unsigned integer value that represents the desired degree of the polynomial - bool checkCoefficients(const unsigned int rows); - /// @param rows unsigned integer value that represents the desired degree of the polynomial in the 1st dimension - /// @param columns unsigned integer value that represents the desired degree of the polynomial in the 2nd dimension - bool checkCoefficients(const unsigned int rows, const unsigned int columns); - /// @param coefficients vector containing the ordered coefficients - /// @param rows unsigned integer value that represents the desired degree of the polynomial - bool checkCoefficients(const Eigen::MatrixXd &coefficients, const unsigned int rows); /// @param coefficients matrix containing the ordered coefficients /// @param rows unsigned integer value that represents the desired degree of the polynomial in the 1st dimension /// @param columns unsigned integer value that represents the desired degree of the polynomial in the 2nd dimension @@ -157,7 +54,8 @@ public: */ /// @param coefficients matrix containing the ordered coefficients /// @param axis unsigned integer value that represents the desired direction of integration - Eigen::MatrixXd integrateCoeffs(const Eigen::MatrixXd &coefficients, int axis); + /// @param times integer value that represents the desired order of integration + Eigen::MatrixXd integrateCoeffs(const Eigen::MatrixXd &coefficients, const int &axis, const int ×); /// Derivative coefficients calculation /** Deriving coefficients for polynomials is done by multiplying the @@ -165,20 +63,23 @@ public: */ /// @param coefficients matrix containing the ordered coefficients /// @param axis unsigned integer value that represents the desired direction of derivation - Eigen::MatrixXd deriveCoeffs(const Eigen::MatrixXd &coefficients, int axis); + /// @param times integer value that represents the desired order of derivation + Eigen::MatrixXd deriveCoeffs(const Eigen::MatrixXd &coefficients, const int &axis, const int ×); public: /// The core functions to evaluate the polynomial /** It is here we implement the different special * functions that allow us to specify certain * types of polynomials. - * The derivative might bee needed during the - * solution process of the solver. It could also - * be a protected function... + * + * Try to avoid many calls to the derivative and integral functions. + * Both of them have to calculate the new coefficients internally, + * which slows things down. Instead, you should use the deriveCoeffs + * and integrateCoeffs functions and store the coefficient matrix + * you need for future calls to evaluate derivative and integral. */ - /// @param coefficients vector containing the ordered coefficients - /// @param x_in double value that represents the current input + /// @param x_in double value that represents the current input in the 1st dimension double evaluate(const Eigen::MatrixXd &coefficients, const double &x_in); /// @param coefficients vector containing the ordered coefficients @@ -186,101 +87,44 @@ public: /// @param y_in double value that represents the current input in the 2nd dimension double evaluate(const Eigen::MatrixXd &coefficients, const double &x_in, const double &y_in); + /// @param coefficients vector containing the ordered coefficients /// @param x_in double value that represents the current input in the 1st dimension /// @param y_in double value that represents the current input in the 2nd dimension - double evaluate(const double &x_in, const double &y_in){return evaluate(this->coefficients, x_in, y_in);} - - /// @param x_in double value that represents the current input in the 1st dimension - double evaluate_x(const double &x_in){return evaluate(x_in, this->y_std);} - - /// @param y_in double value that represents the current input in the 2nd dimension - double evaluate_y(const double &y_in){return evaluate(this->x_std, y_in);} + /// @param axis unsigned integer value that represents the axis to derive for (0=x, 1=y) + double derivative(const Eigen::MatrixXd &coefficients, const double &x_in, const double &y_in, const int &axis); + /// @param coefficients vector containing the ordered coefficients /// @param x_in double value that represents the current input in the 1st dimension /// @param y_in double value that represents the current input in the 2nd dimension - /// @param axis unsigned integer value that represents the axis to solve for (0=x, 1=y) - double derivative(const double &x_in, const double &y_in, int axis); - - /// @param x_in double value that represents the current input in the 1st dimension - /// @param axis unsigned integer value that represents the axis to solve for (0=x, 1=y) - double derivative_x(const double &x_in, int axis = -1){return derivative(x_in, this->y_std, axis);} - - /// @param y_in double value that represents the current input in the 2nd dimension - /// @param axis unsigned integer value that represents the axis to solve for (0=x, 1=y) - double derivative_y(const double &y_in, int axis = -1){return derivative(this->x_std, y_in, axis);} - - /// @param x_in double value that represents the current input in the 1st dimension - /// @param y_in double value that represents the current input in the 2nd dimension - double dzdx(const double &x_in, const double &y_in){return derivative(x_in, y_in, 0);} - - /// @param y_in double value that represents the current input in the 2nd dimension - /// @param y_in double value that represents the current input in the 2nd dimension - double dzdy(const double &x_in, const double &y_in){return derivative(x_in, y_in, 1);} + /// @param axis unsigned integer value that represents the axis to integrate for (0=x, 1=y) + double integral(const Eigen::MatrixXd &coefficients, const double &x_in, const double &y_in, const int &axis); + /// Returns a vector with ALL the real roots of p(x_in,y_in)-z_in + /// @param coefficients vector containing the ordered coefficients /// @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(const double &in, const double &z_in, 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_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) - /// @param z_in double value that represents the current output in the 3rd dimension - double solve_y(const double &x_in, const double &z_in){return solve(x_in, z_in, 1);} + Eigen::VectorXd solve(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const int &axis); + /// Uses the Brent solver to find the roots of p(x_in,y_in)-z_in + /// @param coefficients vector containing the ordered coefficients /// @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 min double value that represents the minimum value /// @param max double value that represents the maximum value /// @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 - /// @param min double value that represents the minimum value in x (1st dimension) - /// @param max double value that represents the maximum value in x (1st 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 - /// @param min double value that represents the minimum value in y (2nd dimension) - /// @param max double value that represents the maximum value in y (2nd 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);} + double solve_limits(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const double &min, const double &max, const int &axis); + /// Uses the Newton solver to find the roots of p(x_in,y_in)-z_in + /// @param coefficients vector containing the ordered coefficients /// @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 guess double value that represents the start value /// @param axis unsigned integer value that represents the axis to solve for (0=x, 1=y) - double solve_guess(const double &in, const double &z_in, const double &guess, 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 - /// @param x_guess double value that represents the start value in x (1st dimension) - double solve_guess_x(const double &y_in, const double &z_in, const double &x_guess){return solve_guess(y_in, z_in, x_guess, 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 - /// @param y_guess double value that represents the start value in y (2nd dimension) - double solve_guess_y(const double &x_in, const double &z_in, const double &y_guess){return solve_guess(x_in, z_in, y_guess, 1);} - + double solve_guess(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const double &guess, const int &axis); 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. @@ -304,23 +148,149 @@ protected: class Poly2DResidual : public FuncWrapper1D { protected: enum dims {iX, iY}; - int targetDim; + Eigen::MatrixXd coefficients; + bool derIsSet; + Eigen::MatrixXd coefficientsDer; + int axis; /// the fixed input != targetDim - double fixed; + double in; /// Object that evaluates the equation Polynomial2D poly; /// Current output value - double output; + double z_in; private: Poly2DResidual(); public: - Poly2DResidual(const Polynomial2D &poly, const double &fixed, const int &targetDim, const double &output); + Poly2DResidual(Polynomial2D &poly, const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const int &axis); virtual ~Poly2DResidual(){}; - double call(double in); - double deriv(double in); + double call(double target); + double deriv(double target); +}; + + +/// A flexible class for polynomials starting at an arbitrary degree != 0, can be negative. +class Polynomial2Dflex : Polynomial2D{ + +public: + /// Constructors + Polynomial2Dflex(){}; + + /// Destructor. No implementation + virtual ~Polynomial2Dflex(){}; + +public: + /// Integration functions + /** Integrating coefficients for polynomials is done by dividing the + * original coefficients by (i+1) and elevating the order by 1 + * through adding a zero as first coefficient. + * Some reslicing needs to be applied to integrate along the x-axis. + * In the brine/solution equations, reordering of the parameters + * avoids this expensive operation. However, it is included for the + * sake of completeness. + */ + /// @param coefficients matrix containing the ordered coefficients + /// @param axis unsigned integer value that represents the desired direction of integration + /// @param times integer value that represents the desired order of integration + /// @param firstExponent integer value that represents the first exponent of the polynomial in axis direction + Eigen::MatrixXd integrateCoeffs(const Eigen::MatrixXd &coefficients, const int &axis, const int ×, const int &firstExponent); + + /// Derivative coefficients calculation + /** Deriving coefficients for polynomials is done by multiplying the + * original coefficients with i and lowering the order by 1. + */ + /// @param coefficients matrix containing the ordered coefficients + /// @param axis unsigned integer value that represents the desired direction of derivation + /// @param times integer value that represents the desired order of derivation + /// @param firstExponent integer value that represents the first exponent of the polynomial in axis direction + Eigen::MatrixXd deriveCoeffs(const Eigen::MatrixXd &coefficients, const int &axis, const int ×, const int &firstExponent); + +public: + /// The core functions to evaluate the polynomial + /** It is here we implement the different special + * functions that allow us to specify certain + * types of polynomials. + * + * Try to avoid many calls to the derivative and integral functions. + * Both of them have to calculate the new coefficients internally, + * which slows things down. Instead, you should use the deriveCoeffs + * and integrateCoeffs functions and store the coefficient matrix + * you need for future calls to evaluate derivative and integral. + */ + /// @param coefficients vector containing the ordered coefficients + /// @param x_in double value that represents the current input in the 1st dimension + /// @param firstExponent integer value that represents the first exponent of the polynomial + double evaluate(const Eigen::MatrixXd &coefficients, const double &x_in, const int &firstExponent); + + /// @param coefficients vector containing the ordered coefficients + /// @param x_in double value that represents the current input in the 1st dimension + /// @param y_in double value that represents the current input in the 2nd dimension + /// @param x_exp integer value that represents the first exponent of the polynomial in the 1st dimension + /// @param y_exp integer value that represents the first exponent of the polynomial in the 2nd dimension + double evaluate(const Eigen::MatrixXd &coefficients, const double &x_in, const double &y_in, const int &x_exp, const int &y_exp); + + /// @param coefficients vector containing the ordered coefficients + /// @param x_in double value that represents the current input in the 1st dimension + /// @param y_in double value that represents the current input in the 2nd dimension + /// @param axis unsigned integer value that represents the axis to derive for (0=x, 1=y) + /// @param x_exp integer value that represents the first exponent of the polynomial in the 1st dimension + /// @param y_exp integer value that represents the first exponent of the polynomial in the 2nd dimension + double derivative(const Eigen::MatrixXd &coefficients, const double &x_in, const double &y_in, const int &axis, const int &x_exp, const int &y_exp); + + /// @param coefficients vector containing the ordered coefficients + /// @param x_in double value that represents the current input in the 1st dimension + /// @param y_in double value that represents the current input in the 2nd dimension + /// @param axis unsigned integer value that represents the axis to integrate for (0=x, 1=y) + /// @param x_exp integer value that represents the first exponent of the polynomial in the 1st dimension + /// @param y_exp integer value that represents the first exponent of the polynomial in the 2nd dimension + double integral(const Eigen::MatrixXd &coefficients, const double &x_in, const double &y_in, const int &axis, const int &x_exp, const int &y_exp); + + /// Returns a vector with ALL the real roots of p(x_in,y_in)-z_in + /// @param coefficients vector containing the ordered coefficients + /// @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) + /// @param x_exp integer value that represents the first exponent of the polynomial in the 1st dimension + /// @param y_exp integer value that represents the first exponent of the polynomial in the 2nd dimension + Eigen::VectorXd solve(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const int &axis, const int &x_exp, const int &y_exp); + + /// Uses the Brent solver to find the roots of p(x_in,y_in)-z_in + /// @param coefficients vector containing the ordered coefficients + /// @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 min double value that represents the minimum value + /// @param max double value that represents the maximum value + /// @param axis unsigned integer value that represents the axis to solve for (0=x, 1=y) + /// @param x_exp integer value that represents the first exponent of the polynomial in the 1st dimension + /// @param y_exp integer value that represents the first exponent of the polynomial in the 2nd dimension + double solve_limits(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const double &min, const double &max, const int &axis, const int &x_exp, const int &y_exp); + + /// Uses the Newton solver to find the roots of p(x_in,y_in)-z_in + /// @param coefficients vector containing the ordered coefficients + /// @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 guess double value that represents the start value + /// @param axis unsigned integer value that represents the axis to solve for (0=x, 1=y) + /// @param x_exp integer value that represents the first exponent of the polynomial in the 1st dimension + /// @param y_exp integer value that represents the first exponent of the polynomial in the 2nd dimension + double solve_guess(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const double &guess, const int &axis, const int &x_exp, const int &y_exp); + +}; + +class Poly2DflexResidual : public Poly2DResidual { +protected: + int x_exp, y_exp; + /// Object that evaluates the equation + Polynomial2Dflex poly; + +private: + Poly2DflexResidual(); + +public: + Poly2DflexResidual(Polynomial2Dflex &poly, const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const int &axis, const int &x_exp, const int &y_exp); + virtual ~Poly2DResidual(){}; }; diff --git a/src/PolyMath.cpp b/src/PolyMath.cpp index f79c10c4..955e1860 100644 --- a/src/PolyMath.cpp +++ b/src/PolyMath.cpp @@ -20,30 +20,8 @@ namespace CoolProp{ /// Basic checks for coefficient vectors. /** Starts with only the first coefficient dimension - * and checks the matrix size against the parameters rows and columns. */ -/// @param rows unsigned integer value that represents the desired degree of the polynomial -bool Polynomial2D::checkCoefficients(const unsigned int rows){ - return this->checkCoefficients(this->coefficients,rows); -} -/// @param rows unsigned integer value that represents the desired degree of the polynomial in the 1st dimension -/// @param columns unsigned integer value that represents the desired degree of the polynomial in the 2nd dimension -bool Polynomial2D::checkCoefficients(const unsigned int rows, const unsigned int columns){ - return this->checkCoefficients(this->coefficients,rows, columns); -} -/// @param coefficients vector containing the ordered coefficients -/// @param rows unsigned integer value that represents the desired degree of the polynomial -bool Polynomial2D::checkCoefficients(const Eigen::MatrixXd &coefficients, const unsigned int rows){ - if (coefficients.cols() != 1) { - throw ValueError(format("You have a 2D coefficient matrix (%d,%d), please use the 2D checks. ",coefficients.rows(),coefficients.cols())); - } - if (coefficients.rows() == rows){ - return true; - } else { - throw ValueError(format("The number of coefficients %d does not match with %d. ",coefficients.rows(),rows)); - } - return false; -} - + * and checks the matrix size against the parameters rows and columns. + */ /// @param coefficients matrix containing the ordered coefficients /// @param rows unsigned integer value that represents the desired degree of the polynomial in the 1st dimension /// @param columns unsigned integer value that represents the desired degree of the polynomial in the 2nd dimension @@ -71,26 +49,32 @@ bool Polynomial2D::checkCoefficients(const Eigen::MatrixXd &coefficients, const * sake of completeness. */ /// @param coefficients matrix containing the ordered coefficients -/// @param axis unsigned integer value that represents the desired direction of integration -Eigen::MatrixXd Polynomial2D::integrateCoeffs(const Eigen::MatrixXd &coefficients, int axis = -1){ - std::size_t r = coefficients.rows(), c = coefficients.cols(); - Eigen::MatrixXd newCoefficients; +/// @param axis integer value that represents the desired direction of integration +/// @param times integer value that represents the desired order of integration +Eigen::MatrixXd Polynomial2D::integrateCoeffs(const Eigen::MatrixXd &coefficients, const int &axis = -1, const int × = 1){ + Eigen::MatrixXd oldCoefficients; + Eigen::MatrixXd newCoefficients(coefficients); + std::size_t r, c, i, j; 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.); + for (int k = 0; k < times; k++){ + oldCoefficients = Eigen::MatrixXd(newCoefficients); + r = oldCoefficients.rows(), c = oldCoefficients.cols(); + newCoefficients = Eigen::MatrixXd::Zero(r+1,c); + newCoefficients.block(1,0,r,c) = oldCoefficients.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.); + for (int k = 0; k < times; k++){ + oldCoefficients = Eigen::MatrixXd(newCoefficients); + r = oldCoefficients.rows(), c = oldCoefficients.cols(); + newCoefficients = Eigen::MatrixXd::Zero(r,c+1); + newCoefficients.block(0,1,r,c) = oldCoefficients.block(0,0,r,c); + for (i = 0; i < r; ++i) { + for (j = 0; j < c; ++j) newCoefficients(i,j+1) /= (j+1.); } } break; @@ -107,26 +91,35 @@ Eigen::MatrixXd Polynomial2D::integrateCoeffs(const Eigen::MatrixXd &coefficient * original coefficients with i and lowering the order by 1. */ /// @param coefficients matrix containing the ordered coefficients -/// @param axis unsigned integer value that represents the desired direction of derivation -Eigen::MatrixXd Polynomial2D::deriveCoeffs(const Eigen::MatrixXd &coefficients, int axis = -1){ - std::size_t r = coefficients.rows(), c = coefficients.cols(); +/// @param axis integer value that represents the desired direction of derivation +/// @param times integer value that represents the desired order of integration +Eigen::MatrixXd Polynomial2D::deriveCoeffs(const Eigen::MatrixXd &coefficients, const int &axis = -1, const int × = 1){ + if (times < 0) throw ValueError(format("You have to provide a positive order for derivation, %d is not valid. ",times)); + if (times == 0) return coefficients; + // Recursion is also possible, but not recommended + //Eigen::MatrixXd newCoefficients; + //if (times > 1) newCoefficients = deriveCoeffs(coefficients, axis, times-1); + //else newCoefficients = Eigen::MatrixXd(coefficients); Eigen::MatrixXd newCoefficients(coefficients); + std::size_t r, c, i, j; switch (axis) { case 0: - for (size_t i = 1; i < r; ++i) { - for (size_t j = 0; j < c; ++j) { - newCoefficients(i,j) *= i; + for (int k = 0; k < times; k++){ + r = newCoefficients.rows(), c = newCoefficients.cols(); + for (i = 1; i < r; ++i) { + for (j = 0; j < c; ++j) newCoefficients(i,j) *= i; } + removeRow(newCoefficients,0); } - 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; + for (int k = 0; k < times; k++){ + r = newCoefficients.rows(), c = newCoefficients.cols(); + for (i = 0; i < r; ++i) { + for (j = 1; j < c; ++j) newCoefficients(i,j) *= j; } + removeColumn(newCoefficients,0); } - removeColumn(newCoefficients,0); break; default: throw ValueError(format("You have to provide a dimension, 0 or 1, for derivation, %d is not valid. ",axis)); @@ -168,28 +161,29 @@ double Polynomial2D::evaluate(const Eigen::MatrixXd &coefficients, const double } +/// @param coefficients vector containing the ordered coefficients /// @param x_in double value that represents the current input in the 1st dimension /// @param y_in double value that represents the current input in the 2nd dimension -/// @param axis unsigned integer value that represents the axis to solve for (0=x, 1=y) -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; - } - return result; +/// @param axis unsigned integer value that represents the axis to derive for (0=x, 1=y) +double Polynomial2D::derivative(const Eigen::MatrixXd &coefficients, const double &x_in, const double &y_in, const int &axis = -1){ + return this->evaluate(this->deriveCoeffs(coefficients, axis, 1), x_in, y_in); } + + +/// @param coefficients vector containing the ordered coefficients +/// @param x_in double value that represents the current input in the 1st dimension +/// @param y_in double value that represents the current input in the 2nd dimension +/// @param axis unsigned integer value that represents the axis to integrate for (0=x, 1=y) +double Polynomial2D::integral(const Eigen::MatrixXd &coefficients, const double &x_in, const double &y_in, const int &axis = -1){ + return this->evaluate(this->integrateCoeffs(coefficients, axis, 1), x_in,y_in); +} + + +/// @param coefficients vector containing the ordered coefficients /// @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(const double &in, const double &z_in, int axis = -1){ +Eigen::VectorXd Polynomial2D::solve(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const int &axis = -1){ std::size_t r = coefficients.rows(), c = coefficients.cols(); Eigen::VectorXd tmpCoefficients; switch (axis) { @@ -215,18 +209,19 @@ double Polynomial2D::solve(const double &in, const double &z_in, int axis = -1){ std::vector rootsVec; polySolver.realRoots(rootsVec); if (this->do_debug()) std::cout << "Real roots: " << vec_to_string(rootsVec) << std::endl; - return rootsVec[0]; // TODO: implement root selection algorithm + return vec_to_eigen(rootsVec); + //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); +double Polynomial2D::solve_limits(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const double &min, const double &max, const int &axis){ + Poly2DResidual res = Poly2DResidual(*this, coefficients, in, z_in, axis); std::string errstring; double macheps = DBL_EPSILON; double tol = DBL_EPSILON*1e3; - int maxiter = 50; + int maxiter = 10; double result = Brent(res, min, max, macheps, tol, maxiter, errstring); if (this->do_debug()) std::cout << "Brent solver message: " << errstring << std::endl; return result; @@ -236,42 +231,18 @@ double Polynomial2D::solve_limits(const double &in, const double &z_in, const do /// @param z_in double value that represents the current output in the 3rd dimension /// @param guess double value that represents the start value /// @param axis unsigned integer value that represents the axis to solve for (0=x, 1=y) -double Polynomial2D::solve_guess(const double &in, const double &z_in, const double &guess, const int &axis){ - Poly2DResidual res = Poly2DResidual(*this, in, axis, z_in); +double Polynomial2D::solve_guess(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const double &guess, const int &axis){ + Poly2DResidual res = Poly2DResidual(*this, coefficients, in, z_in, axis); std::string errstring; //set_debug_level(1000); double tol = DBL_EPSILON*1e3; - int maxiter = 50; + int maxiter = 10; double result = Newton(res, guess, tol, maxiter, errstring); if (this->do_debug()) std::cout << "Newton solver message: " << errstring << std::endl; return result; } -// 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 * based on the length of the coefficient vector. @@ -321,47 +292,168 @@ double Polynomial2D::baseHorner(std::vector< std::vector > const& coeffi -Poly2DResidual::Poly2DResidual(const Polynomial2D &poly, const double &fixed, const int &targetDim, const double &output){ - this->poly = poly; - this->fixed = fixed; - - switch (targetDim) { +Poly2DResidual::Poly2DResidual(Polynomial2D &poly, const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const int &axis){ + switch (axis) { case iX: case iY: - this->targetDim = targetDim; + this->axis = axis; break; default: - throw ValueError(format("You have to provide a dimension to the solver, %d is not valid. ",targetDim)); + throw ValueError(format("You have to provide a dimension to the solver, %d is not valid. ",axis)); 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.dzdx(in,fixed); - if (targetDim==iY) return poly.dzdy(fixed,in); - return -_HUGE; -} - - -/// Set the coefficient matrix. -/// @param coefficients matrix containing the ordered coefficients -void Polynomial2D::setCoefficients(const Eigen::MatrixXd &coefficients){ + this->poly = poly; this->coefficients = coefficients; - this->coefficientsDerX = this->deriveCoeffs(coefficients,0); - this->coefficientsDerY = this->deriveCoeffs(coefficients,1); + this->derIsSet = false; + this->in = in; + this->z_in = z_in; } -void Polynomial2D::setCoefficients(const std::vector > &coefficients){ - this->setCoefficients(vec_to_eigen(coefficients)); + +double Poly2DResidual::call(double target){ + if (axis==iX) return poly.evaluate(coefficients,target,in)-z_in; + if (axis==iY) return poly.evaluate(coefficients,in,target)-z_in; + return -_HUGE; } +double Poly2DResidual::deriv(double target){ + if (!this->derIsSet) { + this->coefficientsDer = poly.deriveCoeffs(coefficients,axis); + this->derIsSet = true; + } + return poly.evaluate(coefficientsDer,target,in); +} + + + + + + + + + +/// Integration functions +/** Integrating coefficients for polynomials is done by dividing the + * original coefficients by (i+1) and elevating the order by 1 + * through adding a zero as first coefficient. + * Some reslicing needs to be applied to integrate along the x-axis. + * In the brine/solution equations, reordering of the parameters + * avoids this expensive operation. However, it is included for the + * sake of completeness. + */ +/// @param coefficients matrix containing the ordered coefficients +/// @param axis unsigned integer value that represents the desired direction of integration +/// @param times integer value that represents the desired order of integration +/// @param firstExponent integer value that represents the first exponent of the polynomial in axis direction +Eigen::MatrixXd Polynomial2Dflex::integrateCoeffs(const Eigen::MatrixXd &coefficients, const int &axis, const int ×, const int &firstExponent){ + Eigen::MatrixXd oldCoefficients; + Eigen::MatrixXd newCoefficients(coefficients); + std::size_t r, c, i, j; + switch (axis) { + case 0: + for (int k = 0; k < times; k++){ + oldCoefficients = Eigen::MatrixXd(newCoefficients); + r = oldCoefficients.rows(), c = oldCoefficients.cols(); + newCoefficients = Eigen::MatrixXd::Zero(r+1,c); + newCoefficients.block(1,0,r,c) = oldCoefficients.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: + for (int k = 0; k < times; k++){ + oldCoefficients = Eigen::MatrixXd(newCoefficients); + r = oldCoefficients.rows(), c = oldCoefficients.cols(); + newCoefficients = Eigen::MatrixXd::Zero(r,c+1); + newCoefficients.block(0,1,r,c) = oldCoefficients.block(0,0,r,c); + for (i = 0; i < r; ++i) { + for (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; + } + return newCoefficients; +} + + +/// Derivative coefficients calculation +/** Deriving coefficients for polynomials is done by multiplying the + * original coefficients with i and lowering the order by 1. + */ +/// @param coefficients matrix containing the ordered coefficients +/// @param axis unsigned integer value that represents the desired direction of derivation +/// @param times integer value that represents the desired order of derivation +/// @param firstExponent integer value that represents the first exponent of the polynomial in axis direction +Eigen::MatrixXd Polynomial2Dflex::deriveCoeffs(const Eigen::MatrixXd &coefficients, const int &axis, const int ×, const int &firstExponent){ + if (times < 0) throw ValueError(format("You have to provide a positive order for derivation, %d is not valid. ",times)); + if (times == 0) return coefficients; + // Recursion is also possible, but not recommended + //Eigen::MatrixXd newCoefficients; + //if (times > 1) newCoefficients = deriveCoeffs(coefficients, axis, times-1); + //else newCoefficients = Eigen::MatrixXd(coefficients); + Eigen::MatrixXd newCoefficients(coefficients); + std::size_t r, c, i, j; + switch (axis) { + case 0: + for (int k = 0; k < times; k++){ + r = newCoefficients.rows(), c = newCoefficients.cols(); + for (i = 1; i < r; ++i) { + for (j = 0; j < c; ++j) newCoefficients(i,j) *= i; + } + removeRow(newCoefficients,0); + } + break; + case 1: + for (int k = 0; k < times; k++){ + r = newCoefficients.rows(), c = newCoefficients.cols(); + for (i = 0; i < r; ++i) { + for (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; + } + return newCoefficients; +} + + + + + + + +Poly2DflexResidual::Poly2DflexResidual(Polynomial2Dflex &poly, const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const int &axis, const int &x_exp, const int &y_exp) + :Poly2DResidual(poly, coefficients, in, z_in, axis){ + this->x_exp = x_exp; + this->y_exp = y_exp; +} + +double Poly2DflexResidual::call(double target){ + if (axis==iX) return poly.evaluate(coefficients,target,in,x_exp,y_exp)-z_in; + if (axis==iY) return poly.evaluate(coefficients,in,target,x_exp,y_exp)-z_in; + return -_HUGE; +} + +double Poly2DflexResidual::deriv(double target){ + if (!this->derIsSet) { + if (axis==iX) this->coefficientsDer = poly.deriveCoeffs(coefficients,axis,1,x_exp); + if (axis==iY) this->coefficientsDer = poly.deriveCoeffs(coefficients,axis,1,y_exp); + this->derIsSet = true; + } + return poly.evaluate(coefficientsDer,target,in,x_exp,y_exp); +} + + + + + } @@ -394,7 +486,7 @@ void Polynomial2D::setCoefficients(const std::vector > &coef TEST_CASE("Internal consistency checks and example use cases for PolyMath.cpp","[PolyMath]") { - bool PRINT = false; + bool PRINT = true; std::string tmpStr; /// Test case for "SylthermXLT" by "Dow Chemicals" @@ -420,45 +512,39 @@ TEST_CASE("Internal consistency checks and example use cases for PolyMath.cpp"," Eigen::MatrixXd matrix2Dtmp; std::vector > vec2Dtmp; - CoolProp::Polynomial2D poly2D; - - SECTION("Coefficient parsing and setting") { - CHECK_NOTHROW(poly2D.setCoefficients(cHeat2D)); - CHECK_THROWS(poly2D.checkCoefficients(4,5)); - CHECK_NOTHROW(poly2D.checkCoefficients(3,4)); - CHECK_NOTHROW(poly2D.setCoefficients(matrix2D)); - CHECK_THROWS(poly2D.checkCoefficients(4,5)); - CHECK_NOTHROW(poly2D.checkCoefficients(3,4)); - + SECTION("Coefficient parsing") { + CoolProp::Polynomial2D poly; + CHECK_THROWS(poly.checkCoefficients(matrix2D,4,5)); + CHECK( poly.checkCoefficients(matrix2D,3,4) ); } SECTION("Coefficient operations") { Eigen::MatrixXd matrix; CoolProp::Polynomial2D poly; - CHECK_THROWS(poly2D.integrateCoeffs(matrix2D)); + CHECK_THROWS(poly.integrateCoeffs(matrix2D)); - CHECK_NOTHROW(matrix = poly2D.integrateCoeffs(matrix2D, 0)); + CHECK_NOTHROW(matrix = poly.integrateCoeffs(matrix2D, 0)); tmpStr = CoolProp::mat_to_string(matrix2D); if (PRINT) std::cout << tmpStr << std::endl; tmpStr = CoolProp::mat_to_string(matrix); if (PRINT) std::cout << tmpStr << std::endl << std::endl; - CHECK_NOTHROW(matrix = poly2D.integrateCoeffs(matrix2D, 1)); + CHECK_NOTHROW(matrix = poly.integrateCoeffs(matrix2D, 1)); tmpStr = CoolProp::mat_to_string(matrix2D); if (PRINT) std::cout << tmpStr << std::endl; tmpStr = CoolProp::mat_to_string(matrix); if (PRINT) std::cout << tmpStr << std::endl << std::endl; - CHECK_THROWS(poly2D.deriveCoeffs(matrix2D)); + CHECK_THROWS(poly.deriveCoeffs(matrix2D)); - CHECK_NOTHROW(matrix = poly2D.deriveCoeffs(matrix2D, 0)); + CHECK_NOTHROW(matrix = poly.deriveCoeffs(matrix2D, 0)); tmpStr = CoolProp::mat_to_string(matrix2D); if (PRINT) std::cout << tmpStr << std::endl; tmpStr = CoolProp::mat_to_string(matrix); if (PRINT) std::cout << tmpStr << std::endl << std::endl; - CHECK_NOTHROW(matrix = poly2D.deriveCoeffs(matrix2D, 1)); + CHECK_NOTHROW(matrix = poly.deriveCoeffs(matrix2D, 1)); tmpStr = CoolProp::mat_to_string(matrix2D); if (PRINT) std::cout << tmpStr << std::endl; tmpStr = CoolProp::mat_to_string(matrix); @@ -468,12 +554,12 @@ TEST_CASE("Internal consistency checks and example use cases for PolyMath.cpp"," SECTION("Evaluation and test values"){ Eigen::MatrixXd matrix = CoolProp::vec_to_eigen(cHeat); - CoolProp::Polynomial2D poly(matrix); + CoolProp::Polynomial2D poly; double acc = 0.0001; double T = 273.15+50; - double c = poly.evaluate(T, 0.0); + double c = poly.evaluate(matrix, T, 0.0); double d = 1834.746; { @@ -486,7 +572,7 @@ TEST_CASE("Internal consistency checks and example use cases for PolyMath.cpp"," } c = 2.0; - c = poly.solve_x(0.0, d); + c = poly.solve(matrix, 0.0, d, 0)[0]; { CAPTURE(T); CAPTURE(c); @@ -495,7 +581,7 @@ TEST_CASE("Internal consistency checks and example use cases for PolyMath.cpp"," } c = 2.0; - c = poly.solve_limits_x(0.0, d, -50, 750); + c = poly.solve_limits(matrix, 0.0, d, -50, 750, 0); { CAPTURE(T); CAPTURE(c); @@ -504,7 +590,7 @@ TEST_CASE("Internal consistency checks and example use cases for PolyMath.cpp"," } c = 2.0; - c = poly.solve_guess_x(0.0, d, 350); + c = poly.solve_guess(matrix, 0.0, d, 350, 0); { CAPTURE(T); CAPTURE(c); @@ -528,16 +614,15 @@ TEST_CASE("Internal consistency checks and example use cases for PolyMath.cpp"," SECTION("Integration and derivation tests") { + CoolProp::Polynomial2D poly; + Eigen::MatrixXd matrix(matrix2D); - CoolProp::Polynomial2D poly(matrix); + Eigen::MatrixXd matrixInt = poly.integrateCoeffs(matrix, 1); + Eigen::MatrixXd matrixDer = poly.deriveCoeffs(matrix, 1); + Eigen::MatrixXd matrixInt2 = poly.integrateCoeffs(matrix, 1, 2); + Eigen::MatrixXd matrixDer2 = poly.deriveCoeffs(matrix, 1, 2); - Eigen::MatrixXd matrixInt = poly.integrateCoeffs(matrix, 1); - CoolProp::Polynomial2D polyInt(matrixInt); - - Eigen::MatrixXd matrixDer = poly.deriveCoeffs(matrix, 1); - CoolProp::Polynomial2D polyDer(matrixDer); - - CHECK_THROWS( poly2D.evaluate(matrix,0.0) ); + CHECK_THROWS( poly.evaluate(matrix,0.0) ); double x = 0.3, y = 255.3, val1, val2, val3, val4; @@ -548,10 +633,10 @@ TEST_CASE("Internal consistency checks and example use cases for PolyMath.cpp"," double acc = 0.001; for (double T = Tmin; T