gsl_astep.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_ASTEP_H
#define O2SCL_GSL_ASTEP_H

#include <gsl/gsl_math.h>
#include <gsl/gsl_odeiv.h>

#include <o2scl/adapt_step.h>

#ifndef DOXYGENP
namespace o2scl {
#endif

  /** 
      \brief Control structure for gsl_astep

      This class implements both the "standard" and "scaled" 
      step control methods from GSL. The standard control method
      is the default. To use the scaled controle, set \ref standard
      to <tt>false</tt> and set the scale for each component using
      set_scale().
      
      The control object is a four parameter heuristic based on
      absolute and relative errors \ref eps_abs and \ref eps_rel, and
      scaling factors \ref a_y and \ref a_dydt for the system state
      \f$ y(t) \f$ and derivatives \f$ y^{\prime}(t) \f$ respectively.

      The step-size adjustment procedure for this method begins by
      computing the desired error level \f$ D_i \f$ for each
      component. In the unscaled version,
      \f[
      D_i = \mathrm{eps\_abs}+\mathrm{eps\_rel} \times
      \left( \mathrm{a\_y} | y_i| + \mathrm{a\_dydt}~h 
      | y_i^{\prime}| \right)
      \f]
      while in the scaled version the user specifies the scale for
      each component, \f$ s_i \f$,
      \f[
      D_i = \mathrm{eps\_abs}~s_i+\mathrm{eps\_rel} \times
      \left( \mathrm{a\_y} | y_i| + \mathrm{a\_dydt}~h 
      | y_i^{\prime}| \right)
      \f]

      The desired error level \f$ D_i \f$ is compared to then observed
      error \f$ E_i = |\mathrm{yerr}_i| \f$. If the observed error \f$
      E \f$ exceeds the desired error level \f$ D \f$ by more than 10
      percent for any component then the method reduces the step-size
      by an appropriate factor,
      \f[
      h_{\mathrm{new}} = S~h_{\mathrm{old}} 
      \left(\frac{E}{D}\right)^{-1/q}
      \f]
      where \f$ q \f$ is the consistency order of the method (e.g. \f$
      q=4 \f$ for 4(5) embedded RK), and \f$ S \f$ is a safety factor
      of 0.9. The ratio \f$ E/D \f$ is taken to be the maximum of the
      ratios \f$ E_i/D_i \f$.
      
      If the observed error E is less than 50 percent of the desired
      error level \f$ D \f$ for the maximum ratio \f$ E_i/D_i \f$ then
      the algorithm takes the opportunity to increase the step-size to
      bring the error in line with the desired level,
      \f[
      h_{\mathrm{new}} = S~h_{\mathrm{old}} 
      \left(\frac{E}{D}\right)^{-1/(q+1)}
      \f]
      This encompasses all the standard error scaling methods. To
      avoid uncontrolled changes in the stepsize, the overall
      scaling factor is limited to the range 1/5 to 5.
  */
  template<class vec_t> class gsl_ode_control {

  protected:

    /// Number of scalings
    size_t sdim;
    
    /// Scalings
    double *scale_abs;

  public:

    /// Absolute precision (default \f$ 10^{-6} \f$)
    double eps_abs;
    /// Relative precision (default 0)
    double eps_rel;
    /// Function scaling factor (default 1)
    double a_y;
    /// Derivative scaling factor (default 0)
    double a_dydt;
    /// Use standard or scaled algorithm (default true)
    bool standard;

    gsl_ode_control() {
      eps_abs=1.0e-6;
      eps_rel=0.0;
      a_y=1.0;
      a_dydt=0.0;
      standard=true;

      sdim=0;
    }
    
    virtual ~gsl_ode_control() {
      if (sdim!=0) delete[] scale_abs;
    }      

    /** 
        \brief Set the scaling for each differential equation
    */
    template<class svec_t> int set_scale(size_t nscal, const svec_t &scale) {
      if (nscal>=sdim) {
        if (sdim!=0) delete[] scale_abs;
        scale_abs=new double[nscal];
        sdim=nscal;
      }
      for(size_t i=0;i<nscal;i++) {
        scale_abs[i]=scale[i];
      }
      return 0;
    }

    /// Automatically adjust step-size
    virtual int hadjust(size_t dim, unsigned int ord, const vec_t &y,
                        vec_t &yerr, vec_t &yp, double *h) {
      
      const double S = 0.9;
      const double h_old = *h;
      
      double rmax = DBL_MIN;
      size_t i;
        
      for(i=0; i<dim; i++) {
        double D0;
        if (standard || sdim==0) {
          D0=eps_rel*(a_y*fabs(y[i])+a_dydt*fabs(h_old*yp[i]))+eps_abs;
        } else {
          D0=eps_rel*(a_y*fabs(y[i])+a_dydt*fabs(h_old*yp[i]))+
            eps_abs*scale_abs[i%sdim];
        }
        const double r=fabs(yerr[i]) / fabs(D0);
        rmax = GSL_MAX_DBL(r, rmax);
      }
        
      if(rmax > 1.1) {

        /* decrease step, no more than factor of 5, but a fraction S more
           than scaling suggests (for better accuracy) */

        double r =  S / pow(rmax, 1.0/ord);
          
        if (r < 0.2) {
          r = 0.2;
        }
          
        *h = r * h_old;
          
        return GSL_ODEIV_HADJ_DEC;

      } else if (rmax < 0.5) {

        /* increase step, no more than factor of 5 */
        double r = S / pow(rmax, 1.0/(ord+1.0));
          
        if (r > 5.0) {
          r = 5.0;
        }
          
        /* don't allow any decrease caused by S<1 */
        if (r < 1.0) {
          r = 1.0;
        }
          
        *h = r * h_old;

        return GSL_ODEIV_HADJ_INC;

      } else {

        /* no change */
        return GSL_ODEIV_HADJ_NIL;

      }

    }

  };

