ガウス求積のソースを表示

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

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

{{Indent|<math> I = \int_{-1}^1 f(x)\,dx = \sum_{i=1}^n w_i f(x_i)</math>}}

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

実は、{{mvar|n}} 次の[[ルジャンドル多項式]]の {{mvar|n}} 個の零点（これらは全て {{math|[−1, 1]}} 内にある）を積分点として選び、{{mvar|w{{sub|i}}}} を適切に選ぶと、{{math|''f''(''x'')}} が {{math|2''n'' − 1}} 次以下の多項式であれば上記の式が厳密に成立することが示せる。この場合、{{mvar|w{{sub|i}}}} は {{math|''f''(''x'')}} によらず一意的に定まる。この方法を {{mvar|n}} 次の'''ガウス・ルジャンドル (Gauss–Legendre) 公式'''と呼び、通常はガウス求積またはガウスの数値積分公式と言えばこの方法を指している<ref name="Mori_Natori_Torii_1982">森・名取・鳥居　『数値計算』、岩波書店〈情報科学 18〉、1982年、pp. 130–132.</ref>。

{{math|''f''(''x'')}} が {{math|2''n'' − 1}}次を超える多項式関数の場合、または多項式関数でない場合には、上記の公式は厳密には成立しないが、{{math|''f''(''x'')}} が {{math|2''n'' − 1}} 次以下の多項式関数で精度よく近似できる場合には、上記の公式を {{math|''f''(''x'')}} に対して適用することにより、その {{math|[−1, 1]}} における定積分値を精度よく得ることが期待できる。それ以外のたとえば、[[特異点]]のある関数の積分にはこの公式をそのまま適用することはできないが、被積分関数を {{math|1=''f''(''x'') = ''W''(''x'') ''g''(''x'')}} と表すことができて、{{math|''g''(''x'')}} が多項式で近似できて、{{math|''W''(''x'')}} が既知の関数（重み関数、通常は正値関数）であれば、それに対応する適切な離散的重み {{mvar|w{{sub|i}}}} を使って次のように表せる。

{{Indent|<math>\int_{-1}^1 f(x)\,dx = \int_{-1}^1 W(x) g(x)\,dx \approx \sum_{i=1}^n w_i g(x_i).</math>}}

典型的な重み関数としては、<math>W(x)=(1-x^2)^{-1/2}</math>（ガウス–チェビシェフ）や <math>W(x)=\exp(-x^2)</math>（ガウス–エルミート）がある。この場合の {{mvar|n}} 個の積分点 {{mvar|x{{sub|i}}}} はルジャンドル多項式と同様に、ある[[直交多項式]]のクラスに属する {{mvar|n}} 次多項式の根である<ref name=Press1988>{{citation | last1=Press | first1=William H. | last2=Flannery | first2=Brian P. |last3=Teukolsky | first3=Saul A. | author3-link= | last4=Vetterling | first4= William T. | chapter=&sect;4.5: Gaussian Quadratures and Orthogonal Polynomials | title=Numerical Recipes in C | edition=2nd | date=1988年 | publisher=Cambridge University Press | isbn=978-0-521-43108-8 }}</ref><ref name=Stoer2002>{{citation | last1=Stoer | first1=Josef | last2=Bulirsch | first2=Roland | date=2002年 | title=Introduction to Numerical Analysis | edition=3rd | publisher=Springer | isbn=978-0-387-95452-3<!-- 0-387-95452-X --> }}</ref>。

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

