github/CoolProp
1
1
Fork 0
mirror of https://github.com/CoolProp/CoolProp.git synced 2026-04-01 03:00:13 -04:00
Code Issues Packages Projects Releases Wiki Activity

Started to work on flexible polynomial solvers for fractional exponents

Browse Source
This commit is contained in:
jowr
2014-06-13 13:56:08 +02:00
parent 4a02014439
commit 4e2aa1c404
4 changed files with 456 additions and 371 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
include/CoolPropTools.h
View File

@@ -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);

16
include/MatrixMath.h
View File

@@ -35,16 +35,16 @@ namespace CoolProp{
* it is still needed.
*/
/// @param coefficients matrix containing the ordered coefficients
template <typename T> std::vector<T> eigen_to_vec1D(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> &coefficients, int axis = 1){
template <typename T> std::vector<T> eigen_to_vec1D(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> &coefficients, int axis = 0){
std::vector<T> 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 <class T> Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> vec_to_eigen(c
return result;
}
/// @param coefficients matrix containing the ordered coefficients
template <class T> Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> vec_to_eigen(const std::vector<T> &coefficients, int axis = 1){
template <class T> Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> vec_to_eigen(const std::vector<T> &coefficients, int axis = 0){
size_t nRows = num_rows(coefficients);
Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> 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;
}

360
include/PolyMath.h
View File

@@ -11,135 +11,32 @@
//#include <numeric> // inner_product
//#include <sstream>
//#include "float.h"
#include "MatrixMath.h"
#include <unsupported/Eigen/Polynomials>
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<double> &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<double> &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<double> &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<double> integrateCoeffs(const std::vector<double> &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<double> deriveCoeffs(const std::vector<double> &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<std::vector<double> > &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<double> &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<std::vector<double> > &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<std::vector<double> > &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 &times);
/// 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 &times);
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 &times, 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 &times, 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(){};
};

447
src/PolyMath.cpp
View File

@@ -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 &times = 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 &times = 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<double> 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<double> > 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<std::vector<double> > &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 &times, 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 &times, 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<std::vector<double> > &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<std::vector<double> > 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<Tmax; T+=Tinc) {
val1 = poly.evaluate(x, T-deltaT);
val2 = poly.evaluate(x, T+deltaT);
val1 = poly.evaluate(matrix, x, T-deltaT);
val2 = poly.evaluate(matrix, x, T+deltaT);
val3 = (val2-val1)/2/deltaT;
val4 = polyDer.evaluate(x, T);
val4 = poly.evaluate(matrixDer, x, T);
CAPTURE(T);
CAPTURE(val3);
CAPTURE(val4);
@@ -563,10 +648,25 @@ TEST_CASE("Internal consistency checks and example use cases for PolyMath.cpp","
}
for (double T = Tmin; T<Tmax; T+=Tinc) {
val1 = polyInt.evaluate(x, T-deltaT);
val2 = polyInt.evaluate(x, T+deltaT);
val1 = poly.evaluate(matrixDer, x, T-deltaT);
val2 = poly.evaluate(matrixDer, x, T+deltaT);
val3 = (val2-val1)/2/deltaT;
val4 = poly.evaluate(x, T);
val4 = poly.evaluate(matrixDer2, x, T);
CAPTURE(T);
CAPTURE(val3);
CAPTURE(val4);
tmpStr = CoolProp::mat_to_string(matrixDer);
CAPTURE(tmpStr);
tmpStr = CoolProp::mat_to_string(matrixDer2);
CAPTURE(tmpStr);
CHECK( check_abs(val3,val4,acc) );
}
for (double T = Tmin; T<Tmax; T+=Tinc) {
val1 = poly.evaluate(matrixInt, x, T-deltaT);
val2 = poly.evaluate(matrixInt, x, T+deltaT);
val3 = (val2-val1)/2/deltaT;
val4 = poly.evaluate(matrix, x, T);
CAPTURE(T);
CAPTURE(val3);
CAPTURE(val4);
@@ -578,8 +678,23 @@ TEST_CASE("Internal consistency checks and example use cases for PolyMath.cpp","
}
for (double T = Tmin; T<Tmax; T+=Tinc) {
val1 = poly.evaluate(x, T);
val2 = polyInt.derivative(x, T, 1);
val1 = poly.evaluate(matrixInt2, x, T-deltaT);
val2 = poly.evaluate(matrixInt2, x, T+deltaT);
val3 = (val2-val1)/2/deltaT;
val4 = poly.evaluate(matrixInt, x, T);
CAPTURE(T);
CAPTURE(val3);
CAPTURE(val4);
tmpStr = CoolProp::mat_to_string(matrixInt2);
CAPTURE(tmpStr);
tmpStr = CoolProp::mat_to_string(matrixInt);
CAPTURE(tmpStr);
CHECK( check_abs(val3,val4,acc) );
}
for (double T = Tmin; T<Tmax; T+=Tinc) {
val1 = poly.evaluate(matrix, x, T);
val2 = poly.derivative(matrixInt, x, T, 1);
CAPTURE(T);
CAPTURE(val1);
CAPTURE(val2);
@@ -587,8 +702,8 @@ TEST_CASE("Internal consistency checks and example use cases for PolyMath.cpp","
}
for (double T = Tmin; T<Tmax; T+=Tinc) {
val1 = poly.derivative(x, T, 1);
val2 = polyDer.evaluate(x, T);
val1 = poly.derivative(matrix, x, T, 1);
val2 = poly.evaluate(matrixDer, x, T);
CAPTURE(T);
CAPTURE(val1);
CAPTURE(val2);