Integrator.h

Go to the documentation of this file.

#ifndef _utl_Integrator_h_
#define _utl_Integrator_h_

#include <utl/PolynomialInterpolation.h>
#include <utl/ErrorLogger.h>

#include <vector>


namespace utl {

  template<class Functor>
  class Integrator {

  public:
    Integrator(Functor& functor, const double accuracy = 1e-5)
      : fFunctor(functor), fAccuracy(accuracy) { }

    void SetAccuracy(const double accuracy) { fAccuracy = accuracy; }

    double GetIntegral(const double a, const double b) const
    { return GetRombergIntegral(a, b); }

    double
    GetRombergIntegral(const double a, const double b,
                       const int order = 5, const int maxIterations = 20)
      const
    {
      const double delta = b - a;

      std::vector<double> tInt(maxIterations+1);
      std::vector<double> tStep(maxIterations+1);
      double oldInt = tInt[0] = GetBoxAverage(a, delta, 2);
      tStep[0] = delta;
      int i;
      for (i = 1; i < order; ++i) {
        oldInt = tInt[i] = GetTrapezoidalAverage(oldInt, a, delta, i);
        tStep[i] = delta / (1 << (2*i));  // h^2
      }

      const double eps = 0.5*fAccuracy;
      double extrapolation = 1e10;
      for ( ; i <= maxIterations; ++i) {
        oldInt = tInt[i] = GetTrapezoidalAverage(oldInt, a, delta, i);
        tStep[i] = delta / (1 << (2*i));  // extrapolate h^2 dependence for h->0
        const int first = i - order + 1;
        double residual = 0;
        extrapolation =
          PolynomialInterpolation(order, &tStep[first], &tInt[first], 0, residual);
        if (fabs(residual) <= eps*fabs(extrapolation))
          return delta*extrapolation;
      }
      if (fabs(extrapolation) > fAccuracy)
        WARNING("Maximal level limit reached, result is not accurate.");
      return delta*extrapolation;
    }

    double
    GetSimpsonIntegral(const double a, const double b,
                       const int minLevel = 4, const int maxLevel = 20)
      const
    {
      const double delta = b - a;
      double trapezAverage = GetBoxAverage(a, delta, 2);

      int level;
      for (level = 1; level < minLevel; ++level)
        trapezAverage = GetTrapezoidalAverage(trapezAverage, a, delta, level);

      double trapezOld = trapezAverage;
      trapezAverage = GetTrapezoidalAverage(trapezAverage, a, delta, level++);
      double simpsonOld = (4./3)*trapezAverage - (1./3)*trapezOld;

      for ( ; level <= maxLevel; ++level) {
        trapezOld = trapezAverage;
        trapezAverage = GetTrapezoidalAverage(trapezAverage, a, delta, level);
        const double simpsonAverage = (4./3)*trapezAverage - (1./3)*trapezOld;
        if (fabs(simpsonAverage - simpsonOld) < fAccuracy*fabs(simpsonAverage) ||
            (!simpsonAverage && !simpsonOld))
          return delta*simpsonAverage;
        simpsonOld = simpsonAverage;
      }
      if (fabs(simpsonOld) > fAccuracy)
        WARNING("Maximal level limit reached, result is not accurate.");
      return delta*simpsonOld;
    }

    double
    GetTrapezoidalIntegral(const double a, const double b,
                           const int minLevel = 4, const int maxLevel = 20)
      const
    {
      const double delta = b - a;
      double average = GetBoxAverage(a, delta, 2);
      int level;
      for (level = 1; level <= minLevel; ++level)
        average = GetTrapezoidalAverage(average, a, delta, level);
      for ( ; level <= maxLevel; ++level) {
        const double old = average;
        average = GetTrapezoidalAverage(average, a, delta, level);
        if (fabs(average - old) < fAccuracy*fabs(average) || (!average && !old))
          return delta*average;
      }
      if (fabs(average) > fAccuracy)
        WARNING("Maximal level limit reached, result is not accurate.");
      return delta*average;
    }

