rombergIntegrator.h

Go to the documentation of this file.

// -*- Mode: C++; tab-width: 2; -*-
// vi: set ts=2:
//
// $Id: rombergIntegrator.h,v 1.12 2003/08/26 08:04:22 oliver Exp $
//

#ifndef BALL_MATHS_ROMBERGINTEGRATOR_H
#define BALL_MATHS_ROMBERGINTEGRATOR_H

#ifndef BALL_MATHS_NUMERICALINTERGRATOR_H
# include <BALL/MATHS/numericalIntegrator.h>
#endif

namespace BALL
{
  template <typename Function, typename DataType>
  class RombergIntegrator: public NumericalIntegrator<Function, DataType>
  {
    public:

    BALL_CREATE(RombergIntegrator)

    


    RombergIntegrator(float epsilon=1E-5, Size nsteps=1000);
  
    RombergIntegrator(const RombergIntegrator& romint);

    ~RombergIntegrator();



    const RombergIntegrator& operator = (const RombergIntegrator& romint);
    
    virtual void clear();

    void setEpsilon(float eps);

    void setMaxNumSteps(Size mns);


    
    bool operator == (const RombergIntegrator& romint) const;


    
    DataType integrate(DataType from, DataType to);
    

    DataType trapezoid(DataType h, DataType from, DataType to);
    
    
    protected:

    float epsilon_;
    Size maxNumSteps_;
    vector<DataType> result_;
  };

  template<typename Function, typename DataType>
  BALL_INLINE
  RombergIntegrator<Function, DataType>::RombergIntegrator(float eps, Size nsteps): NumericalIntegrator<Function, DataType>(), epsilon_(eps), maxNumSteps_(nsteps)
  {
    result_.reserve(maxNumSteps_ / 10);
  }

  template<typename Function, typename DataType>
  BALL_INLINE
  RombergIntegrator<Function, DataType>::RombergIntegrator(const RombergIntegrator<Function, DataType>& romint):NumericalIntegrator<Function, DataType>(romint)
  {
  }

  template<typename Function, typename DataType>
  BALL_INLINE
  RombergIntegrator<Function, DataType>::~RombergIntegrator()
  {
  }

  template<typename Function, typename DataType>
  BALL_INLINE
  const RombergIntegrator<Function, DataType>&
  RombergIntegrator<Function, DataType>::operator =
  (const RombergIntegrator<Function, DataType>& romint)
  {
    function_ = romint.function_;
    maxNumSteps_ = romint.maxNumSteps_;
    epsilon_ = romint.epsilon_;
    result_ = romint.result_;
    return *this;
  }

  template<typename Function, typename DataType>
  BALL_INLINE
  void RombergIntegrator<Function, DataType>::clear() 
  {
    // Problem: function_.clear() might not exist... function_.clear();
  }

  template<typename Function, typename DataType>
  BALL_INLINE
  void RombergIntegrator<Function, DataType>::setEpsilon(float eps)
  {
    epsilon_ = eps;
  }

  template<typename Function, typename DataType>
  BALL_INLINE
  void RombergIntegrator<Function, DataType>::setMaxNumSteps(Size nsteps)
  {
    maxNumSteps_ = nsteps;
    result_.reserve(maxNumSteps_ / 10); // we hope that we do not need more than 1/10 - th of the allowed operations
  }
  
  template<typename Function, typename DataType>
  BALL_INLINE
  bool RombergIntegrator<Function, DataType>::operator ==
  (const RombergIntegrator<Function, DataType> &romint) const
    
  {
    return ((function_ == romint.function_)
        && (epsilon_  == romint.epsilon_ )
        && (maxNumSteps_ == romint.maxNumSteps_)
        && (result_      == romint.result_     ));
  }

  template<typename Function, typename DataType>
  BALL_INLINE
  DataType RombergIntegrator<Function, DataType>::trapezoid(DataType h, DataType from, DataType to)
  {
    DataType sum=0;
    DataType helper = (to - from);
    int count;
    
    Size nsteps = (Size) (sqrt((helper*helper)/(h*h)));
    for (count=1; count<nsteps-1; count++)
    {
      sum +=function_(from+(count*h));
    }

    sum+=function_(from)+function_(to);
    sum*=h;

    return sum;
  }


  template<typename Function, typename DataType>
  BALL_INLINE
  DataType RombergIntegrator<Function, DataType>::integrate
  (DataType from, DataType to)
  {
    float h = 1.;
      result_.push_back(trapezoid(h, from, to)); // this is the zeroth approximation
    
    int i=1;
    int j=0;
    int count = 0;
    DataType dev;

    do 
    {
      result_.push_back(trapezoid(h/((float) i+1), from, to));
    
      for (j=1; j <= i; j++) 
      {
        result_.push_back(result_[(i*(i+1))/2 + (j-1)] + 1. / (pow(4, j) - 1) * (result_[(i*(i+1))/2 + j-1 - result_[((i-1)*i)/2+j-1]));
        count++;
      };
      i++;
      dev = result_[((i-2)*(i-1))/2+(i-2)] - result_[((i-1)*(i))/2+(i-1)];
    } while ( (sqrt(dev*dev) > epsilon_) && (count < maxNumSteps_));

    return (result_[((i-1)*(i))/2 + (i-1)]);
  }
} // namespace BALL
  
#endif // BALL_MATHS_ROMBERGINTEGRATOR_H