「ガウス求積」の版間の差分

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

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

この場合の積分点実は、{{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}} 次の'''ガウス・ルジャンドル・ガウス (Legendre-GaussGauss–Legendre) 公式'''と呼び、通常はガウス求積またはガウスの数値積分公式と言えばこの方法を指している<ref name="Mori_Natori_Torii_1982">森・名取・鳥居　『数値計算』、岩波書店〈情報科学 18〉、1982年、pp130-132。pp. 130–132.</ref>。

<{{math>|''f''(''x'')</math>}} が {{math|2''n'' − 1}}次を超える多項式関数の場合、または多項式関数でない場合におは、上記の公式は厳密には成立しないてもが、2n-{{math|''f''(''x'')}} が {{math|2''n'' − 1}} 次以下の多項式関数で精度よく近似できる場合には、上記の公式を <{{math>|''f''(''x'')</math>}} に対して適用することにより、その {{math|[-1−1, 1]}} における定積分値を精度よく得ることが期待できる。それ以外の例たとえば、[[特異点]]のある関数の積分にはこの公式をそのまま適用することはできないが、対象の被積分関数を <{{math>|1=''f''(''x'') = ''W''(''x'') ''g''(''x'')</math>}} と表すことができて、{{math|''g''(''x'')}} が近似多項式で 近似できて、{{math|''W''(''x'')}} が既知の関数（重み関数、通常は正値関数）であれば、それを代替に対応する適切な離散的重み <math>w_i</math>{{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)=e^{\exp(-x^2})</math>（ガウス-–エルミート）がある。この場合の {{mvar|n}} 個の積分点 <math>x_i</math>{{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>。

! 点の個数 ''{{mvar|n''}} !! 点 ''{{mvar|x''_''i''}} !! 重み ''{{mvar|w''_''i''}}

| {{math|1}} || {{math|0}} || {{math|2}}

| {{math|2}} || <math>\pm\sqrt{1/3}</math> || {{math|1}}

| rowspan="2" | {{math|3}} || {{math|0}} || ^{{math|8&frasl;₉}}

| <math>\pm\sqrt{3/5}</math> || ^{{math|5&frasl;₉}}

| rowspan="2" | {{math|4}} || <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&frasl;₂₂₅}}

== 区間の変更 ==

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

+ \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).

正の重み関数 {{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 />にある式番号である。

! 区間 !! {{math|''ω''(''x'')}} !! 直交多項式 !! A & S !! 解説など

| {{math|[−1, 1]}} || <{{math>|1\,</math>}} || [[ルジャンドル多項式]] || 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{{仮リンク|チェビシェフ・ガウス求積法]]|en|Chebyshev–Gauss quadrature}}

| {{math|[−1, 1]}} || <math>\sqrt{1 - x^2}</math> || チェビシェフ多項式（第二種） || 25.4.40 || [[:en:Chebyshev–Gauss quadrature|チェビシェフ・ガウス求積法]]

| {{math|[0, &infin;)}} || <math> e^{\exp(-x}\,) </math> || [[ラゲール多項式]] || 25.4.45 || [[:en:Gauss–Laguerre quadrature{{仮リンク|ガウス・ラゲール求積法]]|en|Gauss–Laguerre quadrature}}

| {{math|(−&infin;, &infin;)}} || <math> e^{\exp(-x^2}) </math> || [[エルミート多項式]] || 25.4.46 || [[:en:Gauss–Hermite quadrature{{仮リンク|ガウス・エルミート求積法]]|en|Gauss–Hermite quadrature}}

=== 基礎的となる定理 ===

<math>p_n</math>{{mvar|p{{sub|n}}}} が自明でない ''{{mvar|n''}} 次の多項式で、次のように表されるとする。

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

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

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

例えば、<math>p_n</math>{{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>。

<math>x_i</math>{{mvar|x{{sub|i}}}} が直交多項式 <math>p_n</math>{{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>

\tilde{P}={}^t[p_0 (x_j),p_1 (x_j),...,p_{n-1}(x_j)]^{T}

である。そして、<math>{{mvar|J</math>}} はいわゆるヤコビ行列である。

:<math>

\mathbfboldsymbol{J}=\left(

\begin{array}{llllllpmatrix}

\end{arraypmatrix}

重みとノードを求めるには、要素が <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> は[[行列の相似|相似]]なので、固有値（ノード）も同じになる。重みは、行列 <math>{{mvar|J</math>}} から計算できる。<math>\phi^{(j)}</math> が固有値 <math>x_j</math>{{mvar|x{{sub|j}}}} に対応する正規化固有ベクトル（すなわち、ユークリッドノルムが1の固有ベクトル）であるとき、固有ベクトルの第一成分から次のように重みが計算できる。

:<math>

w_j=\mu_0 \left(\phi_1^{(j)}\right)^2.

:<math>

\mu_0=\int_a^b w(x) dx.

=== 誤差の見積もり ===

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

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

である。重要な特別な場合 {{math|1=''&omega;''(''x'')&nbsp; = 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>。

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

区間 {{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>は、ガウス求積法の <math>{{mvar|n</math>}} 個の点に <{{math>|''n'' + 1</math>}} 個の点を追加し、求積法としての次数を <{{math>3n|2''n'' + 1</math>}} にするものである。これにより、低次の近似で使う関数値を高次の近似の計算に再利用できる。通常のガウス求積法とクロンロッドの拡張による近似の差分が誤差の見積もりによく利用される。

* [http://www.sitmo.com/eqcat/13 Gaussian Quadrature table] at sitmo.com (リンク切れ。閲覧する場合は [[インターネット・アーカイブ|Internet Archive]] によるキャッシュ ([httphttps://web.archive.org/web/20100413154934/http://www.sitmo.com/eqcat/13 こちら]) を参照。)