Accurate multidimensional integral using boost odeint

Question

What is the recommended way to calculate a multidimensional integral using boost odeint with high accuracy? The following code integrates f=x*y from -1 to 2 but the error relative to an analytic solution is over 1 % (gcc 4.8.2, -std=c++0x):

    #include "array"
    #include "boost/numeric/odeint.hpp"
    #include "iostream"
    using integral_type = std::array<double, 1>;
    int main() {
      integral_type outer_integral{0};

      double current_x = 0;
      boost::numeric::odeint::integrate(
        [&](
          const integral_type&,
          integral_type& dfdx,
          const double x
        ) {
          integral_type inner_integral{0};

          boost::numeric::odeint::integrate(
            [&current_x](
              const integral_type&,
              integral_type& dfdy,
              const double y
            ) {
              dfdy[0] = current_x * y;
            },
            inner_integral,
            -1.0,
            2.0,
            1e-3
          );

          dfdx[0] = inner_integral[0];
        },
        outer_integral,
        -1.0,
        2.0,
        1e-3,
        [&current_x](const integral_type&, const double x) {
          current_x = x; // update x in inner integrator
        }
      );
      std::cout
        << "Exact: 2.25, numerical: "
        << outer_integral[0]
        << std::endl;
      return 0;
    }

prints:

    Exact: 2.25, numerical: 2.19088

Should I just use more stringent stopping condition in the inner integrals or is there a faster/more accurate way to do this? Thanks!

Answer 1

Firstly, there might be better numerical methods to compute high-dimensional integrals than an ODE scheme (http://en.wikipedia.org/wiki/Numerical_integration), but then I think this is a neat application of odeint, I find.

However, the problem with your code is that it assumes that you use an observer in the outer integrate to update the x-value for the inner integration. However, the integrate function uses a dense-output stepper internally which means that the actual time steps and the observer calls are not in synchrony. So the x for the inner integration is not updated at the right moments. A quick fix would be to use an integrate_const with a runge_kutta4 stepper, which uses constant step size and ensure that the observer calls, and thus x-updates, are called after each step of the outer loop. However, this is a bit of a hack relying on some internal details of the integrate routines. A better way would be to design the program in such a way that the state is indeed 2-dimensional, but where each integration works only on one of the two variables.