\(\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_mul.cpp

View page source

Multi-Dimensional Romberg Integration: Example and Test

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


namespace {

   class TestFun {
   private:
      const CppAD::vector<size_t> deg;
   public:
      // constructor
      TestFun(const CppAD::vector<size_t> deg_)
      : deg(deg_)
      { }

      // function F(x) = x[0]^deg[0] * x[1]^deg[1]
      double operator () (const CppAD::vector<double> &x)
      {  size_t i;
         double   f = 1;
         for(i = 0; i < deg[0]; i++)
            f *= x[0];
         for(i = 0; i < deg[1]; i++)
            f *= x[1];
         return f;
      }
   };

}

bool RombergMul(void)
{  bool ok = true;
   size_t i;
   size_t k;

   CppAD::vector<size_t> deg(2);
   deg[0] = 5;
   deg[1] = 3;
   TestFun F(deg);

   CppAD::RombergMul<
      TestFun              ,
      CppAD::vector<size_t>,
      CppAD::vector<double>,
      2                    > RombergMulTest;

   // arugments to RombergMul
   CppAD::vector<double> a(2);
   CppAD::vector<double> b(2);
   CppAD::vector<size_t> n(2);
   CppAD::vector<size_t> p(2);
   for(i = 0; i < 2; i++)
   {  a[i] = 0.;
      b[i] = 1.;
   }
   n[0] = 4;
   n[1] = 3;
   double r, e;

   // int_a1^b1 dx1 int_a0^b0 F(x0,x1) dx0
   //    = [ b0^(deg[0]+1) - a0^(deg[0]+1) ] / (deg[0]+1)
   //    * [ b1^(deg[1]+1) - a1^(deg[1]+1) ] / (deg[1]+1)
   double bpow = 1.;
   double apow = 1.;
   for(i = 0; i <= deg[0]; i++)
   {  bpow *= b[0];
      apow *= a[0];
   }
   double check = (bpow - apow) / double(deg[0]+1);
   bpow = 1.;
   apow = 1.;
   for(i = 0; i <= deg[1]; i++)
   {  bpow *= b[1];
      apow *= a[1];
   }
   check *= (bpow - apow) / double(deg[1]+1);

   double step = (b[1] - a[1]) / exp(log(2.)*double(n[1]-1));
   double spow = 1;
   for(k = 0; k <= n[1]; k++)
   {  spow = spow * step * step;
      double bnd = 3 * double(deg[1] + 1) * spow;

      for(i = 0; i < 2; i++)
         p[i] = k;
      r    = RombergMulTest(F, a, b, n, p, e);

      ok  &= e < bnd;
      ok  &= CppAD::NearEqual(check, r, 0., e);

   }

   return ok;
}