gsl_root_stef.h

/*
  -------------------------------------------------------------------
  
  Copyright (C) 2006, 2007, 2008, Andrew W. Steiner
  
  This file is part of O2scl.
  
  O2scl 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 3 of the License, or
  (at your option) any later version.
  
  O2scl 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 O2scl. If not, see <http://www.gnu.org/licenses/>.

  -------------------------------------------------------------------
*/
#ifndef O2SCL_GSL_ROOT_STEF_H
#define O2SCL_GSL_ROOT_STEF_H

#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>
#include <o2scl/root.h>

#ifndef DOXYGENP
namespace o2scl {
#endif

  /** 
      \brief Steffenson equation solver (GSL)

      This class finds a root of a function a derivative. If the
      derivative is not analytically specified, it is most likely
      preferable to use of the alternatives, gsl_root_brent,
      cern_root, or cern_mroot_root. The function solve_de() performs
      the solution automatically, and a lower-level GSL-like interface
      with set() and iterate() is also provided.

      By default, this solver compares the present value of the root
      (\f$ \mathrm{root} \f$) to the previous value (\f$ \mathrm{x} \f$), 
      and returns
      success if \f$ | \mathrm{root} - \mathrm{x} | < \mathrm{tol} \f$, 
      where \f$ \mathrm{tol}=\mathrm{tolx}+\mathrm{tolf2}~\mathrm{root} \f$ .

      If \ref test_residual is set to true, then the solver 
      additionally requires that the absolute value of the function
      is less than \ref root::tolf.

      The original variable \c x_2 has been removed as it was unused
      in the original GSL code.

      \future There's some extra copying here which can probably
      be removed.
      \future Compare directly to GSL.
      
  */
  template<class param_t, class func_t, class dfunc_t> class gsl_root_stef : 
  public o2scl::root_de<param_t,func_t,dfunc_t> {
      
  protected:
      
    /// Function value
    double f;

    /// Derivative value
    double df;

    /// Previous value of root
    double x_1;

    /// Root
    double x;

    /// Number of iterations
    int count;

    /// The function to solve
    func_t *fp;

    /// The derivative
    dfunc_t *dfp;
      
    /// The function parameters
    param_t *params;
      
  public:
  
    gsl_root_stef() {
      test_residual=false;
      tolf2=1.0e-12;
    }
    
    /// Return the type, \c "gsl_root_stef".
    virtual const char *type() { return "gsl_root_stef"; }

    /// The present solution estimate
    double root;

    /** \brief The relative tolerance for subsequent solutions 
        (default \f$ 10^{-12} \f$)
    */
    double tolf2;

    /**
       \brief Perform an iteration

       After a successful iteration, \ref root contains the
       most recent value of the root. 
    */
    int iterate() {
      
      double x_new, f_new, df_new;
        
      double x_1t = x_1;
      double xt = x;
        
      if (df == 0.0) {
        this->last_conv=gsl_ezerodiv;
        O2SCL_CONV_RET("Derivative is zero in gsl_root_stef::iterate().",
                       o2scl::gsl_ezerodiv,this->err_nonconv);
      }
        
      x_new = xt - (f / df);

      // It is important that the derivative be evaluated first here,
      // because we want the last function evaluation to be an
      // evaluation for the returned root
      (*dfp)(x_new,df_new,*params);
      (*fp)(x_new,f_new,*params);
    
      x_1 = xt;
      x = x_new;
        
      f = f_new;
      df = df_new;
    
      if (!finite (f_new)) {
        O2SCL_ERR_RET("Function not continuous in gsl_root_stef::iterate().",
                      o2scl::gsl_ebadfunc);
      }
        
      if (count < 3) {

        root = x_new;
        count++;

      } else {

        double u = (xt - x_1t);
        double v = (x_new - 2 * xt + x_1t);
        
        if (v == 0) {
          root = x_new;  /* avoid division by zero */
        } else {
          root = x_1t - u * u / v;  /* accelerated value */
        }
      }
      
      if (!finite(df_new)) {
        O2SCL_ERR_RET
          ("Function not differentiable in gsl_root_stef::iterate().",
           o2scl::gsl_ebadfunc);
      }
      
      return o2scl::gsl_success;
    }

    /** \brief Solve \c func using \c x as an initial
        guess using derivatives \c df.
    */
    virtual int solve_de(double &xx, param_t &pa, func_t &fun, dfunc_t &dfun) {
      
      int status1, status2=gsl_continue, iter=0;
        
      status1=set(fun,dfun,xx,pa);

      this->last_conv=0;
      while (status1==gsl_success && status2==gsl_continue && 
             iter<this->ntrial) {
        iter++;

        status1=iterate();

        // Compare present value to previous value
        status2=gsl_root_test_delta(root,xx,this->tolx,tolf2);

        if (test_residual && status2==gsl_success) {
          double y;
          fun(root,y,pa);
          if (fabs(y)>=this->tolf) status2=gsl_continue;
        }

        if (this->verbose>0) {
          double fval;
          fun(root,fval,pa);
          print_iter(root,fval,iter,fabs(root-xx),this->tolx*root,
                     "gsl_root_stef");
        }
        xx=root;
      }
      
      this->last_ntrial=iter;

      if (status1!=gsl_success || status2!=gsl_success) {
        int ret=o2scl::err_hnd->get_errno();
        return ret;
      }
      if (iter>=this->ntrial) {
        if (this->last_conv==0) this->last_conv=gsl_emaxiter;
        O2SCL_CONV_RET("solve_de() exceeded maximum number of iterations.",
                       gsl_emaxiter,this->err_nonconv);
      }

      return o2scl::gsl_success;
    }
    
    /// True if we should test the residual also (default false)
    bool test_residual;

    /** 
        \brief Set the information for the solver

        Set the function, the derivative, the initial guess and
        the parameters.
    */
    int set(func_t &fun, dfunc_t &dfun, double guess, param_t &pa) {
    
      fp=&fun;
      dfp=&dfun;
      params=&pa;
      root=guess;

      dfun(root,df,pa);
      fun(root,f,pa);

      x=root;
      x_1=0.0;
      count=1;
        
      return o2scl::gsl_success;
        
    }

  };
  
#ifndef DOXYGENP
  template<> int io_tlate<gsl_root_stef<int,funct<int>,
    funct<int> > >::input(cinput *co, in_file_format *ins, 
                          gsl_root_stef<int, funct<int>,
                          funct<int> > *ro);
  template<> int io_tlate<gsl_root_stef<int,funct<int>,
    funct<int> > >::output(coutput *co, out_file_format *outs, 
                           gsl_root_stef<int, funct<int>,
                           funct<int> > *ro);
  template<> const char *io_tlate<gsl_root_stef<int,
    funct<int>, funct<int> > >::type();
#endif
  
  typedef io_tlate<gsl_root_stef<int,funct<int>, funct<int> > > 
    gsl_root_stef_io_type;
 
#ifndef DOXYGENP
}
#endif

#endif

---

Documentation generated with Doxygen and provided under the GNU Free Documentation License. See License Information for details.

Project hosting provided by , O2scl Sourceforge Project Page