【C++】ガウス求積

# ガウス求積とは

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

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

詳しくは[Wikipedia](https://ja.wikipedia.org/wiki/%E3%82%AC%E3%82%A6%E3%82%B9%E6%B1%82%E7%A9%8D)等を参照して下さい。多変数関数（$f:\mathbb{R}^3 \rightarrow \mathbb{R}$など）やベクトル関数（$f:\mathbb{R}^3 \rightarrow \mathbb{R}^3$など）にも拡張できます。

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

# C++での実装

## 実関数の場合

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

```c++
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$を計算する静的関数を宣言し、デフォルトコンストラクタにおいてメンバ変数をデフォルト値で初期化するように定義します。

```c++
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`定義します。

```c++
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`クラスを使うと、ガウス求積は以下のように実行することができます。

```c++
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ライブラリを使って行います。

```c++
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;
  }

  // ...
};
```

## まとめ

<details>
<summary>コード全体</summary><div>

```c++
#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
};
```
</div></details>