ガウス求積

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

ガウス求積（ガウスきゅうせき、英: Gaussian quadrature）は、数値解析における数値積分の一種である。求積法とは、関数の定積分の近似であり、一般に積分区間内の所定の点群における関数値の重み付き総和で表される。カール・フリードリヒ・ガウスに因んで名づけられたn-点ガウス求積法とは、次数 2n - 1 以下の多項式について、点 xi と重み wi を i = 1,...,n について適切に選択することで正確な積分値を得るための求積法である。このときの積分区間は慣例的に [−1, 1] とされ、ガウス求積法は次のように式で表される。

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

ガウス求積法が正確な値を生成するのは、関数 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}（ガウス-エルミート）がある。

評価点は直交多項式のクラスに属する多項式の根である[1][2]。

目次

[編集] 基本的問題における求積

上で述べたような積分問題について、対応する多項式はルジャンドル多項式 Pn(x) である。このときの n次多項式は Pn(1) = 1 となるよう正規化され、i 番目のガウスノード xi は i 番目の Pn の根である。重みは次の式で与えられる[3]。

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\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

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

区間 ω(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) にある[2]。

この多項式 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 アルゴリズムと呼ぶ[4]。

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 を参照されたい[4]。

[編集] 誤差見積もり

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

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

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

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

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

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

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

[編集] 脚注

  1. ^ 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 
  2. ^ a b c d Stoer, Josef; Bulirsch, Roland (2002年), Introduction to Numerical Analysis (3rd ed.), Springer, ISBN 978-0-387-95452-3 
  3. ^ 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 
  4. ^ 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 
  5. ^ Kahaner, David; Moler, Cleve; Nash, Stephen (1989年), Numerical Methods and Software, Prentice-Hall, ISBN 978-0-13-627258-8 

[編集] 外部リンク

「http://ja.wikipedia.org/w/index.php?title=ガウス求積&oldid=46748743」から取得