zeroin.h

Go to the documentation of this file.

/*

freesteam - IAPWS-IF97 steam tables library
Copyright (C) 2004-2005  John Pye

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

*/

#ifndef ZEROINVISITOR_H
#define ZEROINVISITOR_H

#include "designbycontract.h"

#include <iostream>
#include <cmath>

#ifndef DBL_EPSILON
        #define DBL_EPSILON 2e-16
#endif


template < class Subject, class Ordinate = double, class Abscissa = double >
class ZeroIn : public DesignByContract {

    private:
    Abscissa ax;                
    Abscissa bx;                
    Ordinate(Subject::*method) (const Abscissa&);       
    Abscissa tol;               
    Abscissa x;                 
    bool issolved;

    Abscissa a, b, c;           
    Ordinate fa;                
    Ordinate fb;                
    Ordinate fc;                

    Subject *subject;

    public:

    ZeroIn() {
            this->tol = 0.0 * Abscissa();
    }

    void setLowerBound(const Abscissa& ax) {
            this->ax = ax;
    }

    Abscissa getLowerBound() const {
            return ax;
    }

    void setUpperBound(const Abscissa& bx) {
            this->bx = bx;
    }

    Abscissa getUpperBound() const{
            return bx;
    }

    void setTolerance(const Abscissa& tol) {
            this->tol = tol;
    }


    void setMethod(Ordinate(Subject::*method) (const Abscissa&)) {
            this->method = method;
    }


    void visit(Subject *subject) {

            //cerr << "Begin zeroin with lower bound " << K_TO_C(getLowerBound()) << "°C and upper bound " << K_TO_C(getUpperBound()) << endl;

            this->subject = subject;
            this->issolved = false;

            a = ax;
            b = bx;
            fa = ((*subject).*method) (a);
            fb = ((*subject).*method) (b);
            c = a;
            fc = fa;

                if(fa == 0.0 * Ordinate()){
                        fb=0.0 * Ordinate(); // used by getError
                        this->issolved = true; // perfect solution
                        return;
                }

            //  Main iteration loop

            for (;;) {
                    Abscissa prev_step = b - a;    
                    Abscissa tol_act;                
                    Abscissa p;                    
                    double q;                      
                    Abscissa new_step;             

                    if (fabs(fc) < fabs(fb)) {
                            a = b;
                            b = c;
                            c = a;  //  Swap data for b to be the best approximation
                            fa = fb;
                            fb = fc;
                            fc = fa;
                    }

                    // DBL_EPSILON is defined in math.h
                    tol_act = 2.0* DBL_EPSILON * fabs(b) + tol / 2.0;

                    new_step = (c - b) / 2.0;

                    if (fabs(new_step) <= tol_act || fb == 0.0 * Ordinate()) {
                            this->issolved = true;       //  Acceptable approx. is found
                            //cerr << "Best solution is b, f(b=" << b << ")=" << fb <<",f(a="<<a<<")="<< fa <<endl;
                            return;
                    }
                    //  Decide if the interpolation can be tried

                    if (fabs(prev_step) >= tol_act   // If prev_step was large enough and was in true direction,
                        && fabs(fa) > fabs(fb))      // Interpolatiom may be tried
                    {
                            register double t1, t2;
                            Abscissa cb;

                            cb = c - b;
                            if (a == c) {
                                    // If we have only two distinct points
                                    // then only linear interpolation can be applied
                                    t1 = fb / fa;
                                    p = cb * t1;
                                    q = 1.0 - t1;
                            } else {
                                    // Quadric inverse interpolation

                                    q = fa / fc;
                                    t1 = fb / fc;
                                    t2 = fb / fa;
                                    p = t2 * (cb * q * (q - t1) - (b - a) * (t1 - 1.0));
                                    q = (q - 1.0) * (t1 - 1.0) * (t2 - 1.0);
                            }
                            if (p > 0.0 * Abscissa()) {
                                    // p was calculated with the opposite sign; make p positive-
                                    q = -q;     // and assign possible minus to q
                            } else {
                                    p = -p;
                            }

                            if (p < (0.75 * cb * q - fabs(tol_act * q) / 2.0)
                                && p < fabs(prev_step * q / 2.0)
                               ) {
                                    // If b+p/q falls in [b,c] and
                                    // isn't too large it is accepted
                                    new_step = p / q;
                            }
                            // If p/q is too large then the bissection procedure can
                            // reduce [b,c] range to more extent
                    }

                    if (fabs(new_step) < tol_act) {     // Adjust the step to be not less
                            if (new_step > 0.0 * Abscissa())    // than tolerance
                                    new_step = tol_act;
                            else
                                    new_step = -tol_act;
                    }

                    a = b;
                    fa = fb;            // Save the previous approx.
                    b += new_step;
                    fb = ((*subject).*method) (b);      //  Do step to a new approxim.

                    if ((fb > 0.0 * Ordinate() && fc > 0.0 * Ordinate())
                        || (fb < 0.0 * Ordinate() && fc < 0.0 * Ordinate())) {
                            c = a;
                            fc = fa;    //  Adjust c for it to have a sign opposite to that of b
                    }
            }
            // (((we never arrive here)))
    }

    bool isSolved(void) const{
            return issolved;
    }

    bool isSolved(const Ordinate &error) {
            if (fabs(fb) > error) {
                    return false;
            }
            return true;
    }

    Abscissa getSolution(void) const {
            return a;
    }

    Ordinate getError(void) const{
            return fb;
    }
};

#endif

---