gsl_inte_qawf.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_INTE_QAWF_H
#define O2SCL_GSL_INTE_QAWF_H

#include <o2scl/inte.h>
#include <o2scl/gsl_inte_qawo.h>
#include <o2scl/gsl_inte_qagiu.h>

#ifndef DOXYGENP
namespace o2scl {
#endif
  
  /** 
      \brief Adaptive integration for oscillatory integrals (GSL)
      
      The number of subdivisions of the original interval which 
      this class is allowed to make is dictated by the workspace
      size for the integration class, which can be set using
      \ref gsl_inte_table::set_wkspace() . 
      
      \todo Improve documentation a little
  */
  template<class param_t, class func_t> class gsl_inte_qawf_sin : 
  public gsl_inte_qawo_sin<param_t,func_t> {
    
  public:
    
    gsl_inte_qawf_sin() {
    }
    
    virtual ~gsl_inte_qawf_sin() {}
    
    /** \brief Integrate function \c func from \c a to \c b.
     */
    virtual double integ(func_t &func, double a, double b, param_t &pa) {
      double res, err;
      integ_err(func,a,b,pa,res,err);
      this->interror=err;
      return res;
    }

    /** \brief Integrate function \c func from \c a to \c b and place
        the result in \c res and the error in \c err
    */
    virtual int integ_err(func_t &func, double a, double b, 
                          param_t &pa, double &res, double &err2) {
      
      this->otable=gsl_integration_qawo_table_alloc
        (this->omega,1.0,GSL_INTEG_SINE,this->tab_size);
      cyclew=gsl_integration_workspace_alloc(this->wkspace);
      
      int status=qawf(func,a,this->tolx,this->wkspace,&res,&err2,pa);
      
      gsl_integration_qawo_table_free(this->otable);
      gsl_integration_workspace_free(cyclew);
      
      return status;
      
    }

#ifndef DOXYGEN_INTERNAL

  protected:
    
    /// The integration workspace
    gsl_integration_workspace *cyclew;
    
    /** 
        \brief The full GSL integration routine called by integ_err()
    */
    int qawf(func_t &func, const double a, 
             const double epsabs, const size_t limit, 
             double *result, double *abserr, param_t &pa) {

      double area, errsum;
      double res_ext, err_ext;
      double correc, total_error = 0.0, truncation_error;

      size_t ktmin = 0;
      size_t iteration = 0;

      struct extrapolation_table table;

      double cycle;
      //double omega = this->otable->omega;

      const double p = 0.9;
      double factor = 1;
      double initial_eps, eps;
      int error_type = 0;

      /* Initialize results */

      initialise (this->w, a, a);

      *result = 0;
      *abserr = 0;

      if (limit > this->w->limit) {
        std::string estr="Iteration limit exceeds workspace ";
        estr+="in gsl_inte_qawf::qawf().";
        O2SCL_ERR_RET(estr.c_str(),gsl_einval);
      }

      /* Test on accuracy */

      if (epsabs <= 0) {
        std::string estr="The absolute tolerance must be positive ";
        estr+="in gsl_inte_qawf::qawf().";
        O2SCL_ERR_RET(estr.c_str(),gsl_ebadtol);
      }

      if (this->omega == 0.0) {
        if (this->otable->sine == GSL_INTEG_SINE) {
          /* The function sin(w x) f(x) is always zero for w = 0 */

          *result = 0;
          *abserr = 0;

          return GSL_SUCCESS;
        } else {
          /* The function cos(w x) f(x) is always f(x) for w = 0 */

          gsl_inte_qagiu<param_t,func_t> iu;
              
          int status=iu.integ_err(func,a,0.0,pa,*result,*abserr);

          return status;
        }
      }

      if (epsabs > GSL_DBL_MIN / (1 - p)) {
        eps = epsabs * (1 - p);
      } else {
        eps = epsabs;
      }

      initial_eps = eps;

      area = 0;
      errsum = 0;

      res_ext = 0;
      err_ext = GSL_DBL_MAX;
      correc = 0;
      
      cycle = (2 * floor (fabs (this->omega)) + 1) * M_PI / fabs (this->omega);

      gsl_integration_qawo_table_set_length (this->otable, cycle);

      initialise_table (&table);

      for (iteration = 0; iteration < limit; iteration++) {
        double area1, error1, reseps, erreps;

        double a1 = a + iteration * cycle;
        double b1 = a1 + cycle;

        double epsabs1 = eps * factor;

        int status=qawo(func,a1,epsabs1,0.0,limit,cyclew,this->otable,
                        &area1,&error1,pa);
          
        this->append_interval (this->w, a1, b1, area1, error1);
          
        factor *= p;

        area = area + area1;
        errsum = errsum + error1;

        /* estimate the truncation error as 50 times the final term */

        truncation_error = 50 * fabs (area1);

        total_error = errsum + truncation_error;

        if (total_error < epsabs && iteration > 4) {
          goto compute_result;
        }

        if (error1 > correc) {
          correc = error1;
        }

        if (status) {
          eps = GSL_MAX_DBL (initial_eps, correc * (1.0 - p));
        }

        if (status && total_error < 10 * correc && iteration > 3) {
          goto compute_result;
        }

        append_table (&table, area);

        if (table.n < 2) {
          continue;
        }

        qelg (&table, &reseps, &erreps);

        ktmin++;

        if (ktmin >= 15 && err_ext < 0.001 * total_error) {
          error_type = 4;
        }

        if (erreps < err_ext) {
          ktmin = 0;
          err_ext = erreps;
          res_ext = reseps;

          if (err_ext + 10 * correc <= epsabs)
            break;
          if (err_ext <= epsabs && 10 * correc >= epsabs)
            break;
        }

      }

      if (iteration == limit) error_type = 1;

      if (err_ext == GSL_DBL_MAX)
        goto compute_result;

      err_ext = err_ext + 10 * correc;

      *result = res_ext;
      *abserr = err_ext;

      if (error_type == 0) {
        return GSL_SUCCESS ;
      }

      if (res_ext != 0.0 && area != 0.0) {
        if (err_ext / fabs (res_ext) > errsum / fabs (area))
          goto compute_result;
      } else if (err_ext > errsum) {
        goto compute_result;
      } else if (area == 0.0) {
        goto return_error;
      }

      if (error_type == 4) {
        err_ext = err_ext + truncation_error;
      }

      goto return_error;

    compute_result:

      *result = area;
      *abserr = total_error;

    return_error:

      if (error_type > 2)
        error_type--;

      if (error_type == 0) {
        return GSL_SUCCESS;
      } else if (error_type == 1) {
        std::string estr="Number of iterations was insufficient ";
        estr+=" in gsl_inte_qawf::qawf().";
        O2SCL_ERR_RET(estr.c_str(),gsl_emaxiter);
      } else if (error_type == 2) {
        std::string estr="Roundoff error prevents tolerance ";
        estr+="from being achieved in gsl_inte_qawf::qawf().";
        O2SCL_ERR_RET(estr.c_str(),gsl_eround);
      } else if (error_type == 3) {
        std::string estr="Bad integrand behavior ";
        estr+=" in gsl_inte_qawf::qawf().";
        O2SCL_ERR_RET(estr.c_str(),gsl_esing);
      } else if (error_type == 4) {
        std::string estr="Roundoff error detected in extrapolation table ";
        estr+="in gsl_inte_qawf::qawf().";
        O2SCL_ERR_RET(estr.c_str(),gsl_eround);
      } else if (error_type == 5) {
        std::string estr="Integral is divergent or slowly convergent ";
        estr+="in gsl_inte_qawf::qawf().";
        O2SCL_ERR_RET(estr.c_str(),gsl_ediverge);
      } else {
        std::string estr="Could not integrate function in gsl_inte_qawf";
        estr+="::qawf() (it may have returned a non-finite result).";
        O2SCL_ERR_RET(estr.c_str(),gsl_efailed);
      }
    }
    

    /// Add the oscillating part to the integrand
    virtual double transform(func_t &func, double x, param_t &pa) {
      double wx = this->omega * x;
      double sinwx = sin(wx) ;
      double y;
      func(x,y,pa);
      return y*sinwx;
    }

#endif
  
    /// Return string denoting type ("gsl_inte_qawf_sin")
    const char *type() { return "gsl_inte_qawf_sin"; }
  
  };

