Bisection.h

Go to the documentation of this file.

#ifndef Bisection_hpp
#define Bisection_hpp

#include <iostream>
#include <vector>

#include <stdexcept>

#include "Secant.h"

using namespace std;

template< class functor, class real>
void stradle_value( 
                   real& X0,
                   real& X1,
                   real& FX0,
                   real& FX1,
                   real Fsolve,
                   functor& obj, 
                   real (functor::*f)(real),
                   bool initialize_fx0 = true,
                   bool initialize_fx1 = true,
                   int direction = 0,
                   real MinX = 0,
                   real MaxX = 0)
{
  bool is_increasing, check_bounds;

  if (initialize_fx0) FX0 =(obj.*f)( X0);
  if (initialize_fx1) FX1 =(obj.*f)( X1);

  check_bounds = MinX!=0 || MaxX!=0;

  if (X0==X1) throw domain_error( "X0==X1 in stradle function");

  if (X0 > X1) {
    real tmp;
    tmp = X0;    X0 = X1;  X1 = tmp;
    tmp = FX0; FX0 = FX1; FX1 = tmp;
  }

  if (FX0 < Fsolve && Fsolve < FX1) return;

  if (direction == 0) {
    is_increasing = FX1 > FX0;
  } else {
    is_increasing = direction > 0;
  }

  if (check_bounds && (X1 > MaxX || X0 < MinX)) {
    throw domain_error( "X1 > MaxX || X0 < MinX in stradle_value");
  }

  if (is_increasing) {    // solve fx0 < fsolve < fx1
    while (FX0 > Fsolve) {
      if (FX1 < FX0) throw range_error( "FX1 < FX0 for monotonic increasing function  at label A in stradle_value");
      FX1 = FX0;
      X0 *= .5;
      if (check_bounds && X0 < MinX) {
        throw domain_error( "X0 < MinX in stradle_value");
      }
      FX0 = (obj.*f)( X0);
    }
    while (Fsolve > FX1) {
      if (FX1 < FX0) throw range_error( "FX1 < FX0 for monotonic increasing function  at label B in stradle_value");
      FX0 = FX1;
      X1 *= 2;
      if (check_bounds && X1 > MaxX) {
        throw domain_error( "X1 > MaxX in stradle_value");
      }
      FX1 = (obj.*f)( X1);
    }
  } else {  // solve fx0 > fsolve > fx1
    while (FX0 < Fsolve) {
      if (FX1 > FX0) throw range_error( "FX1 > FX0 for monotonic decreasing function at label C in stradle_value");
      FX1 = FX0;
      X0 *= .5;
      if (check_bounds && X0 < MinX) {
        throw domain_error( "X0 < MinX in stradle_value");
      }
      FX0 = (obj.*f)( X0);
    }
    while (Fsolve < FX1) {
      if (FX1 > FX0) throw range_error( "FX1 > FX0 for monotonic decreasing function at label D in stradle_value");
      FX0 = FX1;
      X1 *= 2;
      if (check_bounds && X1 > MaxX) {
        throw domain_error( "X1 > MaxX in stradle_value");
      }
      FX1 = (obj.*f)( X1);
    }
  }
}

template<class functor, class real>
class Bisection {
 protected:

  real _x0;

  real _x1;

  real _fsolve;

  real _delta;

  real _epsilon;

  real _x_lower;

  real _x_upper;

  real _F_x_lower;

  bool _calc_F_x_lower;

  real _F_x_upper;

  bool _calc_F_x_upper;

  real _x_mid;

  real _dx;

  real _F_x_mid;

  real _abs_error;

  bool _is_increasing;

  int _max_iter;

  int i;

  bool _converged;

  functor& _obj;

  real (functor::*_f)(real);

 public:

  Bisection( 
    real X0, 
    real X1,
    real Fsolve,
    real Delta,
    real Epsilon, 
    unsigned int max_iter, 
    functor& obj, 
    real (functor::*f)(real)) :
    _x0(X0), _x1(X1), _fsolve(Fsolve), _delta(Delta), _epsilon(Epsilon), _max_iter(max_iter), _obj(obj), _f(f), _calc_F_x_lower(false), _calc_F_x_upper(false)
  {

  }

