ガウス求積

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

テンプレート:Mvar を正の整数とし、テンプレート:Math を 任意の多項式関数とする。テンプレート:Math の テンプレート:Math に渡る定積分値 テンプレート:Mvar を、

の形でなるべく正確に近似する公式を考える。ここで、テンプレート:Mvar は積分点またはガウス点 (ガウスノード)と呼ばれる テンプレート:Math 内の テンプレート:Mvar 個の点であり、テンプレート:Mvar は重みと呼ばれるn個の実数である。

実は、テンプレート:Mvar 次のルジャンドル多項式の テンプレート:Mvar 個の零点（これらは全て テンプレート:Math 内にある）を積分点として選び、テンプレート:Mvar を適切に選ぶと、テンプレート:Math が テンプレート:Math 次以下の多項式であれば上記の式が厳密に成立することが示せる。この場合、テンプレート:Mvar は テンプレート:Math によらず一意的に定まる。この方法を テンプレート:Mvar 次のガウス・ルジャンドル (Gauss–Legendre) 公式と呼び、通常はガウス求積またはガウスの数値積分公式と言えばこの方法を指している。

テンプレート:Math が テンプレート:Math次を超える多項式関数の場合、または多項式関数でない場合には、上記の公式は厳密には成立しないが、テンプレート:Math が テンプレート:Math 次以下の多項式関数で精度よく近似できる場合には、上記の公式を テンプレート:Math に対して適用することにより、その テンプレート:Math における定積分値を精度よく得ることが期待できる。それ以外のたとえば、特異点のある関数の積分にはこの公式をそのまま適用することはできないが、被積分関数を テンプレート:Math と表すことができて、テンプレート:Math が多項式で近似できて、テンプレート:Math が既知の関数（重み関数、通常は正値関数）であれば、それに対応する適切な離散的重み テンプレート:Mvar を使って次のように表せる。

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

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

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

上述のように テンプレート:Mvar 次のこの方法には、テンプレート:Mvar 次のルジャンドル多項式 テンプレート:Math が対応している。このときの テンプレート:Mvar 次多項式は テンプレート:Math となるよう正規化され、テンプレート:Mvar 番目のガウスノード テンプレート:Mvar は テンプレート:Mvar 番目の テンプレート:Mvar の根である。重みは次の式で与えられる。 テンプレート:Indent

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

点の個数 テンプレート:Mvar	点 テンプレート:Mvar	重み テンプレート:Mvar
テンプレート:Math	テンプレート:Math	テンプレート:Math
テンプレート:Math	$\pm \sqrt{1/3}$	テンプレート:Math
テンプレート:Math	テンプレート:Math	テンプレート:Math
	$\pm \sqrt{3/5}$	テンプレート:Math
テンプレート:Math	$\pm \sqrt{(3-2\sqrt{6/5})/7}$	$\frac{18+\sqrt{30}}{36}$
	$\pm \sqrt{(3+2\sqrt{6/5})/7}$	$\frac{18-\sqrt{30}}{36}$
テンプレート:Math	テンプレート:Math	テンプレート:Math
	$\pm \frac{1}{3}\sqrt{5-2\sqrt{10/7}}$	$\frac{322+13\sqrt{70}}{900}$
	$\pm \frac{1}{3}\sqrt{5+2\sqrt{10/7}}$	$\frac{322-13\sqrt{70}}{900}$

区間の変更

一般的な区間 テンプレート:Math についての積分は、ガウス求積法を適用する前にその区間を標準区間 テンプレート:Math に変更する必要がある。この区間変更は以下のように線型変換で行う。

ガウス求積法を適用すると、以下で積分の近似値が得られる。

他の形式

正の重み関数 テンプレート:Mvar を導入することで、より汎用的な積分問題の表現も可能であり、区間 テンプレート:Math 以外にも適用可能である。すなわち、次の形式の問題である。

テンプレート:Mvar, テンプレート:Mvar, テンプレート:Mvar は適当に選択する。テンプレート:Math, テンプレート:Math, テンプレート:Math のとき、前述の問題と同じ形式になる。それ以外の選択では、別の求積法になる。そのうちの一部を下記の表に示す。"A & S" という欄は、Abramowitz and Stegunにある式番号である。

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

基礎となる定理

テンプレート:Mvar が自明でない テンプレート:Mvar 次の多項式で、次のように表されるとする。

テンプレート:Mvar のn個の零点をノード（分点）として選ぶと、次数が テンプレート:Math 以下の任意の多項式について正確な積分値を与えるn個の重み テンプレート:Mvar を選ぶことができる。さらに、それらのノードには重複がなくすべて開区間 テンプレート:Math にある。

この多項式 テンプレート:Mvar は、テンプレート:Math を重み関数とする次数 テンプレート:Mvar の直交多項式である。

計算

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

例えば、テンプレート:Mvar がモニックな テンプレート:Mvar 次直交多項式（最高次の項の係数が テンプレート:Math の テンプレート:Mvar 次直交多項式）なら、次のような漸化式で関係を表すことができる。

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

テンプレート:Mvar が直交多項式 テンプレート:Mvar の根であるとき、前掲の漸化式を $k=0,1,\dots ,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)]$

である。そして、テンプレート:Mvar はいわゆるヤコビ行列である。

$𝑱=\left(\begin{array}{cccccc}{B}_0& 1& 0& \dots & \dots & \dots \\ A_1& B_1& 1& 0& \dots & \dots \\ 0& A_2& B_2& 1& 0& \dots \\ \dots & \dots & \dots & \dots & \dots & \dots \\ \dots & \dots & \dots & A_{n-2}& B_{n-2}& 1\\ \dots & \dots & \dots & \dots & A_{n-1}& B_{n-1}\end{array}\right).$

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

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

$w_j={\mu }_0{\left({\varphi }_1^{(j)}\right)}^2.$

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

${\mu }_0=\underset{a}{\overset{b}{\int }}w(x)dx.$

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

誤差の見積もり

ガウス求積法の誤差は次のように定式化される。被積分関数が連続な テンプレート:Math 次の導関数を持つときには、

となる。ここで テンプレート:Mvar は区間 テンプレート:Math にあり、テンプレート:Mvar は テンプレート:Mvar 次の直交多項式であり、さらに

である。重要な特別な場合 テンプレート:Math については、次のような誤差見積もりがある。

Stoer and Bulirsch によれば、この誤差見積もりはテンプレート:Math 次の導関数を見積もるのが難しいので実用には不向きであり、さらに言うと実際の誤差はこの見積もりの与える上界よりもずっと小さい。別の手法として、次数の異なるガウス求積法を使って、2つの結果の違いから誤差を見積もる方法もある。それにはガウス＝クロンロッド求積法が便利である。

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

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

外部リンク

• 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.