本文主要目的是尽量通俗地介绍高斯求积公式的来龙去脉，帮助大家更好地理解高斯求积公式。
数值积分基础
先来看一般的数值积分
∫
Ω
f
(
x
)
w
(
x
)
d
x
≈
∑
i
=
1
N
w
i
f
(
x
i
)
,
\int_\Omega f(x)w(x)dx\approx\sum_{i=1}^Nw_if(x_i),
∫Ωf(x)w(x)dx≈i=1∑Nwif(xi),
这被称为N个点的机械求积公式，其中有N个节点，N个权重，总共2N个参数，如何选取这些参数，使得其代数精度最高呢？
首先来看看代数精度的范围
- 最小为0，也就是说对于常函数，一定能选到一些节点和权重使得数值积分精确成立。这很简单，节点任选，权重保证满足条件 ∑ i = 1 N w i = ∫ Ω w ( x ) \sum_{i=1}^Nw_i = \int_{\Omega}w(x) ∑i=1Nwi=∫Ωw(x)即可；
- 有上界，不可能无限大，通俗来讲是因为2N个参数不可能保证太多的方程精确成立。
因此我们假设最大代数精度为
k
m
a
x
k_{max}
kmax，通过求解以下方程组来得到相应的节点和权重，如下
∫
Ω
x
k
w
(
x
)
d
x
=
∑
i
=
1
N
w
i
x
i
k
,
k
=
0
,
1
,
…
,
k
m
a
x
\int_\Omega x^kw(x)dx = \sum_{i=1}^Nw_ix_i^k,\quad k = 0,1,\dots, k_{max}
∫Ωxkw(x)dx=i=1∑Nwixik,k=0,1,…,kmax
这个方程本身是非线性的，难以求解得到一组节点
{
x
i
}
\{x_i\}
{xi}和权重
{
w
i
}
\{w_i\}
{wi}。但仔细观察这个方程组，如果我们先有了一组确定的节点，再求解权重，问题就变得很简单了，如下
[
1
1
.
.
.
1
x
1
x
2
.
.
.
x
N
x
1
2
x
2
2
.
.
.
x
N
2
⋮
⋮
⋮
⋮
x
1
k
m
a
x
x
2
k
m
a
x
.
.
.
x
N
k
m
a
x
]
[
w
1
w
2
w
3
⋮
w
N
]
=
[
∫
Ω
w
(
x
)
d
x
∫
Ω
x
w
(
x
)
d
x
∫
Ω
x
2
w
(
x
)
d
x
⋮
∫
Ω
x
k
m
a
x
w
(
x
)
d
x
]
\begin{bmatrix} 1&1&...&1\\ x_1&x_2&...&x_N\\ x_1^2&x_2^2&...&x_N^2\\ \vdots&\vdots&\vdots&\vdots\\ x_1^{k_{max}}&x_2^{k_{max}}&...&x_N^{k_{max}} \end{bmatrix} \begin{bmatrix} w_1\\ w_2\\ w_3\\ \vdots\\ w_N \end{bmatrix} = \begin{bmatrix} \int_\Omega w(x)dx\\ \int_\Omega xw(x)dx\\ \int_\Omega x^2w(x)dx\\ \vdots\\ \int_\Omega x^{k_{max}}w(x)dx \end{bmatrix}
1x1x12⋮x1kmax1x2x22⋮x2kmax.........⋮...1xNxN2⋮xNkmax
w1w2w3⋮wN
=
∫Ωw(x)dx∫Ωxw(x)dx∫Ωx2w(x)dx⋮∫Ωxkmaxw(x)dx
将以上方程简记为
A
x
=
b
Ax=b
Ax=b，要使其有唯一解，只需要
k
m
a
x
k_{max}
kmax取为N，节点
{
x
i
}
\{x_i\}
{xi}互异，则
A
A
A就是一个Vandermond矩阵，其行列式非0，则方程组有唯一解。
到这里我们发现：
- 确定好一组互异的求积节点后，权重的求解几乎是机械的，就能使得求积公式至少有N次代数精度了；
- 虽然 A A A行列式非0，但并不代表数值上求解 A x = b Ax=b Ax=b就是容易的，求得的解的误差与 A A A的条件数有关，而求积节点的选取直接影响A的条件数。
因此求积节点的选取是至关重要的，下面我们讨论如何选取求积节点使得求积公式代数精度最大。
高斯求积公式
对被积函数
f
f
f，我们在求积节点上对其进行Lagrange插值。
f
(
x
)
=
∑
k
=
1
N
f
(
x
k
)
l
k
(
x
)
+
1
N
!
f
(
N
)
(
ξ
(
x
)
)
e
N
(
x
)
,
ξ
(
x
)
∈
Ω
e
N
(
x
)
=
(
x
−
x
1
)
(
x
−
x
2
)
…
(
x
−
x
N
)
,
f(x)=\sum_{k=1}^Nf(x_k)l_k(x)+\frac1{N!}f^{(N)}(\xi(x))e_{N}(x),\quad\xi(x)\in\Omega \\ e_{N}(x) = (x-x_1)(x-x_2)\dots(x-x_N),
f(x)=k=1∑Nf(xk)lk(x)+N!1f(N)(ξ(x))eN(x),ξ(x)∈ΩeN(x)=(x−x1)(x−x2)…(x−xN),
则
∫
Ω
f
(
x
)
w
(
x
)
d
x
=
∑
k
=
1
N
w
k
f
(
x
k
)
+
1
N
!
∫
Ω
f
(
N
)
(
ξ
(
x
)
)
e
N
(
x
)
w
(
x
)
d
x
w
k
=
∫
Ω
l
k
(
x
)
w
(
x
)
d
x
,
\int_{\Omega}f(x)w(x)\mathrm{d}x=\sum_{k=1}^Nw_kf(x_k)+\frac{1}{N!}\int_{\Omega}f^{(N)}(\xi(x))e_{N}(x)w(x)\mathrm{d}x\\ w_k=\int_{\Omega}l_k\left(x\right)w(x)\mathrm{d}x,
∫Ωf(x)w(x)dx=k=1∑Nwkf(xk)+N!1∫Ωf(N)(ξ(x))eN(x)w(x)dxwk=∫Ωlk(x)w(x)dx,
观察右边的误差项，只要
f
(
N
)
f^{(N)}
f(N)和
e
N
(
x
)
e_N(x)
eN(x)关于权函数
w
(
x
)
w(x)
w(x)正交，则误差项为0。
现在考虑 f f f为任意 k k k次多项式，则 f ( N ) f^{(N)} f(N)为 k − N k-N k−N次多项式，要让 f ( N ) f^{(N)} f(N)和 e N ( x ) e_N(x) eN(x)关于权函数 w ( x ) w(x) w(x)正交，也就是说， e N ( x ) e_N(x) eN(x)要与任意一个 k − N k-N k−N次多项式正交！,则
- 首先， k − N ≤ N − 1 k-N\le N-1 k−N≤N−1,即 k m a x ≤ 2 N − 1 k_{max}\le2N-1 kmax≤2N−1，因为若 k − N ≥ N k-N\ge N k−N≥N，则不妨取 f ( N ) = e N ( x ) f^{(N)} = e_N(x) f(N)=eN(x)，这时候误差一定不是0;
- 其次 k m a x = 2 N − 1 k_{max} = 2N-1 kmax=2N−1，只要取 e N ( x ) e_N(x) eN(x)为关于权函数 w ( x ) w(x) w(x)的 N N N次正交多项式即可。这样，由N次正交多项式与任何低于N次的多项式关于权函数 w ( x ) w(x) w(x)正交即知，误差为0。
总结一下，通过以上分析， e N ( x ) = ( x − x 1 ) ( x − x 2 ) … ( x − x N ) e_{N}(x) = (x-x_1)(x-x_2)\dots(x-x_N) eN(x)=(x−x1)(x−x2)…(x−xN)是关于权函数 w ( x ) w(x) w(x)的 N N N次正交多项式的时候，N个点的数值求积公式达到2N-1次代数精度，此时求积节点的选取也一目了然了，就取在区域 Ω \Omega Ω上带权 w ( x ) w(x) w(x)的N次正交多项式的零点。
因而，对于不同的权函数，有相应的正交多项式，使用其零点能构造出在此区域上具有最大代数精度的求积公式，这便是高斯求积公式。
一些权函数和对应的正交多项式如下表，可以根据下表构造指定区域上的高斯求积公式。
| 权函数 | 区域 | 正交多项式 |
|---|---|---|
| 1 | [ − 1 , 1 ] [-1,1] [−1,1] | Legendre |
| e − 1 2 x 2 e^{-\frac{1}{2}x^2} e−21x2 | [ − ∞ , ∞ ] [-\infty,\infty] [−∞,∞] | Probabilist Hermite |
| e − x e^{-x} e−x | [ 0 , ∞ ] [0,\infty] [0,∞] | Laguerre |
如果积分区域不是上述区域怎么办？做个线性变换即可。
例如，有界区域 [ a , b ] [a,b] [a,b]上的积分，通过积分的变量代换公式变换到区域 [ − 1 , 1 ] [-1,1] [−1,1]在进行计算即可。
2940
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
