Boost C++ Libraries

...one of the most highly regarded and expertly designed C++ library projects in the world. Herb Sutter and Andrei Alexandrescu, C++ Coding Standards

This is the documentation for an old version of Boost. Click here to view this page for the latest version.
PrevUpHomeNext
Calculating an Integral

Similar to the generic derivative example, we can calculate integrals in a similar manner:

template<typename value_type, typename function_type>
inline value_type integral(const value_type a,
                           const value_type b,
                           const value_type tol,
                           function_type func)
{
   unsigned n = 1U;

   value_type h = (b - a);
   value_type I = (func(a) + func(b)) * (h / 2);

   for(unsigned k = 0U; k < 8U; k++)
   {
      h /= 2;

      value_type sum(0);
      for(unsigned j = 1U; j <= n; j++)
      {
         sum += func(a + (value_type((j * 2) - 1) * h));
      }

      const value_type I0 = I;
      I = (I / 2) + (h * sum);

      const value_type ratio     = I0 / I;
      const value_type delta     = ratio - 1;
      const value_type delta_abs = ((delta < 0) ? -delta : delta);

      if((k > 1U) && (delta_abs < tol))
      {
         break;
      }

      n *= 2U;
   }

   return I;
}

The following sample program shows how the function can be called, we begin by defining a function object, which when integrated should yield the Bessel J function:

template<typename value_type>
class cyl_bessel_j_integral_rep
{
public:
   cyl_bessel_j_integral_rep(const unsigned N,
      const value_type& X) : n(N), x(X) { }

   value_type operator()(const value_type& t) const
   {
      // pi * Jn(x) = Int_0^pi [cos(x * sin(t) - n*t) dt]
      return cos(x * sin(t) - (n * t));
   }

private:
   const unsigned n;
   const value_type x;
};

   /* The function can now be called as follows: */
int main(int, char**)
{
   using boost::math::constants::pi;
   typedef boost::multiprecision::cpp_dec_float_50 mp_type;

   const float j2_f =
      integral(0.0F,
      pi<float>(),
      0.01F,
      cyl_bessel_j_integral_rep<float>(2U, 1.23F)) / pi<float>();

   const double j2_d =
      integral(0.0,
      pi<double>(),
      0.0001,
      cyl_bessel_j_integral_rep<double>(2U, 1.23)) / pi<double>();

   const mp_type j2_mp =
      integral(mp_type(0),
      pi<mp_type>(),
      mp_type(1.0E-20),
      cyl_bessel_j_integral_rep<mp_type>(2U, mp_type(123) / 100)) / pi<mp_type>();

   // 0.166369
   std::cout
      << std::setprecision(std::numeric_limits<float>::digits10)
      << j2_f
      << std::endl;

   // 0.166369383786814
   std::cout
      << std::setprecision(std::numeric_limits<double>::digits10)
      << j2_d
      << std::endl;

   // 0.16636938378681407351267852431513159437103348245333
   std::cout
      << std::setprecision(std::numeric_limits<mp_type>::digits10)
      << j2_mp
      << std::endl;

   //
   // Print true value for comparison:
   // 0.166369383786814073512678524315131594371033482453329
   std::cout << boost::math::cyl_bessel_j(2, mp_type(123) / 100) << std::endl;
}


PrevUpHomeNext