Gauss 求积公式

在数值分析中，Gauss 求积公式是一种数值积分方法。

基本原理

设有积分$\int_a^b \rho(x) f(x) \mathrm{d}x$，其中，$\rho(x)$是已知（可积）函数，$f(x)$是可积函数空间上的一个函数，在$[a, b]$上取节点$a \leqslant x_0 < x_1 < \cdots < x_n \leqslant b$，考虑求积公式

$$\int _{a}^{b}\rho (x)f(x)\mathrm {d} x\approx \sum _{k=0}^{n}A_{k}f(x_{k}).$$

这里有$2n+2$个未知量需要确定。

令$p_k (x) = x^k, k = 0,1,\cdots,2n+1$代入可求解，得到的公式具有$2n+1$次代数精度。这样的节点称为 Gauss 点，公式称为 Gauss 型求积公式。

充要条件

$x_k, k = 0,1,\cdots,n$为 Gauss 点的充要条件是多项式$\omega(x) = \prod_{k=0}^n (x - x_k)$与任意次数不大于$n$的多项式$p(x)$带权正交，即

$$\forall p(x) \in \R_n [x], \int_a^b \rho(x) \omega(x) p(x) \mathrm{d}x = 0.$$

注意到正交多项式族$\{ \varphi_k (x) \}_{k=0}^\infty$有性质：任意次数不大于$n$的多项式$p(x)$必与$\varphi_{n+1} (x)$正交。因此若取$\omega(x) = \varphi_{n+1}(x)$，那么$\varphi_{n+1} (x)$的根就是$n+1$个 Gauss 点。

误差估计

假设同上，设$P(x)$是节点$\{ x_k \}_{k=0}^n$的 Hermite 插值多项式，则有$f(x_k) = P(x_k), k = 0,1,2,\cdots,n$，那么误差

$$\begin{align} R[f] & = \int_a^b \rho(x) f(x) \mathrm{d}x - \sum_{k=0}^n A_k f(x_k) \\ & = \int_a^b \rho(x) f(x) \mathrm{d}x - \sum_{k=0}^n A_k P(x_k) \\ & = \int_a^b \rho(x) f(x) \mathrm{d}x - \int_a^b \rho(x) P(x) \mathrm{d}x \\ & = \int_a^b \rho(x) [f(x) - P(x)] \mathrm{d}x \\ & = \dfrac{f^{(2n+2)} (\xi)}{(2n+2)!} \int_a^b \rho(x) \omega^2 (x) \mathrm{d}x. \end{align}$$

第三个等号是因为求积公式#A1是插值型的，代数精度是$2n+1$，因此对$2n+1$的多项式$P(x)$精确成立，最后一个等号是 Lagrange 插值的误差。

还可以证明，求积公式#A1的系数均为正数，因此求积公式在有限区间$[a, b]$上是数值稳定的。

例子

Gauss 求积公式一般先要构造正交的函数序列，这依赖于权函数$\rho(x)$的选择，例如，当$[a, b] = [-1, 1], \rho(x) = 1$时正交函数对应于 Legendre 多项式；当$[a, b] = [-1, 1], \rho(x) = \dfrac{1}{\sqrt{1-x^2}}$时正交函数对应于第一类 Chebyshev 多项式，此时$x_k = \cos\dfrac{(2k-1)\pi}{2n+2}, k = 1,2,\cdots,n+1$，此时

$$\int_{-1}^1 \dfrac{f(x) \mathrm{d}x}{\sqrt{1-x^2}} \approx \dfrac{\pi}{n+1} \sum_{k=0}^n f\!\left( \cos\dfrac{(2k-1)\pi}{2n+2} \right).$$

对于一般情形，则要设法构造多项式。假设有带权内积$(\varphi, \psi) = \int_a^b \rho(x) \varphi(x) \psi(x) \mathrm{d}x$，假设$\int_a^b \rho(x) f(x) \approx \sum_{k=0}^n A_k f(x_k).$我们需要的首先是确定 Gauss 点$x_0, x_1, \cdots, x_n$，为此可以使用 Gram-Schmidt 正交化对$p_k (x) = x^k, k = 0,1,\cdots,n+1$，得到$\varphi_{n+1}$，再解出它的零点就是 Gauss 点，然后将其代入$p_k(x), k = 0,1,\cdots,n$和求积公式即可确定$A_k.$这样仅需解一个代数方程和线性方程组即可求解答案。

参考资料

1. 黄云清, 《数值计算方法》, 科学出版社, 北京, 2012-06, ISBN 978-7-0302-3428-5.