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"

std::vector<std::vector<double> > FuncWrapperND::Jacobian(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, std::vector<double> x0, double tol, int maxiter, std::string *errstring)
{
        int iter=0;
        *errstring=std::string("");
        std::vector<double> f0,v,negative_f0;
        std::vector<std::vector<double> > J;
        double error = 999;
        while (iter==0 || fabs(error)>tol){
                f0 = f->call(x0);
                J = f->Jacobian(x0);
                
                // Negate f0
                negative_f0 = f0;
                for (unsigned int i = 0; i<f0.size(); i++){ negative_f0[i] *= -1;}
                // find v from J*v = -f
                v = linsolve(J, negative_f0);
                // Update the guess
        for (unsigned int i = 0; i<v.size(); i++){ x0[i] += v[i];}
                error = root_sum_square(f0);
                if (iter>maxiter){
                        *errstring=std::string("reached maximum number of iterations");
                        x0[0]=_HUGE;
                }
                iter++;
        }
        return x0;
}

double Newton(FuncWrapper1D *f, double x0, double ftol, int maxiter, std::string *errstring)
{
        double x, dx, fval=999;
    int iter=1;
        (*errstring).empty();
        x = x0;
    while ((iter < 2 || fabs(fval) > ftol) && iter < maxiter)
    {
                fval = f->call(x);
                dx = -fval/f->deriv(x);
                x += dx;
                
                if (fabs(dx/x) < 10*DBL_EPSILON)
                {
                        return x;
                }

                if (iter>maxiter)
                {
                        *errstring=std::string("reached maximum number of iterations");
                        throw SolutionError(format("Newton reached maximum number of iterations"));
                }
        iter=iter+1;
    }
    return x;
}

double Secant(FuncWrapper1D *f, double x0, double dx, double tol, int maxiter, std::string *errstring)
{
        double x1=0,x2=0,x3=0,y1=0,y2=0,x,fval=999;
    int iter=1;
        *errstring=std::string("");
        
        if (fabs(dx)==0){ *errstring=std::string("dx cannot be zero"); return _HUGE;}
    while ((iter<=2 || fabs(fval)>tol) && iter<maxiter)
    {
        if (iter==1){x1=x0; x=x1;}
        if (iter==2){x2=x0+dx; x=x2;}
        if (iter>2) {x=x2;}
                        fval=f->call(x);
        if (iter==1){y1=fval;}
        if (iter>1)
        {
                        double deltax = x2-x1;
                        if (fabs(deltax)<1e-14) 
                        { 
                                if (fabs(fval) < tol*10)
                                {
                                        return x; 
                                }
                                else
                                {
                                        throw ValueError("Step is small but not solved to tolerance");
                                }
                        }
            y2=fval;
            x3=x2-y2/(y2-y1)*(x2-x1);
            y1=y2; x1=x2; x2=x3;
                        
        }
                if (iter>maxiter)
                {
                        *errstring=std::string("reached maximum number of iterations");
                        throw SolutionError(format("Secant reached maximum number of iterations"));
                }
        iter=iter+1;
    }
    return x3;
}

double BoundedSecant(FuncWrapper1D *f, double x0, double xmin, double xmax, double dx, double tol, int maxiter, std::string *errstring)
{
        double x1=0,x2=0,x3=0,y1=0,y2=0,x,fval=999;
    int iter=1;
        *errstring=std::string("");
        
        if (fabs(dx)==0){ *errstring=std::string("dx cannot be zero"); return _HUGE;}
    while ((iter<=3 || fabs(fval)>tol) && iter<100)
    {
        if (iter==1){x1=x0; x=x1;}
        if (iter==2){x2=x0+dx; x=x2;}
        if (iter>2) {x=x2;}
                        fval=f->call(x);
        if (iter==1){y1=fval;}
        if (iter>1)
        {
            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)
                {
                        *errstring=std::string("reached maximum number of iterations");
                        throw SolutionError(format("BoundedSecant reached maximum number of iterations"));
                }
        iter=iter+1;
    }
    return x3;
}

double Brent(FuncWrapper1D *f, double a, double b, double macheps, double t, int maxiter, std::string *errstr)
{
        int iter;
        *errstr=std::string("");
        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 (fabs(fb) < t) { return b;}
        if (!ValidNumber(fb)){
                throw ValueError(format("Brent's method f(b) is NAN for b = %g",b).c_str());
        }
        if (fabs(fa) < t) { return a;}
        if (!ValidNumber(fa)){
                throw ValueError(format("Brent's method f(a) is NAN for a = %g",a).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 (fabs(fc)<fabs(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*fabs(b)+t;
        while (fabs(m)>tol && fb!=0){
        // See if a bisection is forced
                if (fabs(e)<tol || fabs(fa) <= fabs(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-fabs(tol*q) || p<fabs(0.5*s*q)){
                d=p/q;
                        }
                        else{
                m=0.5*(c-b);
                d=e=m;
                        }
                }
        a=b;
        fa=fb;
                if (fabs(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 (fb*fc>0){
            // Goto int: from Brent ALGOL code
            c=a;
            fc=fa;
            d=e=b-a;
                }
                if (fabs(fc)<fabs(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*fabs(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(std::string("Brent's method reached maximum number of steps of %d ", maxiter));}
        }
    return b;
}