ガウス求積

ガウス求積（ガウスきゅうせき、英: Gaussian quadrature）またはガウスの数値積分公式とは、カール・フリードリヒ・ガウスに因んで名づけられた数値解析における数値積分法の一種であり、実数のある閉区間(慣例的に [−1, 1] に標準化される)で定義された実数値関数のその閉区間に渡る定積分値を、比較的少ない演算で精度良く求めることができるアルゴリズムである。

n を正の整数とし、f(x) を 任意の多項式関数とする。f(x) の [−1, 1] に渡る定積分値 I を、

$I=\int _{-1}^{1}f(x)\,dx=\sum _{i=1}^{n}w_{i}f(x_{i})$

の形でなるべく正確に近似する公式を考える。ここで、x_i は積分点またはガウス点 (ガウスノード)と呼ばれる [−1, 1] 内の n 個の点であり、w_i は重みと呼ばれるn個の実数である。

実は、n 次のルジャンドル多項式の n 個の零点（これらは全て [−1, 1] 内にある）を積分点として選び、w_i を適切に選ぶと、f(x) が 2n − 1 次以下の多項式であれば上記の式が厳密に成立することが分かっている。この場合、w_i は f(x) によらず一意的に定まる。この方法を n 次のガウス・ルジャンドル (Gauss–Legendre) 公式と呼び、通常はガウス求積またはガウスの数値積分公式と言えばこの方法を指している。

f(x) が 2n − 1次を超える多項式関数の場合、または多項式関数でない場合には、上記の公式は厳密には成立しないが、f(x) が 2n − 1 次以下の多項式関数で精度よく近似できる場合には、上記の公式を f(x) に対して適用することにより、その [−1, 1] における定積分値を精度よく得ることが期待できる。それ以外の例えば、特異点のある関数の積分にはこの公式をそのまま適用することはできないが、対象の関数を f(x) = W(x) g(x) と表すことができ、g(x) が近似多項式で W(x) が既知であれば、それを代替する重み w_i を使って次のように表せる。

$\int _{-1}^{1}f(x)\,dx=\int _{-1}^{1}W(x)g(x)\,dx\approx \sum _{i=1}^{n}w_{i}g(x_{i}).$

典型的な重み関数としては、$W(x)=(1-x^{2})^{-1/2}$（ガウス–チェビシェフ）や $W(x)=e^{-x^{2}}$（ガウス–エルミート）がある。この場合の n 個の積分点 x_i はルジャンドル多項式と同様に、ある直交多項式のクラスに属する n 次多項式の根である。

ガウス・ルジャンドル公式による求積

上述のように n 次のこの方法には、n 次のルジャンドル多項式 P_n(x) が対応している。このときの n 次多項式は P_n(1) = 1 となるよう正規化され、i 番目のガウスノード x_i は i 番目の P_n の根である。重みは次の式で与えられる。

