\(\newcommand{\W}[1]{ \; #1 \; }\) \(\newcommand{\R}[1]{ {\rm #1} }\) \(\newcommand{\B}[1]{ {\bf #1} }\) \(\newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} }\) \(\newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} }\) \(\newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} }\) \(\newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }\)

romberg_one.cpp

View page source

One Dimensional Romberg Integration: Example and Test

# include <cppad/utility/romberg_one.hpp>
# include <cppad/utility/vector.hpp>
# include <cppad/utility/near_equal.hpp>

namespace {
   class Fun {
   private:
      const size_t degree;
   public:
      // constructor
      Fun(size_t degree_) : degree(degree_)
      { }

      // function F(x) = x^degree
      template <class Type>
      Type operator () (const Type &x)
      {  size_t i;
         Type   f = 1;
         for(i = 0; i < degree; i++)
            f *= x;
         return f;
      }
   };
}

bool RombergOne(void)
{  bool ok = true;
   size_t i;

   size_t degree = 4;
   Fun F(degree);

   // arguments to RombergOne
   double a = 0.;
   double b = 1.;
   size_t n = 4;
   size_t p;
   double r, e;

   // int_a^b F(x) dx = [ b^(degree+1) - a^(degree+1) ] / (degree+1)
   double bpow = 1.;
   double apow = 1.;
   for(i = 0; i <= degree; i++)
   {  bpow *= b;
      apow *= a;
   }
   double check = (bpow - apow) / double(degree+1);

   // step size corresponding to r
   double step = (b - a) / exp(log(2.)*double(n-1));
   // step size corresponding to error estimate
   step *= 2.;
   // step size raised to a power
   double spow = 1;

   for(p = 0; p < n; p++)
   {  spow = spow * step * step;

      r = CppAD::RombergOne(F, a, b, n, p, e);

      ok  &= e < double(degree+1) * spow;
      ok  &= CppAD::NearEqual(check, r, 0., e);
   }

   return ok;
}