ガウス求積法
前節で紹介した台形則やシンプソン則は、積分区間を等間隔に分割する手法(ニュートン・コーツの公式)でした。これに対し、ガウス求積法 (Gaussian Quadrature)、特にガウス・ルジャンドル積分 (Gauss-Legendre Integration) は、積分点(分点)の配置と重みを最適化することで、より少ない点数で極めて高い精度を実現する手法です。
この手法は、被積分関数が滑らかである(高階微分可能である)場合に圧倒的な威力を発揮します。
原理
基本公式
ガウス・ルジャンドル積分では、積分区間を に正規化した上で、以下の近似式を用います。
ここで、 は積分点の数(次数)、 は積分点、 は重み係数です。 ガウス求積法の最大の特徴は、 個の点を用いるだけで、 次以下の任意の多項式を厳密に積分できる という点です。
例えば、 の場合、5次関数までの積分が厳密に一致します。シンプソン則(3点等間隔)が3次関数までしか厳密でないことと比較すると、その効率の良さが分かります。
積分点と重み
積分点 は、ルジャンドル多項式 (Legendre Polynomials) の零点(根)として選ばれます。 また、重み は以下の式で決定されます。
これらは解析的に計算可能であり、多くの数値計算ライブラリではあらかじめ計算された値(テーブル)が用いられます。
例: n=5 の場合の積分点と重み
| i | x_i (積分点) | w_i (重み) |
|---|---|---|
| 1, 5 | ± 0.9061798459… | 0.2369268850… |
| 2, 4 | ± 0.5384693101… | 0.4786286704… |
| 3 | 0.0000000000… | 0.5688888888… |
積分点は原点に対して対称に配置され、端点( )よりも内側に配置されるのが特徴です。
積分区間の変換
一般的な積分区間 での積分計算を行うためには、変数変換を行って積分区間を に写像する必要があります。
変数変換:
これにより、公式は以下のように変形されます。
Rustによる実装
ここでは、
の場合のガウス・ルジャンドル積分を実装します。積分点と重みは定数としてハードコードします。
より高い次数が必要な場合は、ルジャンドル多項式の根を数値的に求めるアルゴリズム(ニュートン法など)を実装するか、gauss-quad のような外部クレートを利用するのが一般的です。
/// ガウス・ルジャンドル積分(5点)
///
/// 5点のガウス求積法を用いて、区間 [a, b] で関数 f(x) を積分します。
/// 9次以下の多項式に対して厳密解を与えます。
fn gauss_legendre_5<F>(f: F, a: f64, b: f64) -> f64
where
F: Fn(f64) -> f64,
{
// n=5 の積分点と重み
// 対称性を利用して定義することも可能ですが、ここでは列挙します
const NODES: [f64; 5] = [
0.0,
0.5384693101056831,
-0.5384693101056831,
0.9061798459386640,
-0.9061798459386640,
];
const WEIGHTS: [f64; 5] = [
0.5688888888888889,
0.4786286704993665,
0.4786286704993665,
0.2369268850561891,
0.2369268850561891,
];
let mid = 0.5 * (a + b);
let half_len = 0.5 * (b - a);
let mut sum = 0.0;
for i in 0..5 {
let x = mid + half_len * NODES[i];
sum += WEIGHTS[i] * f(x);
}
sum * half_len
}
fn main() {
// 積分対象: f(x) = sin(x)
// 区間: [0, PI]
// 解析解: 2.0
use std::f64::consts::PI;
let f = |x: f64| x.sin();
let a = 0.0;
let b = PI;
let exact = 2.0;
let result = gauss_legendre_5(f, a, b);
println!("Gaussian Quadrature (n=5)");
println!("Result: {:.16}", result);
println!("Error: {:.2e}", (result - exact).abs());
// 比較: シンプソン則 (n=10, 評価点数は11点)
// ガウス積分(n=5)は評価点数が5点なので、それより多くの点数を使う条件で比較
println!("\n(参考) Simpson's Rule (n=10)");
// シンプソン則の実装は省略(前節参照)し、結果のオーダーのみ記述
// 前節の結果から、n=10のとき誤差は約 6.8e-6 程度
}
実行結果(例):
Gaussian Quadrature (n=5)
Result: 2.0000001102844727
Error: 1.10e-7
(参考) Simpson's Rule (n=10)
Error: ~ 6.8e-6
たった5回の関数評価で、 という高い精度が得られていることがわかります。シンプソン則で同等の精度を得ようとすると、分割数をかなり増やす必要があります。
長所と短所
長所
- 圧倒的な効率: 滑らかな関数に対しては、非常に少ない評価点数で高精度が得られる。
- 計算コストの削減: 関数評価( の計算)に時間がかかる場合、評価回数を減らせるメリットが大きい。
短所
- 汎用性の低さ: 積分点 と重み が事前に必要であり、等間隔データ(実験データなど)には適用できない。
- 特異点への弱さ: 端点に特異点がある場合や、関数が滑らかでない(不連続点や折れ曲がりがある)場合、精度が急激に低下する。
- 計算の複雑さ: 任意の に対して積分点と重みを求めるには、特殊関数の計算が必要となる。
発展:分点の自動生成
任意の次数 に対する積分点と重みは、ルジャンドル多項式 の漸化式と、ニュートン法による求根アルゴリズムを組み合わせることで計算できます。
ルジャンドル多項式の漸化式:
また、導関数の関係式:
これらを用いて、
から
の間で
となる解を探すプログラムを作成すれば、任意の次数のガウス求積法を実装可能です。Rustのエコシステムには gauss-quad などのクレートがあり、これらを利用するのも良い選択肢です。
まとめ
- ガウス求積法は、積分点を最適化することで、滑らかな関数の積分において最高の効率を誇る。
- 点で 次の多項式を厳密に積分できる。
- 実験データではなく、数式として与えられた関数の積分(特に物理シミュレーションにおける相互作用計算など)に適している。
次節では、より柔軟に積分精度を制御できる「適応型積分」について学びます。