== ガウス・ルジャンドル公式による求積 == 
上述のように {{mvar|n}} 次のこの方法には、{{mvar|n}} 次のルジャンドル多項式 {{math|''P''<sub>''n''</sub>(''x'')}} が対応している。このときの {{mvar|n}} 次多項式は {{math|1=''P''<sub>''n''</sub>(1) = 1}} となるよう正規化され、{{mvar|i}} 番目のガウスノード {{mvar|x<sub>i</sub>}} は {{mvar|i}} 番目の {{mvar|P<sub>n</sub>}} の根である。重みは次の式で与えられる<ref name=Abramowitz1972>{{citation | editor1-last=Abramowitz | editor1-first=Milton | editor2-last=Stegun | editor2-first=Irene A. | chapter=&sect;25.4, Integration | title=[[Abramowitz and Stegun|Handbook of Mathematical Functions (with Formulas, Graphs, and Mathematical Tables)]] | date=1972年 | publisher=[[ドーヴァー出版|Dover]] | isbn=978-0-486-61272-0 }}</ref>。
{{Indent|<math> w_i = \frac{2}{\left( 1-x_i^2 \right) [P'_n(x_i)]^2}.</math>}}

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

{| class="wikitable" style="margin:auto; text-align:center; background:white;"
! 点の個数 {{mvar|n}} !! 点 {{mvar|x<sub>i</sub>}} !! 重み {{mvar|w<sub>i</sub>}}
|-
| {{math|1}} || {{math|0}} || {{math|2}}
|-
| {{math|2}} || <math>\pm\sqrt{1/3}</math> || {{math|1}}
|-
| rowspan="2" | {{math|3}} || {{math|0}} || {{math|8/9}}
|- 
| <math>\pm\sqrt{3/5}</math> || {{math|5/9}}
|-
| rowspan="2" | {{math|4}} || <math>\pm\sqrt{\Big( 3 - 2\sqrt{6/5} \Big)/7}</math> || <math>\tfrac{18+\sqrt{30}}{36}</math>
|-
| <math>\pm\sqrt{\Big( 3 + 2\sqrt{6/5} \Big)/7}</math> || <math>\tfrac{18-\sqrt{30}}{36}</math>
|-
| rowspan="3" | {{math|5}} || {{math|0}} || {{math|128/225}}
|-
| <math>\pm\tfrac13\sqrt{5-2\sqrt{10/7}}</math> || <math>\tfrac{322+13\sqrt{70}}{900}</math>
|-
| <math>\pm\tfrac13\sqrt{5+2\sqrt{10/7}}</math> || <math>\tfrac{322-13\sqrt{70}}{900}</math>
|}

== 区間の変更 ==
一般的な区間 {{math|[''a'', ''b'']}} についての積分は、ガウス求積法を適用する前にその区間を標準区間 {{math|[−1, 1]}} に変更する必要がある。この区間変更は以下のように線型変換で行う。

{{Indent|<math>
\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 .
</math>}}

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

{{Indent|<math>
\frac{b-a}{2} \sum_{i=1}^n w_i f\left(\frac{b-a}{2}x_i + \frac{a+b}{2}\right).
</math>}}

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

{{Indent|<math> \int_a^b \omega(x)\,f(x)\,dx .</math>}}

{{mvar|a}}, {{mvar|b}}, {{mvar|&omega;}} は適当に選択する。{{math|1=''a'' = −1}}, {{math|1=''b'' = 1}}, {{math|1=''ω''(''x'') = 1}} のとき、前述の問題と同じ形式になる。それ以外の選択では、別の求積法になる。そのうちの一部を下記の表に示す。"A & S" という欄は、[[Abramowitz and Stegun]]<ref name=Abramowitz1972 />にある式番号である。

{| class="wikitable" style="margin:auto; text-align:center; background:white;"
! 区間 !! {{math|''ω''(''x'')}} !! 直交多項式 !! A & S !! 解説など
|-
| {{math|[−1, 1]}}             || {{math|1}} || [[ルジャンドル多項式]] || 25.4.29 || 本項（上）で解説
|-
| {{math|(−1, 1)}}             || <math>(1-x)^\alpha (1+x)^\beta,\quad \alpha, \beta > -1\,</math> || [[ヤコビ多項式]] || 25.4.33 (<math>\beta=0</math>) ||
|-
| {{math|(−1, 1)}}             || <math>\frac{1}{\sqrt{1 - x^2}}</math> || [[チェビシェフ多項式]]（第一種） || 25.4.38 || {{仮リンク|チェビシェフ・ガウス求積法|en|Chebyshev–Gauss quadrature}}
|-
| {{math|[−1, 1]}}             || <math>\sqrt{1 - x^2}</math> || チェビシェフ多項式（第二種） || 25.4.40 || チェビシェフ・ガウス求積法
|-
| {{math|[0, &infin;)}}        || <math> \exp(-x) </math>|| [[ラゲール多項式]] || 25.4.45 || {{仮リンク|ガウス・ラゲール求積法|en|Gauss–Laguerre quadrature}} 
|-
| {{math|(−&infin;, &infin;)}} || <math> \exp(-x^2) </math>|| [[エルミート多項式]] || 25.4.46 || {{仮リンク|ガウス・エルミート求積法|en|Gauss–Hermite quadrature}}
|}

=== 基礎となる定理 ===
{{mvar|p{{sub|n}}}} が自明でない {{mvar|n}} 次の多項式で、次のように表されるとする。

{{Indent|<math>
\int_a^b \omega(x) \, x^k p_n(x) \, dx = 0, \quad \text{for all }k=0,1,\ldots,n-1.
</math>}}

 {{mvar|p{{sub|n}}}} のn個の[[零点]]をノード（分点）として選ぶと、次数が {{math|2''n'' − 1}} 以下の任意の多項式について正確な積分値を与えるn個の重み {{mvar|w{{sub|i}}}} を選ぶことができる。さらに、それらのノードには重複がなくすべて開区間 {{math|(''a'', ''b'')}} にある<ref name=Stoer2002 />。

この多項式 {{mvar|p{{sub|n}}}} は、{{math|''&omega;''(''x'')}} を重み関数とする次数 {{mvar|n}} の直交多項式である。

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

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

{{Indent|<math>p_{n+1}(x)+(B_n-x)p_n (x)+A_n p_{n-1}(x)=0, \qquad n=1,2,\ldots.</math>}}

このことから、対応する行列の固有値および固有ベクトルからノードと重みを計算することができる。これを一般に Golub–Welsch アルゴリズムと呼ぶ<ref name=Gil2007>{{citation | last1=Gil | first1=Amparo | last2=Segura | first2=Javier |last3=Temme | first3=Nico M. | chapter=&sect;5.3: Gauss quadrature| title=Numerical Methods for Special Functions | date=2007年 | publisher=SIAM | isbn=978-0-898716-34-4 }}</ref><ref>Walter Gautschi:"A Software Repository for Gaussian Quadratures and Christoffel Functions",SIAM,ISBN978-1611976342,(2020).</ref>。

{{mvar|x{{sub|i}}}} が直交多項式 {{mvar|p{{sub|n}}}} の根であるとき、前掲の漸化式を <math>k=0,1,\ldots, n-1</math> について用い、<math>p_n (x_j)=0</math> であることを踏まえると、次が成り立つことがわかる。

:<math>
J\tilde{P}=x_j \tilde{P}.
</math>

ここで
:<math>
\tilde{P}={}^t[p_0 (x_j),p_1 (x_j),...,p_{n-1}(x_j)]
</math>
である。そして、{{mvar|J}} はいわゆるヤコビ行列である。

:<math>
\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}
.
</math>

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

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

:<math>
w_j=\mu_0 \left(\phi_1^{(j)}\right)^2.
</math>

ここで <math>\mu_0</math> は重み関数の積分である。

:<math>
\mu_0=\int_a^b w(x) dx.
</math>

詳しくは Gil, Segura & Temme 2007 を参照されたい<ref name=Gil2007 />。

=== 誤差の見積もり ===
ガウス求積法の誤差は次のように定式化される<ref name=Stoer2002 />。被積分関数が連続な {{math|2''n''}} 次の導関数を持つときには、

{{Indent|<math> \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) </math>}}

となる。ここで {{mvar|ξ}} は区間 {{math|(''a'', ''b'')}} にあり、{{mvar|p<sub>n</sub>}} は {{mvar|n}} 次の直交多項式であり、さらに

{{Indent|<math> (f,g) = \int_a^b \omega(x) f(x) g(x) \, dx \,\!</math>}}

である。重要な特別な場合 {{math|1=''&omega;''(''x'') = 1}} については、次のような誤差見積もりがある<ref name=Kahaner1989>{{citation | last1=Kahaner | first1=David | last2=Moler | first2=Cleve | author2-link= | last3=Nash | first3=Stephen | title=Numerical Methods and Software | date=1989年 | publisher=Prentice-Hall | isbn=978-0-13-627258-8 }}</ref>。

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

Stoer and Bulirsch<ref name=Stoer2002 /> によれば、この誤差見積もりは{{math|2''n''}} 次の導関数を見積もるのが難しいので実用には不向きであり、さらに言うと実際の誤差はこの見積もりの与える上界よりもずっと小さい。別の手法として、次数の異なるガウス求積法を使って、2つの結果の違いから誤差を見積もる方法もある。それにはガウス＝クロンロッド求積法が便利である。

=== ガウス＝クロンロッド求積法 ===
{{main|ガウス＝クロンロッド求積法}}
区間 {{math|[''a'', ''b'']}} を分割すると、各部分区間のガウス評価点は元の区間での評価点とは一致せず（奇数の場合の0を除く）、従って、新たに評価点を求める必要がある。[[ガウス＝クロンロッド求積法]]<ref>Notaris, S. E. (2016). Gauss–Kronrod quadrature formulae–a survey of fifty years of research. Electron. Trans. Numer. Anal, 45, 371-404.</ref><ref>Gauss-Kronrod quadrature formula. Encyclopedia of Mathematics. URL: http://www.encyclopediaofmath.org/index.php?title=Gauss-Kronrod_quadrature_formula&oldid=22491</ref>は、ガウス求積法の {{mvar|n}} 個の点に {{math|''n'' + 1}} 個の点を追加し、求積法としての次数を {{math|2''n'' + 1}} にするものである。これにより、低次の近似で使う関数値を高次の近似の計算に再利用できる。通常のガウス求積法とクロンロッドの拡張による近似の差分が誤差の見積もりによく利用される。

== 脚注 ==
{{Reflist}}

== 外部リンク ==
* [http://www.alglib.net/integral/gq/ ALGLIB] 数値積分のためのアルゴリズムが多数掲載されている（言語は C# / C++ / Delphi / Visual Basic / その他）
* [https://www.gnu.org/software/gsl/ GNU Scientific Library] - [[C言語]]版の QUADPACK アルゴリズムを含む（[[GNU Scientific Library]]も参照）
* [http://nm.mathforcollege.com/topics/gauss_quadrature.html Gaussian Quadrature Rule of Integration - Notes, PPT, Matlab, Mathematica, Maple, Mathcad] at ''Holistic Numerical Methods Institute''
* [http://www.sitmo.com/eqcat/13 Gaussian Quadrature table] at sitmo.com (リンク切れ。閲覧する場合は [[インターネットアーカイブ|Internet Archive]] によるキャッシュ ([https://web.archive.org/web/20100413154934/http://www.sitmo.com/eqcat/13 こちら]) を参照。)
* [http://mathworld.wolfram.com/Legendre-GaussQuadrature.html Legendre-Gauss Quadrature at MathWorld]
* [http://demonstrations.wolfram.com/GaussianQuadrature/ Gaussian Quadrature] by Chris Maes and Anton Antonov, [[ウルフラム・リサーチ|The Wolfram Demonstrations Project]].

{{integral}}

{{DEFAULTSORT:かうすきゆうせき}}
[[Category:数値積分]]
[[Category:カール・フリードリヒ・ガウス]]
[[Category:数学に関する記事]]