  /** \brief Adaptive integration a function with finite limits of 
      integration (GSL)

      The number of subdivisions of the original interval which 
      this class is allowed to make is dictated by the workspace
      size for the integration class, which can be set using
      \ref gsl_inte_table::set_wkspace() . 
      
      \todo Verbose output has been setup for this class, but this
      needs to be done for the other GSL-like integrators
  */
  template<class param_t, class func_t> class gsl_inte_qawf_cos : 
  public gsl_inte_qawf_sin<param_t,func_t> {
    
  public:

    gsl_inte_qawf_cos() {
    }

    virtual ~gsl_inte_qawf_cos() {}
      
    /** \brief Integrate function \c func from \c a to \c b and place
        the result in \c res and the error in \c err
    */
    virtual int integ_err(func_t &func, double a, double b, 
                          param_t &pa, double &res, double &err2) {

      this->otable=gsl_integration_qawo_table_alloc
        (this->omega,b-a,GSL_INTEG_COSINE,this->tab_size);
      this->cyclew=gsl_integration_workspace_alloc(this->wkspace);
      
      int status=qawf(func,a,this->tolx,this->wkspace,&res,&err2,pa);
      
      gsl_integration_qawo_table_free(this->otable);
      gsl_integration_workspace_free(this->cyclew);
      
      return status;
      
    }

#ifndef DOXYGEN_INTERNAL

  protected:
    
    /// Add the oscillating part to the integrand
    virtual double transform(func_t &func, double x, param_t &pa) {
      double wx = this->omega * x;
      double coswx = cos(wx) ;
      double y;
      func(x,y,pa);
      return y*coswx;
    }

#endif
  
    /// Return string denoting type ("gsl_inte_qawf_cos")
    const char *type() { return "gsl_inte_qawf_cos"; }
  
  };

#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