ガウス求積とは

実数の標準化された閉区間$[-1, 1]$で定義された実関数$f:\mathbb{R} \rightarrow \mathbb{R}$の定積分を数値的に求める手法です。

I = \int_{-1}^{1} f(x)\, \mathrm{d}x = \sum_{i=1}^N w_i f(x_i)

詳しくはWikipedia等を参照して下さい。多変数関数（$f:\mathbb{R}^3 \rightarrow \mathbb{R}$など）やベクトル関数（$f:\mathbb{R}^3 \rightarrow \mathbb{R}^3$など）にも拡張できます。

今回はガウス求積のC++での実装例を紹介します。

C++での実装

実関数の場合

関数値を計算する点$x_i$とその重み$w_i$は、積分点数Nをテンプレートパラメータとするgaussian_quadratureテンプレートクラスのメンバ変数として定義します。

template <size_t N>
class gaussian_quadrature {
 public:
  /// ...

 private:
  std::array<double, N> x_;  /// Abscissas
  std::array<double, N> w_;  /// Weights
};

各Nにおけるデフォルトの$x_i$および$w_i$を計算する静的関数を宣言し、デフォルトコンストラクタにおいてメンバ変数をデフォルト値で初期化するように定義します。

template <typename size_t N>
class gaussian_quadrature {
 public:
  static const std::array<double, N>& default_abscissas() noexcept;
  static const std::array<double, N>& default_weights() noexcept;

  gaussian_quadrature() : x_{default_abscissas()}, w_{default_weights()} {}

  // ...
};

default_abscissasおよびdefault_weights関数の実装は、ソースファイルで定義します。

積分を実行する関数を、関数$f$をパラメータとするメンバ関数テンプレートとして定義します。$f:\mathbb{R} \rightarrow \mathbb{R}$の場合の関数をintegrate_1dとして、$f:\mathbb{R}^2 \rightarrow \mathbb{R}$の場合の関数をintegrate_2d定義します。

template <typename size_t N>
class gaussian_quadrature {
 public:
  // ...

  /// 積分
  template <typename F>
  double integrate_1d(F&& f) {
    double I = 0;
    for (size_t i = 0; i < N; ++i) {
      I += w_[i] * f(x_[i]);
    }
    return I;
  }

  /// 二重積分
  template <typename F>
  double integrate_2d(F&& f) {
    double I = 0;
    for (size_t i = 0; i < N; ++i) {
      for (size_t j = 0; j < N; ++j {
        I += w_[i] * w_[j] * f(x_[i], x_[j]);
      }
    }
    return I;
  }

  // ...
};

このgaussian_quadratureクラスを使うと、ガウス求積は以下のように実行することができます。

gaussian_quadrature<5> integrator;

const auto I1 = integrator.integrate_1d([](double x) { return x + 1; });
assert(std::abs(I1 - 2.0) < 1e-6);

const auto I2 = integrator.integrate_2d([](double x, double y) { return x + y + 1; });
assert(std::abs(I2 - 4.0) < 1e-6);

ベクトル関数の場合

例えば２次元ベクトル場の面積分を求める場合、上のintegrate_2dの関数パラメータFの戻り値による関数オーバーロード（戻り値がスカラーかベクトルか）が必要となります。C++では戻り値によるオーバーロードができないので、SFINAEを使って実現します。関数の戻り値はC++17ならばstd::invoke_result_tを使って得ることが可能です。

ここでは、ベクトル計算はEigenライブラリを使って行います。

template <typename size_t N>
class gaussian_quadrature {
 public:
  // ...

  /// 二重積分、実関数バージョン
  /// std::invoke_result_tはC++17から使用可能
  template <typename F>
  auto integrate_2d(F&& f) -> std::enable_if_t<
      std::is_same_v<std::invoke_result_t<F, double, double>, double>,
      double> {
    double I = 0;
    for (size_t i = 0; i < N; ++i) {
      for (size_t j = 0; j < N; ++j {
        I += w_[i] * w_[j] * f(x_[i], x_[j]);
      }
    }
    return I;
  }

  /// 二重積分、ベクトル関数バージョン
  template <typename F>
  auto integrate_2d(F&& f) -> std::enable_if_t<
      std::is_same_v<std::invoke_result_t<F, double, double>, Eigen::Vector2d>,
      Eigen::Vector2d> {
    Eigen::Vector2d I = Eigen::Vector2d::Zero();
    for (size_t i = 0; i < N; ++i) {
      for (size_t j = 0; j < N; ++j {
        I += w_[i] * w_[j] * f(x_[i], x_[j]);
      }
    }
    return I;
  }

  // ...
};

まとめ

コード全体

# include <Eigen/Core>
# include <array>
# include <type_traits>

template <size_t N>
class gaussian_quadrature {
 public:
  static const std::array<double, N>& default_abscissas() noexcept;
  static const std::array<double, N>& default_weights() noexcept;

  gaussian_quadrature() : x_{default_abscissas()}, w_{default_weights()} {}

  /// 積分
  template <typename F>
  double integrate_1d(F&& f) {
    double I = 0;
    for (size_t i = 0; i < N; ++i) {
      I += w_[i] * f(x_[i]);
    }
    return I;
  }

  /// 二重積分、実関数バージョン
  /// std::invoke_result_tはC++17から使用可能
  template <typename F>
  auto integrate_2d(F&& f) -> std::enable_if_t<
      std::is_same_v<std::invoke_result_t<F, double, double>, double>,
      double> {
    double I = 0;
    for (size_t i = 0; i < N; ++i) {
      for (size_t j = 0; j < N; ++j {
        I += w_[i] * w_[j] * f(x_[i], x_[j]);
      }
    }
    return I;
  }

  /// 二重積分、ベクトル関数バージョン
  template <typename F>
  auto integrate_2d(F&& f) -> std::enable_if_t<
      std::is_same_v<std::invoke_result_t<F, double, double>, Eigen::Vector2d>,
      Eigen::Vector2d> {
    Eigen::Vector2d I = Eigen::Vector2d::Zero();
    for (size_t i = 0; i < N; ++i) {
      for (size_t j = 0; j < N; ++j {
        I += w_[i] * w_[j] * f(x_[i], x_[j]);
      }
    }
    return I;
  }

 private:
  std::array<double, N> x_;  /// Abscissas
  std::array<double, N> w_;  /// Weights
};