gsl_chebapp.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_CHEBAPP_H
#define O2SCL_GSL_CHEBAPP_H

#include <o2scl/collection.h>
#include <o2scl/funct.h>
#include <gsl/gsl_chebyshev.h>
#include <cmath>

#ifndef DOXYGENP
namespace o2scl {
#endif

  /** 
      \brief Chebyshev approximation (GSL)
      
      Approximate a function using a Chebyshev series:
      \f[
      f(x) = \sum_n c_n T_n(x) \qquad \mathrm{where}
      \qquad T_n(x)=\cos(n \arccos x)
      \f]

      \future Implement eval_err(), eval_n() and eval_n_err() methods
      from GSL.
  */
  template<class param_t, class func_t> class gsl_chebapp {
  public:
  
    gsl_chebapp() {
      order=5;
      init_called=false;
      gcs=gsl_cheb_alloc(order);
    }
    
    /** \brief Initialize a Chebyshev approximation of the function
        \c func over the interval from \c a to \c b

        The interval must be specified so that \f$ a < b \f$ .
    */
    int init(func_t &func, double a, double b, param_t &vp) {
      size_t k, j;
  
      if(a>=b) {
        O2SCL_ERR_RET("Interval a>=b in gsl_chebapp::init().",gsl_einval);
      }
      gcs->a=a;
      gcs->b=b;
  
      double bma=0.5*(gcs->b-gcs->a);
      double bpa=0.5*(gcs->b+gcs->a);
      double fac=2.0/(gcs->order+1.0);
  
      for(k=0;k<=gcs->order;k++) {
        double y=cos(M_PI*(k+0.5)/(gcs->order+1));
        gcs->f[k]=func(y*bma+bpa,vp);
      }
  
      for(j=0;j<=gcs->order;j++) {
        double sum=0.0;
        for(k=0;k<=gcs->order; k++) {
          sum+=gcs->f[k]*cos(M_PI*j*(k+0.5)/(gcs->order+1));
        }
        gcs->c[j]=fac*sum;
      }

      init_called=true;

      return 0;
    }

    /** 
        \brief Set the order (default 5)
        
        The function init() must be called after calling set_order()
        to reinitialize the series for the new order.
     */
    int set_order(size_t od) {
      order=od;
      init_called=false;
      gsl_cheb_free(gcs);
      gcs=gsl_cheb_alloc(order);
      return 0;
    }

    /** \brief Evaluate the approximation
     */
    double eval(double x) {
      if (init_called==false) {
        O2SCL_ERR("Series not initialized in gsl_chebapp::eval().",
                gsl_einval);
        return 0.0;
      }
      return gsl_cheb_eval(gcs,x);
    }
    
    /** 
        \brief Return a pointer to an approximation to the derivative
        
        The new gsl_chebapp object is allocated by new, and the memory
        should be deallocated using delete by the user.
     */
    gsl_chebapp *deriv() {
      if (init_called==false) {
        O2SCL_ERR("Series not initialized in gsl_chebapp::deriv().",
                gsl_einval);
        return 0;
      }
      gsl_chebapp *gc=new gsl_chebapp;
      gc->set_order(order);
      gc->init_called=true;
      gsl_cheb_calc_deriv(gc->gcs,gcs);
      return gc;
    }

    /** 
        \brief Return a pointer to an approximation to the integral

        The new gsl_chebapp object is allocated by new, and the memory
        should be deallocated using delete by the user.
     */
    gsl_chebapp *inte() {
      if (init_called==false) {
        O2SCL_ERR("Series not initialized in gsl_chebapp::inte().",
                gsl_einval);
        return 0;
      }
      gsl_chebapp *gc=new gsl_chebapp;
      gc->set_order(order);
      gc->init_called=true;
      gsl_cheb_calc_integ(gc->gcs,gcs);
      return gc;
    }

    /** 
        \brief Get the coefficient

        Legal values of the argument are 0 to \c order (inclusive)
    */
    double get_coefficient(size_t ix) {
      if (ix<order+1) {
        return gcs->c[ix];
      }
      O2SCL_ERR
        ("Requested invalid coefficient in gsl_chebapp::get_coefficient()",
         gsl_einval);
      return 0.0;
    }

#ifndef DOXYGENP

  protected:

    /// The GSL representation of the series
    gsl_cheb_series *gcs;

    /// The order
    size_t order;

    /// True if init has been called
    bool init_called;

#endif
  
  };

#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