  Bisection( 
    real X0, 
    real X1,
    real FX0,
    real FX1,
    real Fsolve,
    real Delta,
    real Epsilon, 
    unsigned int max_iter, 
    functor& obj, 
    real (functor::*f)(real)) :
    _x0(X0), _x1(X1), _F_x_lower(FX0), _calc_F_x_lower(true), 
    _F_x_upper(FX1), _calc_F_x_upper(true), _fsolve(Fsolve), 
    _delta(Delta), _epsilon(Epsilon), _max_iter(max_iter), _obj(obj), _f(f)
  {

  }

  
  void do_iteration(  bool print_results)
    {

      _converged = false;

      _x_lower = _x0;
      _x_upper = _x1;
      if (_calc_F_x_lower) _F_x_lower =  (_obj.*_f)( _x_lower); 
      if (_calc_F_x_upper) _F_x_upper =  (_obj.*_f)( _x_upper); 

      _is_increasing = (_F_x_upper > _F_x_lower);

      if (print_results) {
        cout << "x_lower = " << _x_lower << endl;
        cout << "x_upper = " << _x_upper << endl;

        cout << "F_x_lower = " << _F_x_lower << endl;
        cout << "F_x_upper = " << _F_x_upper << endl;

        if (_is_increasing) {
          cout << "F is increasing" << endl;
        } else {
          cout << "F is decreasing" << endl;
        }
      }

      _converged = false;
    
      for (i=0; i<_max_iter; i++) {
        _x_mid = .5 * (_x_lower + _x_upper);
        _dx = _x_upper - _x_mid;
        _F_x_mid = (_obj.*_f)( _x_mid);
        _abs_error = _F_x_mid - _fsolve;
        if (_abs_error < 0) _abs_error = -_abs_error;

        if (print_results) {
          cerr << "i = " << i 
               << " x_lower = " << _x_lower 
               << " x_upper = " << _x_upper 
               << " x_mid = " << _x_mid 
               << " _F_x_mid = " << _F_x_mid << endl;

        /*cout << "x_lower = " << _x_lower << endl;
          cout << "x_upper = " << _x_upper << endl;
          cout << "F_x_lower = " << _F_x_lower << endl;
          cout << "F_x_mid   = " << _F_x_mid << endl;
          cout << "F_x_upper = " << _F_x_upper << endl; */

        }

        if (_abs_error < _epsilon || _dx < _delta) {
          _converged = true;
          break;
        }

        if (_F_x_mid > _fsolve) { 
          if (_is_increasing) {
            _x_upper = _x_mid;
            _F_x_upper = _F_x_mid;
          } else {
            _x_lower = _x_mid;
            _F_x_lower = _F_x_mid;
          }
        }

        if (_F_x_mid < _fsolve) {
          if (_is_increasing) {
            _x_lower = _x_mid;
            _F_x_lower = _F_x_mid;
          } else {
            _x_upper = _x_mid;
            _F_x_upper = _F_x_mid;
          }
        }
      }
    }
                                 

  real get_x_mid() const {
    return _x_mid;
  }

  bool get_converged() const {
    return _converged;
  }

  real get_df_dx() const {
    return (_F_x_upper - _F_x_mid) / (_x_upper - _x_mid);
  }
};

template<class functor, class real>
class Bisection_Secant : public Bisection< functor, real >  {
 public:

  Bisection_Secant( 
    real X0, 
    real X1,
    real Fsolve,
    real Delta,
    real Epsilon, 
    unsigned int max_iter, 
    functor& obj, 
    real (functor::*f)(real)) :
    Bisection< functor, real > ( X0, X1, Fsolve, Delta, Epsilon, max_iter, obj, f )
  {
  }

  Bisection_Secant( 
    real X0, 
    real X1,
    real FX0,
    real FX1,
    real Fsolve,
    real Delta,
    real Epsilon, 
    unsigned int max_iter, 
    functor& obj, 
    real (functor::*f)(real)) :
    Bisection< functor, real > ( X0, X1, FX0, FX1, Fsolve, Delta, Epsilon, max_iter, obj, f )
  {
  }

  void do_iteration( bool print_results) {
    unsigned int i;
    real df_dx;

    this->_converged = false;

    this->_x_lower = this->_x0;
    this->_x_upper = this->_x1;
    if (this->_calc_F_x_lower) this->_F_x_lower =  (this->_obj.*this->_f)( this->_x_lower); 
    if (this->_calc_F_x_upper) this->_F_x_upper =  (this->_obj.*this->_f)( this->_x_upper); 
    
    if (print_results) {
      cout << "x_lower = " << this->_x_lower << endl;
      cout << "x_upper = " << this->_x_upper << endl;

      cout << "F_x_lower = " << this->_F_x_lower << endl;
      cout << "F_x_upper = " << this->_F_x_upper << endl;
    }

    this->_is_increasing = this->_F_x_upper > this->_F_x_lower;
    
    for (i=0; i<this->_max_iter; i++) {
      if (i % 2 == 1) {
        df_dx = (this->_F_x_upper - this->_F_x_lower) / (this->_x_upper - this->_x_lower);
        this->_x_mid = this->_x_upper - (this->_F_x_upper - this->_fsolve) / df_dx;
        if (this->_x_mid < this->_x_lower || this->_x_mid > this->_x_upper) {
          this->_x_mid = .5 * (this->_x_lower + this->_x_upper);
        }
      } else {
        this->_x_mid = .5 * (this->_x_lower + this->_x_upper);
      }

      this->_dx = this->_x_upper - this->_x_mid;
      this->_F_x_mid = (this->_obj.*this->_f)( this->_x_mid);
      this->_abs_error = this->_F_x_mid - this->_fsolve;
      if (this->_abs_error < 0) this->_abs_error = -this->_abs_error;

      if (print_results) {
        cerr << "i = " << i 
             << " x_lower = " << this->_x_lower 
             << " x_upper = " << this->_x_upper 
             << " x_mid = " << this->_x_mid 
             << " this->_F_x_mid = " << this->_F_x_mid << endl;
      }

      if (this->_abs_error < this->_epsilon || this->_dx < this->_delta) {
        this->_converged = true;
        break;
      }

      if (this->_F_x_mid > this->_fsolve) { 
        if (this->_is_increasing) {
          this->_x_upper = this->_x_mid;
          this->_F_x_upper = this->_F_x_mid;
        } else {
          this->_x_lower = this->_x_mid;
          this->_F_x_lower = this->_F_x_mid;
        }
      }

      if (this->_F_x_mid < this->_fsolve) {
        if (this->_is_increasing) {
          this->_x_lower = this->_x_mid;
          this->_F_x_lower = this->_F_x_mid;
        } else {
          this->_x_upper = this->_x_mid;
          this->_F_x_upper = this->_F_x_mid;
        }
      }
    }
  }
};

#endif

---