ガウス求積

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

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


${\displaystyle I=\int _{-1}^{1}f($x)\,dx=\sum _{i=1}^{n}w_{i}f(x_{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を使って次のように表せる。


${\displaystyle \int _{-1}^{1}f($x)\,dx=\int _{-1}^{1}W(x)g(x)\,dx\approx \sum _{i=1}^{n}w_{i}g(x_{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}).$

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

重み関数と指定区間に付随するn次の直交多項式を考え、それの区間内にあるn個の零点を分点にとして被積分関数f(x)をHermite補間公式で近似したものを考えると、直交多項式の重み関数に対する直交性から、f(x)に重み関数を掛けて積分したものは、直交関数のn個の零点に於けるf(x)の関数値それぞれに重みをかけたものの和で近似される︵結果的にf(x)の各分点における導関数値は積分の近似値には寄与しない︶。このようにして重み関数に対応するガウス型の数値積分公式を導くことができて、分点がnであるときには被積分関数が2n−1次以下の任意の多項式に対して正確な積分値を与えるということが示せる.

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

上述のように 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] に変更する必要がある。この区間変更は以下のように線型変換で行う。


${\displaystyle \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] 以外にも適用可能である。すなわち、次の形式の問題である。


${\displaystyle \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, ∞)	$\exp(-x)$	ラゲール多項式	25.4.45	ガウス・ラゲール求積法
(−∞, ∞)	$\exp(-x^{2})$	エルミート多項式	25.4.46	ガウス・エルミート求積法

基礎となる定理

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


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

pn のn個の零点をノード（分点）として選ぶと、次数が 2n − 1 以下の任意の多項式について正確な積分値を与えるn個の重み wi を選ぶことができる。さらに、それらのノードには重複がなくすべて開区間 (a, b) にある。

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

計算

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

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


${\displaystyle 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 はいわゆるヤコビ行列である。

${\boldsymbol {J}}={\begin{pmatrix}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{pmatrix}}.$ 

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

重みとノードを求めるには、要素が ${\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から計算できる。${\displaystyle \phi ^{($j)}}   が固有値 x_jに対応する正規化固有ベクトル︵すなわち、ユークリッドノルムが1の固有ベクトル︶であるとき、固有ベクトルの第一成分から次のように重みが計算できる。

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

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

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

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

誤差の見積もり

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


${\displaystyle \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次の直交多項式であり、さらに


${\displaystyle (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 個の点を追加し、求積法としての次数を 2n+ 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.