demo.cpp

   1 #include <cmath>
   2 #include <iostream>
   3
   4 #include "cubaint.hpp"
   5
   6 int test_integrand (const int *ndim, const double x[],
   7                     const int *ncomp, double f[], void *userdata)
   8 {
   9   for (int icomp=0; icomp<*ncomp; icomp++) {
  10     f[icomp] = 1;
  11     for (int idim=0; idim<*ndim; idim++)
  12       f[icomp] *= exp(-pow(x[idim],2));
  13   }
  14
  15   return 0;
  16 }
  17
  18 int test_integrand_c (const int *ndim, const double x[],
  19                       const int *ncomp, std::complex<double> f[], void *userdata)
  20 {
  21   for (int icomp=0; icomp<*ncomp; icomp++) {
  22     f[icomp] = 1.;
  23     for (int idim=0; idim<*ndim; idim++)
  24       f[icomp] *= exp(-pow(x[idim],2));
  25   }
  26
  27   return 0;
  28 }
  29
  30 int main ()
  31 {
  32   int nregions, neval, fail;
  33   double error, prob, integral;
  34
  35   std::cout << "int_0^1 dx e-x^2: " << std::endl;
  36   for (int itype=0; itype<4; itype++) {
  37     cubaint::integrate (&test_integrand, cubaint::intmethod(itype),
  38                         {{0,1}}, 1,
  39                         nregions,
  40                         neval, fail, &integral,
  41                         &error, &prob, 0);
  42     std::cout << integral << " +- " << error << " (" << itype << ")" << std::endl;
  43   }
  44
  45   std::cout << "int_0^INF dx e-x^2: " << std::endl;
  46   for (int itype=0; itype<4; itype++) {
  47     cubaint::integrate (&test_integrand, cubaint::intmethod(itype),
  48                         {{0,cubaint::INF}}, 1,
  49                         nregions,
  50                         neval, fail, &integral,
  51                         &error, &prob, 0);
  52     std::cout << integral << " +- " << error << " (" << itype << ")" << std::endl;
  53   }
  54
  55   std::cout << "int_-INF^0 dx e-x^2: " << std::endl;
  56   for (int itype=0; itype<4; itype++) {
  57     cubaint::integrate (&test_integrand, cubaint::intmethod(itype),
  58                         {{cubaint::INF,0}}, 1,
  59                         nregions,
  60                         neval, fail, &integral,
  61                         &error, &prob, 0);
  62     std::cout << integral << " +- " << error << " (" << itype << ")" << std::endl;
  63   }
  64
  65   std::cout << "int_-INF^INF dx e-x^2: " << std::endl;
  66   for (int itype=0; itype<4; itype++) {
  67     cubaint::integrate (&test_integrand, cubaint::intmethod(itype),
  68                         {{cubaint::INF,cubaint::INF}}, 1,
  69                         nregions,
  70                         neval, fail, &integral,
  71                         &error, &prob, 0);
  72     std::cout << integral << " +- " << error << " (" << itype << ")" << std::endl;
  73   }
  74
  75   std::cout << "int_1^INF dx e-x^2: " << std::endl;
  76   for (int itype=0; itype<4; itype++) {
  77     cubaint::integrate (&test_integrand, cubaint::intmethod(itype),
  78                         {{1,cubaint::INF}}, 1,
  79                         nregions,
  80                         neval, fail, &integral,
  81                         &error, &prob, 0);
  82     std::cout << integral << " +- " << error << " (" << itype << ")" << std::endl;
  83   }
  84
  85   std::cout << "int_-INF^1 dx e-x^2: " << std::endl;
  86   for (int itype=0; itype<4; itype++) {
  87     cubaint::integrate (&test_integrand, cubaint::intmethod(itype),
  88                         {{cubaint::INF,1}}, 1,
  89                         nregions,
  90                         neval, fail, &integral,
  91                         &error, &prob, 0);
  92     std::cout << integral << " +- " << error << " (" << itype << ")" << std::endl;
  93   }
  94
  95   std::complex<double> error_c, prob_c, integral_c;
  96   std::cout << "int_0^1 dx (e-x^2+i0): " << std::endl;
  97   for (int itype=0; itype<4; itype++) {
  98     cubaint::integrate (&test_integrand_c, cubaint::intmethod(itype),
  99                         {{0,1}}, 1,
 100                         nregions,
 101                         neval, fail, &integral_c,
 102                         &error_c, &prob_c, 0);
 103     std::cout << integral_c << " +- " << error_c << " (" << itype << ")" << std::endl;
 104   }
 105
 106 }