  /** 
      \brief Adaptive ODE stepper (GSL)
      
      To modify the ODE stepper which is used, use the
      adapt_step::set_step().

      Default template arguments
      - \c param_t - \ref multi_funct
      - \c func_t - \ref ode_funct
      - \c vec_t - \ref ovector_view
      - \c alloc_vec_t - \ref ovector
      - \c alloc_t - \ref ovector_alloc

      \future Fix so that memory allocation/deallocation is performed only
      as necessary
      \future Allow user to find out how many steps were taken, etc.

  */
  template<class param_t, class func_t=ode_funct<param_t>, 
    class vec_t=ovector_view, class alloc_vec_t=ovector, 
    class alloc_t=ovector_alloc> class gsl_astep : 
  public adapt_step<param_t,func_t,vec_t,alloc_vec_t,alloc_t> {
    
#ifndef DOXYGEN_INTERNAL
    
    protected:
    
    /// Memory allocator for objects of type \c alloc_vec_t
    alloc_t ao;
    
    /// Temporary storage for yout
    alloc_vec_t yout;

    /// Temporary storage for dydx_out
    alloc_vec_t dydx_out;
    
    /// The size of the last step
    double last_step;
    /// The number of steps (?)
    unsigned long int count;
    /// The number of failed steps
    unsigned long int failed_steps;

    /// Apply the evolution for the next adaptive step
    int evolve_apply(double &t, double &h, double t1, size_t nvar,
                     vec_t &y, vec_t &dydx, vec_t &yout2, vec_t &yerr,
                     vec_t &dydx_out2, param_t &pa, func_t &derivs) {
      
      double t0 = t;
      double h0 = h;
      int step_status;
      int final_step = 0;
      /* remaining time, possibly less than h */
      double dt = t1 - t0;  
        
      if ((dt < 0.0 && h0 > 0.0) || (dt > 0.0 && h0 < 0.0)) {
        O2SCL_ERR("Step direction must match interval direction", 
                gsl_einval);
      }
      
      bool try_step_again=true;
      while(try_step_again) {

        if ((dt >= 0.0 && h0 > dt) || (dt < 0.0 && h0 < dt)) {
          h0 = dt;
          final_step = 1;
        } else{
          final_step = 0;
        }
      
        step_status=this->stepp->step(t0,h0,nvar,y,dydx,yout2,yerr,
                                      dydx_out2,pa,derivs);
      
        /* Check for stepper internal failure */
        if (step_status != gsl_success) {
          // Notify user which step-size caused the failure
          h=h0;
          return step_status;
        }
  
        count++;
        last_step = h0;
        if (final_step) {
          t = t1;
        } else{
          t = t0 + h0;
        }
  
        /* Check error and attempt to adjust the step. */
        const int hadjust_status=con.hadjust(nvar,this->stepp->get_order(),
                                             yout2,yerr,dydx_out2,&h0);
        
        if (hadjust_status == GSL_ODEIV_HADJ_DEC) {
          
          /* Step was decreased. Undo and go back to try again. */
          failed_steps++;
          try_step_again=true;
        } else {
          try_step_again=false;
        }

      }
        
      /* suggest step size for next time-step */
      h = h0;  
  
      return step_status;
    }
      
#endif
    
    public:
      
    /// Control specification
    gsl_ode_control<vec_t> con;
    
    gsl_astep() {
      this->verbose=0;
    }
      
    virtual ~gsl_astep() {
    }
    
    /** 
        \brief Make an adaptive integration step of the system 
        \c derivs with derivatives
          
        This attempts to take a step of size \c h from the point \c
        x of an \c n-dimensional system \c derivs starting with \c y
        and given the initial derivatives \c dydx. On exit, \c x, \c
        y and \c dydx contain the new values at the end of the step,
        \c h contains the size of the step, \c dydx
        contains the derivative at the end of the step, and \c yerr
        contains the estimated error at the end of the step.
    */
    virtual int astep_derivs(double &x, double &h, double xmax,
                             size_t n, vec_t &y, vec_t &dydx, 
                             vec_t &yerr, param_t &pa, 
                             func_t &derivs) {
      count=0;
      failed_steps=0;
      last_step=0.0;

      ao.allocate(yout,n);
      ao.allocate(dydx_out,n);

      int ret=evolve_apply(x,h,xmax,n,y,dydx,yout,yerr,dydx_out,
                           pa,derivs);

      for(size_t i=0;i<n;i++) {
        y[i]=yout[i];
        dydx[i]=dydx_out[i];
      }

      if (this->verbose>0) {
        std::cout << x << " ";
        for(size_t j=0;j<n;j++) std::cout << y[j] << " ";
        std::cout << std::endl;
      }

      ao.free(yout);
      ao.free(dydx_out);

      return ret;
    }
      
    /** 
        \brief Make an adaptive integration step of the system 
        \c derivs 
          
        This attempts to take a step of size \c h from the point \c
        x of an \c n-dimensional system \c derivs starting with \c
        y. On exit, \c x and \c y contain the new values at the end
        of the step, \c h contains the size of the step, \c dydx_out
        contains the derivative at the end of the step, and \c yerr
        contains the estimated error at the end of the step.
    */
    virtual int astep(double &x, double &h, double xmax,
                      size_t n, vec_t &y, vec_t &u_dydx_out,
                      vec_t &yerr, param_t &pa, 
                      func_t &derivs) {
      count=0;
      failed_steps=0;
      last_step=0.0;

      alloc_vec_t dydx;
      ao.allocate(yout,n);
      ao.allocate(dydx,n);
      
      int ret=derivs(x,n,y,dydx,pa);
      if (ret!=0) {
        O2SCL_ADD_ERR_RET("Initial derivative evalution failed in astep().",
                    gsl_efailed);
      }
      
      ret=evolve_apply(x,h,xmax,n,y,dydx,yout,yerr,u_dydx_out,
                       pa,derivs);
      
      for(size_t i=0;i<n;i++) {
        y[i]=yout[i];
      }

      if (this->verbose>0) {
        std::cout << x << " ";
        for(size_t j=0;j<n;j++) std::cout << y[j] << " ";
        std::cout << std::endl;
      }

      ao.free(yout);
      ao.free(dydx);

      return ret;
    }
      
  };
  
#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