3

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

内部積分でより厳しい停止条件を使用する必要がありますか、これを行うためのより高速で正確な方法はありますか? ありがとう！

4

1 に答える 1

2

まず、ODE スキーム ( http://en.wikipedia.org/wiki/Numerical_integration )よりも高次元の積分を計算するためのより優れた数値的方法があるかもしれませんが、これは odeint の適切なアプリケーションだと思います。

ただし、コードの問題は、外側の統合でオブザーバーを使用して内側の統合の x 値を更新することを前提としていることです。しかしintegrate関数は高密度出力ステッパーを内部で使用します。これは、実際の時間ステップとオブザーバー呼び出しが同期していないことを意味します。したがって、内部統合の x は適切なタイミングで更新されません。簡単な修正は、定数のステップ サイズを使用し、外側のループの各ステップの後にオブザーバー呼び出し、つまり x-updates が呼び出されるようにする runge_kutta4 ステッパーで integrate_const を使用することです。ただし、これは統合ルーチンの内部の詳細に依存するちょっとしたハックです。より良い方法は、状態が実際に 2 次元になるようにプログラムを設計することですが、各統合は 2 つの変数のうちの 1 つに対してのみ機能します。

于 2014-04-03T14:37:02.627 に答える