c++11 - 使用 boost odeint 的精确多维积分

Question

使用 boost odeint 以高精度计算多维积分的推荐方法是什么？以下代码将 f=x*y 从 -1 积分到 2，但相对于解析解的误差超过 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;
    }

印刷：

    Exact: 2.25, numerical: 2.19088

我应该在内部积分中使用更严格的停止条件还是有更快/更准确的方法来做到这一点？谢谢！

Answer 1

首先，可能有比 ODE 方案（http://en.wikipedia.org/wiki/Numerical_integration）更好的数值方法来计算高维积分，但我认为这是 odeint 的一个巧妙应用，我发现。

但是，您的代码的问题在于它假定您在外部积分中使用观察者来更新内部积分的 x 值。但是，那integrate函数内部使用密集输出步进器，这意味着实际时间步长和观察者调用不同步。所以内部积分的 x 不会在正确的时刻更新。一个快速的解决方法是使用带有 runge_kutta4 步进器的集成常量，它使用恒定步长并确保在外循环的每个步骤之后调用观察者调用，从而调用 x-updates。然而，这有点依赖于集成例程的一些内部细节。更好的方法是以这样一种方式设计程序，即状态确实是二维的，但每个积分仅适用于两个变量之一。