ガウス求積 （ガウスきゅうせき、英 : Gaussian quadrature ）またはガウスの数値積分公式 とは、カール・フリードリヒ・ガウス に因んで名づけられた数値解析 における数値積分 法の一種であり、実数 のある閉区間(慣例的に [−1, 1] に標準化される)で定義された実数値関数のその閉区間に渡る定積分値を、比較的少ない演算で精度良く求めることができるアルゴリズム である。
n を正の整数 とし、f (x )多項式関数 とする。f (x )[−1, 1] に渡る定積分値 I を、
I
=
∫ ∫ -->
− − -->
1
1
f
(
x
)
d
x
=
∑ ∑ -->
i
=
1
n
w
i
f
(
x
i
)
{\displaystyle I=\int _{-1}^{1}f(x)\,dx=\sum _{i=1}^{n}w_{i}f(x_{i})}
の形でなるべく正確に近似する公式を考える。ここで、xi は積分点 またはガウス点 (ガウスノード) と呼ばれる [−1, 1] 内の n 個の点であり、wi は重み と呼ばれるn 個の実数である。
実は、n 次のルジャンドル多項式 の n 個の零点（これらは全て [−1, 1] 内にある）を積分点として選び、wi を適切に選ぶと、f (x )2n − 1 次以下の多項式であれば上記の式が厳密に成立することが示せる。この場合、wi は f (x )n 次のガウス・ルジャンドル (Gauss–Legendre) 公式 と呼び、通常はガウス求積またはガウスの数値積分公式と言えばこの方法を指している[ 1]
f (x )2n − 1 次を超える多項式関数の場合、または多項式関数でない場合には、上記の公式は厳密には成立しないが、f (x )2n − 1 次以下の多項式関数で精度よく近似できる場合には、上記の公式を f (x )[−1, 1] における定積分値を精度よく得ることが期待できる。それ以外のたとえば、特異点 のある関数の積分にはこの公式をそのまま適用することはできないが、被積分関数を f (x ) = W (x ) g (x )g (x )W (x )wi を使って次のように表せる。
∫ ∫ -->
− − -->
1
1
f
(
x
)
d
x
=
∫ ∫ -->
− − -->
1
1
W
(
x
)
g
(
x
)
d
x
≈ ≈ -->
∑ ∑ -->
i
=
1
n
w
i
g
(
x
i
)
.
{\displaystyle \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
{\displaystyle W(x)=(1-x^{2})^{-1/2}}
W
(
x
)
=
exp
-->
(
− − -->
x
2
)
{\displaystyle W(x)=\exp(-x^{2})}
n 個の積分点 xi はルジャンドル多項式と同様に、ある直交多項式 のクラスに属する n 次多項式の根である[ 2] [ 3]
重み関数と指定区間に付随するn次の直交多項式を考え、それの区間内にあるn 個の零点を分点にとして被積分関数f (x )をHermite補間公式で近似したものを考えると、直交多項式の重み関数に対する直交性から、f (x )に重み関数を掛けて積分したものは、直交関数のn 個の零点に於けるf (x )の関数値それぞれに重みをかけたものの和で近似される（結果的にf (x )の各分点における導関数値は積分の近似値には寄与しない）。このようにして重み関数に対応するガウス型の数値積分公式を導くことができて、分点がn であるときには被積分関数が2n −1次以下の任意の多項式に対して正確な積分値を与えるということが示せる.
上述のように n 次のこの方法には、n 次のルジャンドル多項式 P n x )n 次多項式は P n i 番目のガウスノード xi は i 番目の Pn の根である。重みは次の式で与えられる[ 4]
w
i
=
2
(
1
− − -->
x
i
2
)
[
P
n
′
(
x
i
)
]
2
.
{\displaystyle w_{i}={\frac {2}{\left(1-x_{i}^{2}\right)[P'_{n}(x_{i})]^{2}}}.}
低次の求積法は次のようになる。
点の個数 n
点 xi
重み wi
1 0 2
2
± ± -->
1
/
3
{\displaystyle \pm {\sqrt {1/3}}}
1
3 0 8/9
± ± -->
3
/
5
{\displaystyle \pm {\sqrt {3/5}}}
5/9
4
± ± -->
(
3
− − -->
2
6
/
5
)
/
7
{\displaystyle \pm {\sqrt {{\Big (}3-2{\sqrt {6/5}}{\Big )}/7}}}
18
+
30
36
{\displaystyle {\tfrac {18+{\sqrt {30}}}{36}}}
± ± -->
(
3
+
2
6
/
5
)
/
7
{\displaystyle \pm {\sqrt {{\Big (}3+2{\sqrt {6/5}}{\Big )}/7}}}
18
− − -->
30
36
{\displaystyle {\tfrac {18-{\sqrt {30}}}{36}}}
5 0 128/225
± ± -->
1
3
5
− − -->
2
10
/
7
{\displaystyle \pm {\tfrac {1}{3}}{\sqrt {5-2{\sqrt {10/7}}}}}
322
+
13
70
900
{\displaystyle {\tfrac {322+13{\sqrt {70}}}{900}}}
± ± -->
1
3
5
+
2
10
/
7
{\displaystyle \pm {\tfrac {1}{3}}{\sqrt {5+2{\sqrt {10/7}}}}}
322
− − -->
13
70
900
{\displaystyle {\tfrac {322-13{\sqrt {70}}}{900}}}
一般的な区間 [a , b ] についての積分は、ガウス求積法を適用する前にその区間を標準区間 [−1, 1] に変更する必要がある。この区間変更は以下のように線型変換で行う。
∫ ∫ -->
a
b
f
(
x
)
d
x
=
b
− − -->
a
2
∫ ∫ -->
− − -->
1
1
f
(
b
− − -->
a
2
x
+
a
+
b
2
)
d
x
.
{\displaystyle \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.}
ガウス求積法を適用すると、以下で積分の近似値が得られる。
b
− − -->
a
2
∑ ∑ -->
i
=
1
n
w
i
f
(
b
− − -->
a
2
x
i
+
a
+
b
2
)
.
{\displaystyle {\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] 以外にも適用可能である。すなわち、次の形式の問題である。
∫ ∫ -->
a
b
ω ω -->
(
x
)
f
(
x
)
d
x
.
{\displaystyle \int _{a}^{b}\omega (x)\,f(x)\,dx.}
a , b , ω は適当に選択する。a = −1b = 1ω (x ) = 1Abramowitz and Stegun [ 4]
区間
ω (x )直交多項式
A & S
解説など
[−1, 1] 1 ルジャンドル多項式 25.4.29
本項（上）で解説
(−1, 1)
(
1
− − -->
x
)
α α -->
(
1
+
x
)
β β -->
,
α α -->
,
β β -->
>
− − -->
1
{\displaystyle (1-x)^{\alpha }(1+x)^{\beta },\quad \alpha ,\beta >-1\,}
ヤコビ多項式 25.4.33 (
β β -->
=
0
{\displaystyle \beta =0}
(−1, 1)
1
1
− − -->
x
2
{\displaystyle {\frac {1}{\sqrt {1-x^{2}}}}}
チェビシェフ多項式 （第一種）25.4.38
チェビシェフ・ガウス求積法 （英語版 ）
[−1, 1]
1
− − -->
x
2
{\displaystyle {\sqrt {1-x^{2}}}}
チェビシェフ多項式（第二種）
25.4.40
チェビシェフ・ガウス求積法
[0, ∞)
exp
-->
(
− − -->
x
)
{\displaystyle \exp(-x)}
ラゲール多項式 25.4.45
ガウス・ラゲール求積法 （英語版 ）
(−∞, ∞)
exp
-->
(
− − -->
x
2
)
{\displaystyle \exp(-x^{2})}
エルミート多項式 25.4.46
ガウス・エルミート求積法 （英語版 ）
pn が自明でない n 次の多項式で、次のように表されるとする。
∫ ∫ -->
a
b
ω ω -->
(
x
)
x
k
p
n
(
x
)
d
x
=
0
,
for all
k
=
0
,
1
,
… … -->
,
n
− − -->
1.
{\displaystyle \int _{a}^{b}\omega (x)\,x^{k}p_{n}(x)\,dx=0,\quad {\text{for all }}k=0,1,\ldots ,n-1.}
pn のn個の零点 をノード（分点）として選ぶと、次数が 2n − 1 以下の任意の多項式について正確な積分値を与えるn個の重み wi を選ぶことができる。さらに、それらのノードには重複がなくすべて開区間 (a , b ) にある[ 3] この多項式 pn は、ω (x )n の直交多項式である。
ガウス求積法のノード xi と重み wi を計算するための基本的ツールは、直交多項式群と対応する重み関数が満たす3項漸化式である。
例えば、pn がモニックな n 次直交多項式（最高次の項の係数が 1 の n 次直交多項式）なら、次のような漸化式 で関係を表すことができる。
p
n
+
1
(
x
)
+
(
B
n
− − -->
x
)
p
n
(
x
)
+
A
n
p
n
− − -->
1
(
x
)
=
0
,
n
=
1
,
2
,
… … -->
.
{\displaystyle 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] [ 6]
xi が直交多項式 pn の根であるとき、前掲の漸化式を
k
=
0
,
1
,
… … -->
,
n
− − -->
1
{\displaystyle k=0,1,\ldots ,n-1}
p
n
(
x
j
)
=
0
{\displaystyle p_{n}(x_{j})=0}
J
P
~ ~ -->
=
x
j
P
~ ~ -->
.
{\displaystyle J{\tilde {P}}=x_{j}{\tilde {P}}.}
ここで
P
~ ~ -->
=
t
[
p
0
(
x
j
)
,
p
1
(
x
j
)
,
.
.
.
,
p
n
− − -->
1
(
x
j
)
]
{\displaystyle {\tilde {P}}={}^{t}[p_{0}(x_{j}),p_{1}(x_{j}),...,p_{n-1}(x_{j})]}
である。そして、J はいわゆるヤコビ行列である。
J
=
(
B
0
1
0
… … -->
… … -->
… … -->
A
1
B
1
1
0
… … -->
… … -->
0
A
2
B
2
1
0
… … -->
… … -->
… … -->
… … -->
… … -->
… … -->
… … -->
… … -->
… … -->
… … -->
A
n
− − -->
2
B
n
− − -->
2
1
… … -->
… … -->
… … -->
… … -->
A
n
− − -->
1
B
n
− − -->
1
)
.
{\displaystyle {\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}}.}
したがって、ガウス求積法のノードは三重対角行列 の固有値として計算できる。
重みとノードを求めるには、要素が
J
i
,
i
=
J
i
,
i
{\displaystyle {\mathcal {J}}_{i,i}=J_{i,i}}
i
=
1
,
… … -->
,
n
{\displaystyle i=1,\ldots ,n}
J
i
− − -->
1
,
i
=
J
i
,
i
− − -->
1
=
J
i
,
i
− − -->
1
J
i
− − -->
1
,
i
,
i
=
2
,
… … -->
,
n
{\displaystyle {\mathcal {J}}_{i-1,i}={\mathcal {J}}_{i,i-1}={\sqrt {J_{i,i-1}J_{i-1,i}}},\,i=2,\ldots ,n}
対称 な三重対角行列
J
{\displaystyle {\mathcal {J}}}
J
{\displaystyle \mathbf {J} }
J
{\displaystyle {\mathcal {J}}}
相似 なので、固有値（ノード）も同じになる。重みは、行列 J から計算できる。
ϕ ϕ -->
(
j
)
{\displaystyle \phi ^{(j)}}
xj に対応する正規化固有ベクトル（すなわち、ユークリッドノルムが1の固有ベクトル）であるとき、固有ベクトルの第一成分から次のように重みが計算できる。
w
j
=
μ μ -->
0
(
ϕ ϕ -->
1
(
j
)
)
2
.
{\displaystyle w_{j}=\mu _{0}\left(\phi _{1}^{(j)}\right)^{2}.}
ここで
μ μ -->
0
{\displaystyle \mu _{0}}
μ μ -->
0
=
∫ ∫ -->
a
b
w
(
x
)
d
x
.
{\displaystyle \mu _{0}=\int _{a}^{b}w(x)dx.}
詳しくは Gil, Segura & Temme 2007 を参照されたい[ 5]
ガウス求積法の誤差は次のように定式化される[ 3] 2n 次の導関数を持つときには、
∫ ∫ -->
a
b
ω ω -->
(
x
)
f
(
x
)
d
x
− − -->
∑ ∑ -->
i
=
1
n
w
i
f
(
x
i
)
=
f
(
2
n
)
(
ξ ξ -->
)
(
2
n
)
!
(
p
n
,
p
n
)
{\displaystyle \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
)
=
∫ ∫ -->
a
b
ω ω -->
(
x
)
f
(
x
)
g
(
x
)
d
x
{\displaystyle (f,g)=\int _{a}^{b}\omega (x)f(x)g(x)\,dx\,\!}
である。重要な特別な場合 ω (x ) = 1[ 7]
(
b
− − -->
a
)
2
n
+
1
(
n
!
)
4
(
2
n
+
1
)
[
(
2
n
)
!
]
3
f
(
2
n
)
(
ξ ξ -->
)
,
a
<
ξ ξ -->
<
b
.
{\displaystyle {\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を除く）、従って、新たに評価点を求める必要がある。ガウス＝クロンロッド求積法 [ 8] [ 9] n 個の点に n + 12n + 1 にするものである。これにより、低次の近似で使う関数値を高次の近似の計算に再利用できる。通常のガウス求積法とクロンロッドの拡張による近似の差分が誤差の見積もりによく利用される。
^ 森・名取・鳥居 『数値計算』、岩波書店〈情報科学 18〉、1982年、pp. 130–132.
^ 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 ^ a b c d Stoer, Josef; Bulirsch, Roland (2002年), Introduction to Numerical Analysis (3rd ed.), Springer, ISBN 978-0-387-95452-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
^ 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
^ Walter Gautschi:"A Software Repository for Gaussian Quadratures and Christoffel Functions",SIAM,ISBN978-1611976342,(2020).
^ Kahaner, David; Moler, Cleve; Nash, Stephen (1989年), Numerical Methods and Software , Prentice-Hall, ISBN 978-0-13-627258-8 ^ Notaris, S. E. (2016). Gauss–Kronrod quadrature formulae–a survey of fifty years of research. Electron. Trans. Numer. Anal, 45, 371-404.
^ Gauss-Kronrod quadrature formula. Encyclopedia of Mathematics. URL: http://www.encyclopediaofmath.org/index.php?title=Gauss-Kronrod_quadrature_formula&oldid=22491
![{\displaystyle w_{i}={\frac {2}{\left(1-x_{i}^{2}\right)[P'_{n}(x_{i})]^{2}}}.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ced974d823852285fd484f2fb39ae82744f0637f)
![{\displaystyle {\tilde {P}}={}^{t}[p_{0}(x_{j}),p_{1}(x_{j}),...,p_{n-1}(x_{j})]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/2f7c26ac75f9535281a7d4b859f8825ab95234d2)
![{\displaystyle {\frac {(b-a)^{2n+1}(n!)^{4}}{(2n+1)[(2n)!]^{3}}}f^{(2n)}(\xi ),\qquad a<\xi <b.\,\!}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b45b889d7e11000ada8c0dfb98431306b267731a)