$w_{i}={\frac {2}{(\left(1-x_{i}^{2}\right)P'_{n}(x_{i}))^{2}}}.$

低次の求積法は次のようになる。

点の個数 n	点 xi	重み wi
1	0	2
2	$\pm {\sqrt {1/3}}$	1
3	0	8/9
	$\pm {\sqrt {3/5}}$	5/9
4	$\pm {\sqrt {{\Big (}3-2{\sqrt {6/5}}{\Big )}/7}}$	${\tfrac {18+{\sqrt {30}}}{36}}$
	$\pm {\sqrt {{\Big (}3+2{\sqrt {6/5}}{\Big )}/7}}$	${\tfrac {18-{\sqrt {30}}}{36}}$
5	0	128/225
	$\pm {\tfrac {1}{3}}{\sqrt {5-2{\sqrt {10/7}}}}$	${\tfrac {322+13{\sqrt {70}}}{900}}$
	$\pm {\tfrac {1}{3}}{\sqrt {5+2{\sqrt {10/7}}}}$	${\tfrac {322-13{\sqrt {70}}}{900}}$

区間変更

区間 [a, b] についての積分は、ガウス求積法を適用する前に [−1, 1] に変更する必要がある。この区間変更は以下のように行う。

$\int _{a}^{b}f(x)\,dx={\frac {b-a}{2}}\int _{-1}^{1}f\left({\frac {b-a}{2}}x+{\frac {a+b}{2}}\right)\,dx.$

ガウス求積法を適用後、以下の近似が得られる。

${\frac {b-a}{2}}\sum _{i=1}^{n}w_{i}f\left({\frac {b-a}{2}}x_{i}+{\frac {a+b}{2}}\right).$

他の形式

正の重み関数 ω を導入することで、より汎用的な積分問題の表現も可能であり、区間 [−1, 1] 以外にも適用可能である。すなわち、次の形式の問題である。

$\int _{a}^{b}\omega (x)\,f(x)\,dx.$

a, b, ω は適当に選択する。a = −1, b = 1, ω(x) = 1 のとき、前述の問題と同じ形式になる。それ以外の選択では、別の求積法になる。そのうちの一部を下記の表に示す。"A & S" という欄は、Abramowitz and Stegunにある式番号である。

区間	ω(x)	直交多項式	A & S	解説など
[−1, 1]	1	ルジャンドル多項式	25.4.29	本項（上）で解説
(−1, 1)	$(1-x)^{\alpha }(1+x)^{\beta },\quad \alpha ,\beta >-1\,$	ヤコビ多項式	25.4.33 ($\beta =0$)
(−1, 1)	${\frac {1}{\sqrt {1-x^{2}}}}$	チェビシェフ多項式（第一種）	25.4.38	チェビシェフ・ガウス求積法
[−1, 1]	${\sqrt {1-x^{2}}}$	チェビシェフ多項式（第二種）	25.4.40	チェビシェフ・ガウス求積法
[0, ∞)	$e^{-x}\,$	ラゲール多項式	25.4.45	ガウス・ラゲール求積法
(−∞, ∞)	$e^{-x^{2}}$	エルミート多項式	25.4.46	ガウス・エルミート求積法

基礎的な定理

p_n が自明でない n 次の多項式で、次のように表されるとする。

$\int _{a}^{b}\omega (x)\,x^{k}p_{n}(x)\,dx=0,\quad {\text{for all }}k=0,1,\ldots ,n-1.$

ノードとして p_n の零点を選ぶなら、全ての 2n − 1 以下の次数の多項式について正確に積分を計算できる重み w_i が存在する。さらに、それらノードは全て開区間 (a, b) にある。

この多項式 p_n は、重み関数 ω(x) に関連する次数 n の直交多項式である。

計算

ガウス求積法のノード x_i と重み w_i を計算するための基本的ツールは、直交多項式群と対応する重み関数が満たす3項漸化式である。

例えば、p_n がモニックな n 次直交多項式（最高次の項の係数が 1 の n 次直交多項式）なら、次のような漸化式で関係を表すことができる。

$p_{n+1}(x)+(B_{n}-x)p_{n}(x)+A_{n}p_{n-1}(x)=0,\qquad n=1,2,\ldots .$

このことから、対応する線型代数問題の固有値および固有ベクトルからノードと重みを計算できる。これを一般に Golub–Welsch アルゴリズムと呼ぶ。

x_i が直交多項式 p_n の根であるとき、前掲の漸化式を $k=0,1,\ldots ,n-1$ について用い、$p_{n}(x_{j})=0$ であることを踏まえると、次が成り立つことがわかる。

$J{\tilde {P}}=x_{j}{\tilde {P}}.$

ここで

${\tilde {P}}={}^{t}[p_{0}(x_{j}),p_{1}(x_{j}),...,p_{n-1}(x_{j})]$

である。そして、J はいわゆるヤコビ行列である。

$\mathbf {J} =\left({\begin{array}{llllll}B_{0}&1&0&\ldots &\ldots &\ldots \\A_{1}&B_{1}&1&0&\ldots &\ldots \\0&A_{2}&B_{2}&1&0&\ldots \\\ldots &\ldots &\ldots &\ldots &\ldots &\ldots \\\ldots &\ldots &\ldots &A_{n-2}&B_{n-2}&1\\\ldots &\ldots &\ldots &\ldots &A_{n-1}&B_{n-1}\end{array}}\right).$

したがって、ガウス求積法のノードは三重対角行列の固有値として計算できる。

重みとノードを求めるには、要素が ${\mathcal {J}}_{i,i}=J_{i,i}$, $i=1,\ldots ,n$ と ${\mathcal {J}}_{i-1,i}={\mathcal {J}}_{i,i-1}={\sqrt {J_{i,i-1}J_{i-1,i}}},\,i=2,\ldots ,n$ から成る対称な三重対角行列 ${\mathcal {J}}$ の方が好ましい。$\mathbf {J}$ と ${\mathcal {J}}$ は相似なので、固有値（ノード）も同じになる。重みは、行列 J から計算できる。$\phi ^{(j)}$ が固有値 x_j に対応する正規化固有ベクトル（すなわち、ユークリッドノルムが1の固有ベクトル）であるとき、固有ベクトルの第一成分から次のように重みが計算できる。

$w_{j}=\mu _{0}\left(\phi _{1}^{(j)}\right)^{2}.$

ここで $\mu _{0}$ は重み関数の積分である。

$\mu _{0}=\int _{a}^{b}w(x)dx.$

詳しくは Gil, Segura & Temme 2007 を参照されたい。

誤差見積もり

ガウス求積法の誤差は次のように定式化される。積分対象の関数が 2n 次の連続導関数を持つとき、

$\int _{a}^{b}\omega (x)\,f(x)\,dx-\sum _{i=1}^{n}w_{i}\,f(x_{i})={\frac {f^{(2n)}(\xi )}{(2n)!}}\,(p_{n},p_{n})$

となり、ξ は (a, b) にあり、p_n は n 次の直交多項式であり、さらに

$(f,g)=\int _{a}^{b}\omega (x)f(x)g(x)\,dx\,\!$

である。ω(x) = 1 となる重要な特殊ケースでは、次のような誤差見積もりがある。

${\frac {(b-a)^{2n+1}(n!)^{4}}{(2n+1)[(2n)!]^{3}}}f^{(2n)}(\xi ),\qquad a<\xi <b.\,\!$

Stoer and Bulirsch によれば、2n 次の導関数を見積もるのが難しいのでこの誤差見積もりは実用には不便であり、さらに言えば実際の誤差は導関数の界よりもずっと小さい。別の手法として、異なる次数のガウス求積法を使い、2つの結果の差分から誤差を見積もる方法もある。この場合、ガウス＝クロンロッド求積法が便利である。

ガウス＝クロンロッド求積法

区間 [a, b] を分割すると、各部分区間のガウス評価点は元の区間での評価点とは一致せず（奇数の場合の0を除く）、従って、新たに評価点を求める必要がある。ガウス＝クロンロッド求積法は、ガウス求積法の n 個の点に n + 1 個の点を追加し、求積法としての次数を 3n + 1 にするものである。これにより、低次の近似で使う関数値を高次の近似の計算に再利用できる。通常のガウス求積法とクロンロッドの拡張による近似の差分が誤差の見積もりによく利用される。

外部リンク

• ALGLIB 数値積分のためのアルゴリズムが多数掲載されている（言語は C# / C++ / Delphi / Visual Basic / その他）
• GNU Scientific Library - C言語版の QUADPACK アルゴリズムを含む（GNU Scientific Libraryも参照）
• Gaussian Quadrature Rule of Integration - Notes, PPT, Matlab, Mathematica, Maple, Mathcad at Holistic Numerical Methods Institute
• Gaussian Quadrature table at sitmo.com (リンク切れ。閲覧する場合は Internet Archive によるキャッシュ (こちら) を参照。)
• Legendre-Gauss Quadrature at MathWorld
• Gaussian Quadrature by Chris Maes and Anton Antonov, The Wolfram Demonstrations Project.