  private:
    double
    GetTrapezoidalAverage(const double previousApproximation,
                          const double a, const double delta, const int level)
      const
    {
      const int n = 1 << (level - 1);
      const double step = delta/n;
      return 0.5*(previousApproximation +
                  GetBoxAverage(a+0.5*step, step, n));
    }

    double
    GetBoxAverage(const double start, const double step, const int n)
      const
    {
      double sum = 0;
      for (int i = 0; i < n; ++i) {
        const double x = start + i*step;
        sum += fFunctor(x);
      }
      return sum/n;
    }

    Functor& fFunctor;
    double fAccuracy;

  };


  template<class Functor>
  inline
  Integrator<Functor>
  MakeIntegrator(Functor& f)
  {
    return Integrator<Functor>(f);
  }


  template<class Functor, unsigned int dim>
  class VectorIntegrator {

  public:
    VectorIntegrator(Functor& functor, const double accuracy = 1e-5)
      : fFunctor(functor), fAccuracy(accuracy) { }

    void SetAccuracy(const double accuracy) { fAccuracy = accuracy; }

    void
    GetIntegral(double* const result, const double a, const double b,
                const int order = 5, const int maxIterations = 20)
      const
    {
      for (unsigned int k = 0; k < dim; ++k)
        result[k] = 0;

      const double delta = b - a;

      std::vector<double> tInt[dim];
      for (unsigned int k = 0; k < dim; ++k)
        tInt[k].resize(maxIterations+1);
      std::vector<double> tStep(maxIterations+1);
      double oldInt[dim];
      GetBoxAverage(oldInt, a, delta, 2);
      for (unsigned int k = 0; k < dim; ++k)
        tInt[k][0] = oldInt[k];
      tStep[0] = delta;
      int i;
      for (i = 1; i < order; ++i) {
        GetTrapezoidalAverage(oldInt, a, delta, i);
        for (unsigned int k = 0; k < dim; ++k)
          tInt[k][i] = oldInt[k];
        tStep[i] = delta / (1 << (2*i));  // h^2
      }

      const double eps = 0.5*fAccuracy;
      for (unsigned int k = 0; k < dim; ++k)
        result[k] = 0;
      for ( ; i <= maxIterations; ++i) {
        GetTrapezoidalAverage(oldInt, a, delta, i);
        for (unsigned int k = 0; k < dim; ++k)
          tInt[k][i] = oldInt[k];
        tStep[i] = delta / (1 << (2*i));  // extrapolate h^2 dependence for h->0
        const int first = i - order + 1;

        unsigned int nAccepted = 0;
        for (unsigned int k = 0; k < dim; ++k) {
          double residual = 0;
          result[k] =
            PolynomialInterpolation(order, &tStep[first], &tInt[k][first], 0, residual);
          if (fabs(residual) <= eps*fabs(result[k]))
            ++nAccepted;
        }
        if (nAccepted == dim) {
          for (unsigned int k = 0; k < dim; ++k)
            result[k] *= delta;
          return;
        }
      }

      WARNING("Maximal level limit reached, result is not accurate.");
      for (unsigned int k = 0; k < dim; ++k)
        result[k] *= delta;
    }

  private:
    void
    GetTrapezoidalAverage(double* const previousApproximation,
                          const double a, const double delta, const int level)
      const
    {
      const int n = 1 << (level - 1);
      const double step = delta/n;
      double add[dim];
      GetBoxAverage(add, a+0.5*step, step, n);
      for (unsigned int k = 0; k < dim; ++k)
        previousApproximation[k] = 0.5*(previousApproximation[k] + add[k]);
    }

    void
    GetBoxAverage(double* const sum, const double start, const double step, const int n)
      const
    {
      for (unsigned int k = 0; k < dim; ++k)
        sum[k] = 0;

      double result[dim];
      for (int i = 0; i < n; ++i) {
        const double x = start + i*step;
        fFunctor(result, x);
        for (unsigned int k = 0; k < dim; ++k)
          sum[k] += result[k];
      }

      for (unsigned int k = 0; k < dim; ++k)
        sum[k] /= n;
    }

    Functor& fFunctor;
    double fAccuracy;

  };

}


#endif