Secant.h

Go to the documentation of this file.

#ifndef Secant_hpp
#define Secant_hpp

#include <iostream>
#include <vector>

using namespace std;

template<class functor, class real>
class SecantSolve0 {
 protected:
  // key attributes
  int _max_iter;

  real _x0;
  real _x1;

  real _Fsolve;
  
  real _epsilon;

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

  functor& _obj;

  // functional attributes
  bool _check_boundary;
  real _min_x;
  real _max_x;


  // derived attributes
  mutable real _x_converge;
  mutable real _f_converge;
  mutable vector<real> _x;
  mutable vector<real> _F_x;
  mutable vector<real> _abs_Err_F_x;
  mutable vector<real> _dF_x;
  mutable bool _converged;
  mutable bool _diverged;
  mutable bool _boundary_solution;
  mutable int _number_boundary_hits;
  mutable bool _lower_boundary_hit;
  mutable bool _upper_boundary_hit;
  mutable bool _boundary_hit;

  int check_boundary( real & x, real df_dx) {

    if (_check_boundary) {
      _lower_boundary_hit = _upper_boundary_hit = false;
      if (x < _min_x) {
        x = _min_x + 2 * _epsilon / df_dx;
        _lower_boundary_hit = true;
        _number_boundary_hits++;      
      } else if (x > _max_x) {
        x = _max_x - 2 * _epsilon / df_dx;
        _upper_boundary_hit = true;
        _number_boundary_hits++;
      }
      _boundary_hit = _lower_boundary_hit || _upper_boundary_hit;
    }
    
    return _number_boundary_hits;
  }

 public:  

  bool get_boundary_solution() const {
    return _boundary_solution;
  }


  void set_check_boundary( bool check_boundary) {
    _check_boundary = check_boundary;
  }

  void set_min_x( real min_x) {
    _min_x = min_x;
  }

  void set_max_x( real max_x) {
    _max_x = max_x;
  }

  real get_x_converge() const {
    return _x_converge;
  }

  real get_f_converge() const {
    return _f_converge;
  }

  functor& get_functor() { return _obj;}

  bool get_converged() const { return _converged;}

  SecantSolve0( 
    real X0, 
    real X1,
    real Fsolve,
    real Epsilon, 
    unsigned int max_iter, 
    functor& obj, 
    real (functor::*f)(real), 
    bool check_boundary = false,
    real min_x = 0,
    real max_x = 0)
    :  _x0( X0), _x1( X1), _Fsolve( Fsolve), _epsilon( Epsilon),  _max_iter(max_iter), _obj(obj), _f(f),  _check_boundary(check_boundary), _min_x(min_x), _max_x(max_x), 
    _x_converge(0), _f_converge(0), 
    _x(max_iter), _F_x(max_iter), _abs_Err_F_x(max_iter), _dF_x(max_iter), _converged(false), _diverged(false), 
    _lower_boundary_hit(false), _upper_boundary_hit(false), _boundary_hit(false)
  {

  }

  void do_iteration( bool print_results) {
    int i;
    real dx, abs_dx, err;

    _converged = _diverged = _boundary_solution = false;

    _x[0] = _x0;
    _x[1] = _x1;
    _F_x[0] = (_obj.*_f)( _x[0]);
    _F_x[1] = (_obj.*_f)( _x[1]);
    _abs_Err_F_x[1] = (_F_x[1] > 0)? _F_x[1] : -_F_x[1];
    _dF_x[0] = 0;
    _dF_x[1] = (_F_x[1] - _F_x[0]) / (_x[1] - _x[0]);

    if (print_results) {
      cout << 0 << " " << _x[0] << " " << _F_x[0] << " not calculated" << endl;
      cout << 1 << " " << _x[1] << " " << _F_x[1] << " " << _dF_x[1] << endl;
    } 

    _number_boundary_hits=0;
    for (i=2; i<_max_iter; i++) {

      _x[i] = _x[i-1] - (_F_x[i-1] - _Fsolve) / _dF_x[i-1];

      if (check_boundary( _x[i], _dF_x[i-1])>1) {
        if (_lower_boundary_hit) {
          _boundary_solution = true;
          _x_converge = _min_x;
          _f_converge = (_obj.*_f)( _min_x);
          break;
        }
        if (_upper_boundary_hit) {
          _boundary_solution = true;
          _x_converge = _max_x;
          _f_converge = (_obj.*_f)( _max_x);
          break;
        }
      } else {
        if (_lower_boundary_hit) {
          _x[i] = _min_x + _epsilon;
        }
        if (_upper_boundary_hit) {
          _x[i] = _max_x - _epsilon;
        }
      }

      _F_x[i] = (_obj.*_f)( _x[i]);
      err = _F_x[i] - _Fsolve;
      _abs_Err_F_x[i] = (err > 0)? err : -err;
      dx = (_x[i] - _x[i-1]);
      abs_dx = (dx > 0)? dx : -dx;
        
      if (_abs_Err_F_x[i] < _epsilon || abs_dx < _epsilon) {
        _x_converge = _x[i];
        _f_converge = _F_x[i];
        _converged = true;
        break;
      }
      if (_abs_Err_F_x[i] > _abs_Err_F_x[i-1]) {
        _diverged = true;
      }

      _dF_x[i] = (_F_x[i] - _F_x[i-1]) / (_x[i] - _x[i-1]);

      if (print_results) {
        cout << i << " " << _x[i] << " " << _F_x[i] << " " << _dF_x[i] << endl;
      } 
    }

    if (_boundary_solution && print_results) {
      cout << i << " " << _x_converge << " " << _f_converge << " boundary_solution" << endl;
    } else if (_converged && print_results) {
      cout << i << " " << _x[i] << " " << _F_x[i] << " answer converged" << endl;
    } else if (_diverged && print_results) {
      cout << "answer diverged" << endl;
    }
  }
};


#endif

---