ガウス求積

出典: フリー百科事典『ウィキペディア(Wikipedia)』
移動: 案内検索

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

nを正の整数とし、f(x) を 2n-1 次以下の任意の多項式関数とすると、f(x) によらずnのみで決まる、積分点またはガウス点 (ガウスノード)と呼ばれる [-1, 1] 内のn個の点 x_i と、重み と呼ばれるn個の実数 w_i が一意的に存在して、f(x)の [-1, 1] における定積分値 I は厳密に、

 I = \int_{-1}^1 f(x)\,dx = \sum_{i=1}^n w_i f(x_i)

で与えられる。

この場合の積分点はn次のルジャンドル多項式のn個の零点に一致する。この方法をn次のルジャンドル・ガウス(Legendre-Gauss)公式と呼び、通常はガウス求積またはガウスの数値積分公式と言えばこの方法を指している[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 次多項式の根である[2][3]

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

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

 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 89
\pm\sqrt{3/5} 59
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 128225
\pm\tfrac13\sqrt{5-2\sqrt{10/7}} \tfrac{322+13\sqrt{70}}{900}
\pm\tfrac13\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

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

区間 ω(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 以下の次数の多項式について正確に積分を計算できる重み wi が存在する。さらに、それらノードは全て開区間 (a, b) にある[3]

この多項式 p_n は、重み関数 \omega (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 アルゴリズムと呼ぶ[5]

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


J\tilde{P}=x_j \tilde{P}

ここで 
\tilde{P}=[p_0 (x_j),p_1 (x_j),...,p_{n-1}(x_j)]^{T}
である。そして、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 を参照されたい[5]

誤差見積もり[編集]

ガウス求積法の誤差は次のように定式化される[3]。積分対象の関数が 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) にあり、pnn次の直交多項式であり、さらに

 (f,g) = \int_a^b \omega(x) f(x) g(x) \, dx \,\!

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

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

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

ガウス=クロンロッド求積法[編集]

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

脚注[編集]

  1. ^ 森・名取・鳥居 『数値計算』、岩波書店〈情報科学 18〉、1982年、pp130-132。
  2. ^ Press, William H.; Flannery, Brian P.; Teukolsky, Saul A.; Vetterling, William T. (1988年), “§4.5: Gaussian Quadratures and Orthogonal Polynomials”, Numerical Recipes in C (2nd ed.), Cambridge University Press, ISBN 978-0-521-43108-8 
  3. ^ a b c d Stoer, Josef; Bulirsch, Roland (2002年), Introduction to Numerical Analysis (3rd ed.), Springer, ISBN 978-0-387-95452-3 
  4. ^ a b Abramowitz, Milton; Stegun, Irene A., eds. (1972年), “§25.4, Integration”, Handbook of Mathematical Functions (with Formulas, Graphs, and Mathematical Tables), Dover, ISBN 978-0-486-61272-0 
  5. ^ a b Gil, Amparo; Segura, Javier; Temme, Nico M. (2007年), “§5.3: Gauss quadrature”, Numerical Methods for Special Functions, SIAM, ISBN 978-0-898716-34-4 
  6. ^ Kahaner, David; Moler, Cleve; Nash, Stephen (1989年), Numerical Methods and Software, Prentice-Hall, ISBN 978-0-13-627258-8 

外部リンク[編集]