Solvers.cpp

Go to the documentation of this file.

#include <vector>
#include "Solvers.h"
#include "math.h"
#include "MatrixMath.h"
#include <iostream>
#include "CoolPropTools.h"
#include <Eigen/Dense>
 
namespace CoolProp {
 
std::vector<std::vector<double>> FuncWrapperND::Jacobian(const std::vector<double>& x) {
    double epsilon;
    std::size_t N = x.size();
    std::vector<double> r, xp;
    std::vector<std::vector<double>> J(N, std::vector<double>(N, 0));
    std::vector<double> r0 = call(x);
    // Build the Jacobian by column
    for (std::size_t i = 0; i < N; ++i) {
        xp = x;
        epsilon = 0.001 * x[i];
        xp[i] += epsilon;
        r = call(xp);
 
        for (std::size_t j = 0; j < N; ++j) {
            J[j][i] = (r[j] - r0[j]) / epsilon;
        }
    }
    return J;
}
 
std::vector<double> NDNewtonRaphson_Jacobian(FuncWrapperND* f, const std::vector<double>& x, double tol, int maxiter, double w) {
    int iter = 0;
    f->errstring.clear();
    std::vector<double> f0, v;
    std::vector<std::vector<double>> JJ;
    std::vector<double> x0 = x;
    Eigen::VectorXd r(x0.size());
    Eigen::MatrixXd J(x0.size(), x0.size());
    double error = 999;
    while (iter == 0 || std::abs(error) > tol) {
        f0 = f->call(x0);
        JJ = f->Jacobian(x0);
 
        for (std::size_t i = 0; i < x0.size(); ++i) {
            r(i) = f0[i];
            for (std::size_t j = 0; j < x0.size(); ++j) {
                J(i, j) = JJ[i][j];
            }
        }
 
        Eigen::VectorXd v = J.colPivHouseholderQr().solve(-r);
 
        // Update the guess
        double max_relchange = -1;
        for (std::size_t i = 0; i < x0.size(); i++) {
            x0[i] += w * v(i);
            double relchange = std::abs(v(i) / x0[i]);
            if (std::abs(x0[i]) > 1e-16 && relchange > max_relchange) {
                max_relchange = relchange;
            }
        }
 
        // Stop if the solution is not changing by more than numerical precision
        double max_abschange = v.cwiseAbs().maxCoeff();
        if (max_abschange < DBL_EPSILON * 100) {
            return x0;
        }
        if (max_relchange < 1e-12) {
            return x0;
        }
        error = root_sum_square(f0);
        if (iter > maxiter) {
            f->errstring = "reached maximum number of iterations";
            x0[0] = _HUGE;
        }
        iter++;
    }
    return x0;
}
 
double Newton(FuncWrapper1DWithDeriv* f, double x0, double ftol, int maxiter) {
    double x, dx, dfdx, fval = 999;
 
    // Initialize
    f->iter = 0;
    f->errstring.clear();
    x = x0;
    while (f->iter < maxiter || std::abs(fval) > ftol) {
        fval = f->call(x);
        dfdx = f->deriv(x);
        dx = -fval / dfdx;
 
        if (!ValidNumber(fval)) {
            f->errstring = "Residual function in newton returned invalid number";
            throw ValueError("Residual function in newton returned invalid number");
        };
 
        if (f->verbosity > 0) {
            std::cout << format("i: %d, x: %0.15g, dx: %g, f: %g, dfdx: %g", f->iter, x, dx, fval, dfdx) << std::endl;
        }
 
        x += dx;
 
        if (std::abs(dx / x) < 1e-11) {
            return x;
        }
 
        if (f->iter > maxiter) {
            f->errstring = "Newton reached maximum number of iterations";
            throw SolutionError(format("Newton reached maximum number of iterations"));
        }
        f->iter++;
    }
    return x;
}
double Halley(FuncWrapper1DWithTwoDerivs* f, double x0, double ftol, int maxiter, double xtol_rel) {
    double x, dx, fval = 999, dfdx, d2fdx2;
 
    // Initialize
    f->iter = 0;
    f->errstring.clear();
    x = x0;
 
    // The relaxation factor (less than 1 for smaller steps)
    double omega = f->options.get_double("omega", 1.0);
 
    while (f->iter < 2 || std::abs(fval) > ftol) {
        if (f->input_not_in_range(x)) {
            std::string msg = format("Input [%0.15g] is out of range", x);
            f->errstring = msg;
            throw ValueError(msg);
        }
 
        fval = f->call(x);
        dfdx = f->deriv(x);
        d2fdx2 = f->second_deriv(x);
 
        if (!ValidNumber(fval)) {
            f->errstring = "Residual function in Halley returned invalid number";
            throw ValueError("Residual function in Halley returned invalid number");
        };
        if (!ValidNumber(dfdx)) {
            f->errstring = "Derivative function in Halley returned invalid number";
            throw ValueError("Derivative function in Halley returned invalid number");
        };
 
        dx = -omega * (2 * fval * dfdx) / (2 * POW2(dfdx) - fval * d2fdx2);
 
        if (f->verbosity > 0) {
            std::cout << format("i: %d, x: %0.15g, dx: %g, f: %g, dfdx: %g, d2fdx2: %g", f->iter, x, dx, fval, dfdx, d2fdx2) << std::endl;
        }
 
        x += dx;
 
        if (std::abs(dx / x) < xtol_rel) {
            return x;
        }
 
        if (f->iter > maxiter) {
            f->errstring = "Halley reached maximum number of iterations";
            throw SolutionError(format("Halley reached maximum number of iterations"));
        }
        f->iter += 1;
    }
    return x;
}
 
double Householder4(FuncWrapper1DWithThreeDerivs* f, double x0, double ftol, int maxiter, double xtol_rel) {
    double x, dx, fval = 999, dfdx, d2fdx2, d3fdx3;
 
    // Initialization
    f->iter = 1;
    f->errstring.clear();
    x = x0;
 
    // The relaxation factor (less than 1 for smaller steps)
    double omega = f->options.get_double("omega", 1.0);
 
    while (f->iter < 2 || std::abs(fval) > ftol) {
        if (f->input_not_in_range(x)) {
            throw ValueError(format("Input [%g] is out of range", x));
        }
 
        fval = f->call(x);
        dfdx = f->deriv(x);
        d2fdx2 = f->second_deriv(x);
        d3fdx3 = f->third_deriv(x);
 
        if (!ValidNumber(fval)) {
            throw ValueError("Residual function in Householder4 returned invalid number");
        };
        if (!ValidNumber(dfdx)) {
            throw ValueError("Derivative function in Householder4 returned invalid number");
        };
        if (!ValidNumber(d2fdx2)) {
            throw ValueError("Second derivative function in Householder4 returned invalid number");
        };
        if (!ValidNumber(d3fdx3)) {
            throw ValueError("Third derivative function in Householder4 returned invalid number");
        };
 
        dx = -omega * fval * (POW2(dfdx) - fval * d2fdx2 / 2.0) / (POW3(dfdx) - fval * dfdx * d2fdx2 + d3fdx3 * POW2(fval) / 6.0);
 
        x += dx;
 
        if (std::abs(dx / x) < xtol_rel) {
            return x;
        }
 
        if (f->iter > maxiter) {
            f->errstring = "reached maximum number of iterations";
            throw SolutionError(format("Householder4 reached maximum number of iterations"));
        }
        f->iter += 1;
    }
    return x;
}
 
double Secant(FuncWrapper1D* f, double x0, double dx, double tol, int maxiter) {
#if defined(COOLPROP_DEEP_DEBUG)
    static std::vector<double> xlog, flog;
    xlog.clear();
    flog.clear();
#endif
 
    // Initialization
    double x1 = 0, x2 = 0, x3 = 0, y1 = 0, y2 = 0, x = x0, fval = 999;
    f->iter = 1;
    f->errstring.clear();
 
    // The relaxation factor (less than 1 for smaller steps)
    double omega = f->options.get_double("omega", 1.0);
 
    if (std::abs(dx) == 0) {
        f->errstring = "dx cannot be zero";
        return _HUGE;
    }
    while (f->iter <= 2 || std::abs(fval) > tol) {
        if (f->iter == 1) {
            x1 = x0;
            x = x1;
        }
        if (f->iter == 2) {
            x2 = x0 + dx;
            x = x2;
        }
        if (f->iter > 2) {
            x = x2;
        }
 
        if (f->input_not_in_range(x)) {
            throw ValueError(format("Input [%g] is out of range", x));
        }
 
        fval = f->call(x);
 
#if defined(COOLPROP_DEEP_DEBUG)
        xlog.push_back(x);
        flog.push_back(fval);
#endif
 
        if (!ValidNumber(fval)) {
            throw ValueError("Residual function in secant returned invalid number");
        };
        if (f->iter == 1) {
            y1 = fval;
        }
        if (f->iter > 1) {
            double deltax = x2 - x1;
            if (std::abs(deltax) < 1e-14) {
                return x;
            }
            y2 = fval;
            double deltay = y2 - y1;
            if (f->iter > 2 && std::abs(deltay) < 1e-14) {
                return x;
            }
            x3 = x2 - omega * y2 / (y2 - y1) * (x2 - x1);
            y1 = y2;
            x1 = x2;
            x2 = x3;
        }
        if (f->iter > maxiter) {
            f->errstring = std::string("reached maximum number of iterations");
            throw SolutionError(format("Secant reached maximum number of iterations"));
        }
        f->iter += 1;
    }
    return x3;
}
 
double BoundedSecant(FuncWrapper1D* f, double x0, double xmin, double xmax, double dx, double tol, int maxiter) {
    double x1 = 0, x2 = 0, x3 = 0, y1 = 0, y2 = 0, x, fval = 999;
    int iter = 1;
    f->errstring.clear();
    if (std::abs(dx) == 0) {
        f->errstring = "dx cannot be zero";
        return _HUGE;
    }
    while (iter <= 3 || std::abs(fval) > tol) {
        if (iter == 1) {
            x1 = x0;
            x = x1;
        } else if (iter == 2) {
            x2 = x0 + dx;
            x = x2;
        } else {
            x = x2;
        }
        fval = f->call(x);
        if (iter == 1) {
            y1 = fval;
        } else {
            y2 = fval;
            x3 = x2 - y2 / (y2 - y1) * (x2 - x1);
            // Check bounds, go half the way to the limit if limit is exceeded
            if (x3 < xmin) {
                x3 = (xmin + x2) / 2;
            }
            if (x3 > xmax) {
                x3 = (xmax + x2) / 2;
            }
            y1 = y2;
            x1 = x2;
            x2 = x3;
        }
        if (iter > maxiter) {
            f->errstring = "reached maximum number of iterations";
            throw SolutionError(format("BoundedSecant reached maximum number of iterations"));
        }
        iter = iter + 1;
    }
    f->errcode = 0;
    return x3;
}
 
double ExtrapolatingSecant(FuncWrapper1D* f, double x0, double dx, double tol, int maxiter) {
#if defined(COOLPROP_DEEP_DEBUG)
    static std::vector<double> xlog, flog;
    xlog.clear();
    flog.clear();
#endif
 
    // Initialization
    double x1 = 0, x2 = 0, x3 = 0, y0 = 0, y1 = 0, y2 = 0, x = x0, fval = 999;
    f->iter = 1;
    f->errstring.clear();
 
    // The relaxation factor (less than 1 for smaller steps)
    double omega = f->options.get_double("omega", 1.0);
 
    if (std::abs(dx) == 0) {
        f->errstring = "dx cannot be zero";
        return _HUGE;
    }
    while (f->iter <= 2 || std::abs(fval) > tol) {
        if (f->iter == 1) {
            x1 = x0;
            x = x1;
        }
        if (f->iter == 2) {
            x2 = x0 + dx;
            x = x2;
        }
        if (f->iter > 2) {
            x = x2;
        }
 
        if (f->input_not_in_range(x)) {
            throw ValueError(format("Input [%g] is out of range", x));
        }
 
        fval = f->call(x);
 
#if defined(COOLPROP_DEEP_DEBUG)
        xlog.push_back(x);
        flog.push_back(fval);
#endif
 
        if (!ValidNumber(fval)) {
            if (f->iter == 1) {
                return x;
            } else {
                return x2 - omega * y1 / (y1 - y0) * (x2 - x1);
            }
        };
        if (f->iter == 1) {
            y1 = fval;
        }
        if (f->iter > 1) {
            double deltax = x2 - x1;
            if (std::abs(deltax) < 1e-14) {
                return x;
            }
            y2 = fval;
            double deltay = y2 - y1;
            if (f->iter > 2 && std::abs(deltay) < 1e-14) {
                return x;
            }
            x3 = x2 - omega * y2 / (y2 - y1) * (x2 - x1);
            y0 = y1;
            y1 = y2;
            x1 = x2;
            x2 = x3;
        }
        if (f->iter > maxiter) {
            f->errstring = std::string("reached maximum number of iterations");
            throw SolutionError(format("Secant reached maximum number of iterations"));
        }
        f->iter += 1;
    }
    return x3;
}
 
double Brent(FuncWrapper1D* f, double a, double b, double macheps, double t, int maxiter) {
    int iter;
    f->errstring.clear();
    double fa, fb, c, fc, m, tol, d, e, p, q, s, r;
    fa = f->call(a);
    fb = f->call(b);
 
    // If one of the boundaries is to within tolerance, just stop
    if (std::abs(fb) < t) {
        return b;
    }
    if (!ValidNumber(fb)) {
        throw ValueError(format("Brent's method f(b) is NAN for b = %g, other input was a = %g", b, a).c_str());
    }
    if (std::abs(fa) < t) {
        return a;
    }
    if (!ValidNumber(fa)) {
        throw ValueError(format("Brent's method f(a) is NAN for a = %g, other input was b = %g", a, b).c_str());
    }
    if (fa * fb > 0) {
        throw ValueError(format("Inputs in Brent [%f,%f] do not bracket the root.  Function values are [%f,%f]", a, b, fa, fb));
    }
 
    c = a;
    fc = fa;
    iter = 1;
    if (std::abs(fc) < std::abs(fb)) {
        // Goto ext: from Brent ALGOL code
        a = b;
        b = c;
        c = a;
        fa = fb;
        fb = fc;
        fc = fa;
    }
    d = b - a;
    e = b - a;
    m = 0.5 * (c - b);
    tol = 2 * macheps * std::abs(b) + t;
    while (std::abs(m) > tol && fb != 0) {
        // See if a bisection is forced
        if (std::abs(e) < tol || std::abs(fa) <= std::abs(fb)) {
            m = 0.5 * (c - b);
            d = e = m;
        } else {
            s = fb / fa;
            if (a == c) {
                //Linear interpolation
                p = 2 * m * s;
                q = 1 - s;
            } else {
                //Inverse quadratic interpolation
                q = fa / fc;
                r = fb / fc;
                m = 0.5 * (c - b);
                p = s * (2 * m * q * (q - r) - (b - a) * (r - 1));
                q = (q - 1) * (r - 1) * (s - 1);
            }
            if (p > 0) {
                q = -q;
            } else {
                p = -p;
            }
            s = e;
            e = d;
            m = 0.5 * (c - b);
            if (2 * p < 3 * m * q - std::abs(tol * q) || p < std::abs(0.5 * s * q)) {
                d = p / q;
            } else {
                m = 0.5 * (c - b);
                d = e = m;
            }
        }
        a = b;
        fa = fb;
        if (std::abs(d) > tol) {
            b += d;
        } else if (m > 0) {
            b += tol;
        } else {
            b += -tol;
        }
        fb = f->call(b);
        if (!ValidNumber(fb)) {
            throw ValueError(format("Brent's method f(t) is NAN for t = %g", b).c_str());
        }
        if (std::abs(fb) < macheps) {
            return b;
        }
        if (fb * fc > 0) {
            // Goto int: from Brent ALGOL code
            c = a;
            fc = fa;
            d = e = b - a;
        }
        if (std::abs(fc) < std::abs(fb)) {
            // Goto ext: from Brent ALGOL code
            a = b;
            b = c;
            c = a;
            fa = fb;
            fb = fc;
            fc = fa;
        }
        m = 0.5 * (c - b);
        tol = 2 * macheps * std::abs(b) + t;
        iter += 1;
        if (!ValidNumber(a)) {
            throw ValueError(format("Brent's method a is NAN").c_str());
        }
        if (!ValidNumber(b)) {
            throw ValueError(format("Brent's method b is NAN").c_str());
        }
        if (!ValidNumber(c)) {
            throw ValueError(format("Brent's method c is NAN").c_str());
        }
        if (iter > maxiter) {
            throw SolutionError(format("Brent's method reached maximum number of steps of %d ", maxiter));
        }
        if (std::abs(fb) < 2 * macheps * std::abs(b)) {
            return b;
        }
    }
    return b;
}
 
}; /* namespace CoolProp */