NewtonRaphson.h

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2009 by Clark Sims                                      *
 *   http://AcumenSoftwareInc.com/WhoWeAre/Clark_Sims.html                 *
 *   ClarkSims@AcumenSoftwareInc.com                                       *
 *                                                                         *
 *   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 NewtonRaphson_hpp
#define NewtonRaphson_hpp

#include <iostream>
#include <vector>

using namespace std;

template<class functor, class real>
class NewtonRaphsonSolve0 {
 protected:
  // key attributes

  int _max_iter;

  real _x0;

  real _epsilon;

  real _fsolve;

  bool _twice_df;

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

  real (functor::*_df)(real);

  real (functor::*_d2f)(real);

  functor& _F;

  // 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_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_d2f( real (functor::*d2f)(real)) {
    _d2f = d2f;
  }

  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 _F;}

   
  NewtonRaphsonSolve0( 
    real X0, 
    real Epsilon, 
    bool Twice_df,
    unsigned int max_iter, 
    functor& F, 
    real (functor::*f)(real), 
    real (functor::*df)(real),  
    real (functor::*d2f)(real) = NULL,
    real fsolve = 0,
    bool check_boundary = false,
    real min_x = 0,
    real max_x = 0)
    :  _x0( X0), _epsilon( Epsilon), _twice_df( Twice_df), _max_iter(max_iter), _F(F), _f(f), _df(df), _d2f(d2f), _fsolve(fsolve), _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_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( ostream* output) {
    int i;

    _converged = _diverged = _boundary_solution = false;

    _x[0] = _x0;

    for (i=_number_boundary_hits=0; ;) {
      _F_x[i]       = (_F.*_f)( _x[i]) - _fsolve;

      _abs_F_x[i] = (_F_x[i] > 0)? _F_x[i] : -_F_x[i];

      if (_abs_F_x[i] < _epsilon) {
        _x_converge = _x[i];
        _f_converge = _F_x[i];
        _converged = true;
        break;
      }

      if (i>0 && _abs_F_x[i] > _abs_F_x[i-1]) {
        _diverged = true;
        break;
      }

      if (output) {
        *output << i << " " << _x[i] << " " << _F_x[i] << " " << _dF_x[i] << endl;
      } 

      _dF_x[i] = (_F.*_df)( _x[i]);

      ++i;
      if (i >= _max_iter) break;

      if (_d2f != 0) {  
        real dx = (_F_x[i-1]) / _dF_x[i-1];
        real df2_dx2 = (_F.*_d2f) ( _x[i-1]);
        _dF_x[i-1] -= .5 * dx * df2_dx2;
      } else if (_twice_df) {
        _x[i] = _x[i-1] - (_F_x[i-1]) / _dF_x[i-1];
        if (check_boundary( _x[i], _dF_x[i-1])>2) {
          goto boundary;
        }

        if (_boundary_hit) {
          _dF_x[i-1] = (_F.*_df)( _x[i]);
        } else {
          real df_dx = (_F.*_df)( _x[i]);
          _dF_x[i-1] = .5 * (_dF_x[i-1] + df_dx);
        }
      }

      _x[i] = _x[i-1] - (_F_x[i-1]) / _dF_x[i-1];
      
      if (check_boundary( _x[i], _dF_x[i-1])>2) {
        goto boundary;
      }

      continue;

    boundary:
      if (_lower_boundary_hit) {
        _boundary_solution = true;
        _x_converge = _min_x;
        _f_converge = (_F.*_f)( _min_x);
        break;
      }

      if (_upper_boundary_hit) {
        _boundary_solution = true;
        _x_converge = _max_x;
        _f_converge = (_F.*_f)( _max_x);
        break;
      }
    }

    if (_diverged && output) {
      *output << "answer diverged" << endl;
    }

    if (_boundary_solution && output) {
      *output << i << " " << _x_converge << " " << _f_converge << " boundary_solution" << endl;
    } else if (_converged && output) {
      *output << i << " " << _x[i] << " " << _F_x[i] << " answer converged" << endl;
    }
  }
};


#endif

---