Gauss数值积分

关于 Gauss 数值积分公式的推导以及相关性质的笔记，注意并不是与正态分布相关的高斯积分 \(\int e^{-x^2}\,dx\)，在 wiki 上 Gauss 数值积分公式又称为 Gauss 求积公式。

Gauss 积分推导

最优代数精度

对于标准积分区间 \([-1,1]\) 的数值积分公式，节点和系数的选取有 2n+2 个自由度，至多可以达到 2n+1 阶的代数精度。

\[ I(f) = \int_{-1}^1 f(x)\,dx \approx I_n(f) = \sum_{i=0}^n A_i f(x_i) \]

断言：无论节点和系数如何选取，上述形式的数值积分公式不可能具有 2n+2 阶的代数精度。

证明断言：构造 2n+2 次多项式 \(p(x)\)

\[ p(x) = \Pi_{i=0}^n (x-x_i)^2 \ge 0 \]

显然有

\[ \begin{aligned} 0 < \int_{-1}^1 p(x)\,dx = I(p) \neq I_n(p) = \sum_{i=0}^n A_i p(x_i) = 0 \end{aligned} \]

因此不可能具有 2n+2 阶代数精度。

那么，是否可以选取合适的 \(\{A_i\},\{x_i\}\) 使得代数精度达到 2n+1 呢？答案是肯定的，我们将这种调整所有节点和系数，以达到最优代数精度的数值积分称为 Gauss 积分。

正向推导 Gauss 积分

考虑由这些节点组成的 n+1 次多项式 \(q(x)\)

\[ q(x) = \Pi_{i=0}^n (x-x_i) \]

对于任意不超过 2n+1 次的多项式 \(f(x)\)，都可以通过 \(q(x)\) 进行带余除法，得到下式

\[ f(x) = p(x)q(x) + r(x) \]

其中商式 \(p(x)\) 和余式 \(r(x)\) 均为至多 n 次的多项式，那么代入数值积分公式得到

\[ \begin{aligned} I(f) &= \int_{-1}^1 f(x)\,dx = \int_{-1}^1 p(x)q(x)\,dx + \int_{-1}^1 r(x)\,dx\\ &= I(pq) + I(r)\\ I_n(f) &= \sum_{i=0}^n A_i f(x_i) = \sum_{i=0}^n A_i\left[p(x_i)q(x_i)+r(x_i)\right] = \sum_{i=0}^n A_i r(x_i)\\ &= I_n(r) \end{aligned} \]

如果我们要求数值积分达到 2n+1 阶代数精度： \(I(f)=I_n(f)\)，那么必须满足两个条件：

第一个条件，因为 \(I_n(f)=I_n(r)\) 不再和 \(p,q\) 相关，我们应当确保 \(I(pq)=0\)；
第二个条件，要求 \(I(r)=I_n(r)\)，即对于至多 n 次的多项式需要精确成立。

第一个条件说明，\(q(x)\) 是和任意 \(P^n\) 多项式都正交的一个特殊的 n+1 次多项式

\[ \begin{aligned} (p,q)_{L^2} &= \int_{-1}^1 p(x)q(x)\,dx = 0\,&\forall p(x)\in P^n \end{aligned} \]

这个条件不足以唯一确定 \(q(x)\)，因为 \(k q(x)\) 也满足要求，这里 \(k\) 取任意非零实数。但是可以唯一确定 \(q(x)\) 的所有零点 \(\{x_i\}\)。（我们还可以保证 \(q(x)\) 不含重根，所有 n+1 个根都落在 \((-1,1)\) 内部，见下文正交多项式）从而我们可以确定数值积分公式的所有节点。

第二个条件，我们可以针对上述 n+1 个节点进行拉格朗日插值：对于不超过 n 次的多项式 \(f(x)\)，插值得到的 n 次多项式等于自身

\[ f(x) = \sum_{i=0}^n f(x_i) \frac{\Pi_{j=0}^n(x-x_j)}{\Pi_{j\neq i}(x_i-x_j)} \]

因此

\[ I(f) = \int_{-1}^1 f(x) = \sum_{i=0}^n f(x_i) \int_{-1}^1 \frac{\Pi_{j=0}^n(x-x_j)}{\Pi_{j\neq i}(x_i-x_j)}\,dx \]

取积分系数 \(A_i\) 为

\[ A_i = \int_{-1}^1 \frac{\Pi_{j=0}^n(x-x_j)}{\Pi_{j\neq i}(x_i-x_j)}\,dx \]

那么我们就可以得到一个完整的在积分区间 \([-1,1]\) 的 Gauss 积分公式

\[ \int_{-1}^1 f(x)\,dx \approx \sum_{i=0}^n A_i f(x_i) \]

其中积分节点 \(\{x_i,i=0,\dots,n\} \subset [-1,1]\) 满足

\[ \begin{aligned} \int_{-1}^1 p(x) \Pi_{i=0}^n (x-x_i)\,dx &= 0\,&\forall p(x) \in P^n \end{aligned} \]

系数为

\[ A_i = \int_{-1}^1 \frac{\Pi_{j=0}^n(x-x_j)}{\Pi_{j\neq i}(x_i-x_j)}\,dx \]

正交多项式

在积分区间 \([-1,1]\) 上，记 \(P^n\) 为不超过 n 次的多项式组成的线性空间，定义两个多项式的内积

\[ (f,g)_{L^2} = \int_{-1}^1 f(x)g(x)\,dx \]

称多项式 \(f\) 和 \(g\) 正交，如果 \((f,g)_{L^2} = 0\)。

我们可以构造 n+1 维空间 \(P^{n}\) 的一组正交基 \(\{g_i,i=0,\dots,n\}\)，称这些多项式为正交基，如果满足

\[ \begin{aligned} (g_i,g_j)_{L^2} &= 0, &i\neq j \end{aligned} \]

我们可以使用 Gram-Schmidt 正交化来逐渐得到 n+1 个两两正交的多项式。或者通过下列递推公式，可以更快地得到一组首一的正交多项式：

如下递归定义的多项式序列是两两正交的，\(p_n\) 为首一的 n 次多项式，进而前 n+1 个可以组成 \(P^n\) 的一组正交基。

\[ \begin{aligned} p_0(x) &= 1\\ p_1(x) &= x - a_1\\ p_n(x) &= (x-a_n)p_{n-1}(x) - b_n p_{n-2}(x)\,&n \ge 2 \\ &\left\{ \begin{aligned} a_n = (x p_{n-1},p_{n-1})_{L^2}/(p_{n-1},p_{n-1})_{L^2}\\ b_n = (x p_{n-1},p_{n-2})_{L^2}/(p_{n-2},p_{n-2})_{L^2} \end{aligned} \right. \end{aligned} \]

可以证明：如果 n+1 次非零多项式 \(p(x)\) 与任意不超过 n 次多项式都正交，那么 \(p(x)\) 所有零点落在 \((-1,1)\)，并且不含重根，这两个命题可以等价描述为：\(p(x)\) 在 \((-1,1)\) 变号 n+1 次。

反证：\(p(x)\) 若在 \((-1,1)\) 只变号 \(m\le n\) 次，在 \(\{t_i,i=1,\dots,m\}\) 变号，这里 \(t_i\) 单增。

\[ -1 = t_0 < t_1 < \cdots < t_m < t_{m+1} = 1 \]

那么定义 m 次多项式 \(g(x)\) 为

\[ g(x) = \Pi_{i=1}^m (x-t_i) \]

在 \((-1,1)\) 满足 \(f(x)g(x)\) 不变号，因此

\[ \int_{-1}^1 f(x)g(x)\,dx \neq 0 \]

但是 \(g(x)\) 是不超过 n 次的多项式，\(f\) 应当与 \(g\) 正交，得到矛盾，假设失败。

这个命题保证了我们可以通过正交多项式 \(p_{n+1}(x)\)，拿到所需的节点 \(\{x_i\}\)，它们全都落在 \((-1,1)\)，两两互异。

在积分区间 \([-1,1]\) 和标准内积定义下，递推公式得到的正交多项式序列称为勒让德多项式。勒让德多项式前几项为

\[ \begin{aligned} p_0(x) &= 1\\ p_1(x) &= x\\ p_2(x) &= x^2 -\frac13\\ p_3(x) &= x^3 - \frac35 x\\ \end{aligned} \]

勒让德多项式还有其它的等价定义，例如

\[ p_n(x) = \frac{1}{2^n n!}\frac{d^n}{dx^n}[(x^2-1)^n] \]

注意到 Legendre 多项式 \(p_2(x)\) 的零点为 \(\pm\frac{\sqrt3}3\)，\(p_3(x)\) 的零点为 \(0,\pm \sqrt{\frac35}\)，我们可以得到两个具体的 Gauss 积分公式 (对应为 n=1，n=2)

\[ \begin{aligned} p_2(x) = x^2 -\frac13\Rightarrow \int_{-1}^1 f(x)\,dx &\approx f(-\frac{\sqrt{3}}3) + f(\frac{\sqrt3}3)\\ p_3(x) = x^3 -\frac35 x \Rightarrow \int_{-1}^1 f(x)\,dx &\approx \frac59 f(-\sqrt{\frac35}) + \frac89 f(0) + \frac59f(\sqrt{\frac35})\\ \end{aligned} \]

Gauss 积分的性质

我们已经得到 \([-1,1]\) 的 Gauss 积分公式

\[ I(f) = \int_{-1}^1 f(x)\,dx \approx I_n(f) = \sum_{i=0}^n A_i f(x_i) \]

已知这个公式达到了最优的 2n+1 阶的代数精度，还有如下的性质：

\(\sum_{i=0}^n A_i=2\)，系数和为 2，直接取 \(f(x)=1\) 即可得到。
系数严格为正 \(A_i >0\)，可以取 \(f(x)=l_i(x)^2\) 为拉格朗日基函数的平方，此时积分大于零，数值积分结果为 \(A_i\)，数值积分精确成立，因此 \(A_i >0\)。
积分节点和系数都具有对称性：\(A_i = A_{n-i}\) ，\(x_i = -x_{n-i}\)，因为作代换 \(x\to -x\)，积分公式形式应保持不变。

带权 Gauss 积分

前文中我们考虑的是如下形式的积分和对应的数值积分

\[ I(f) = \int_{-1}^1 f(x)\,dx \]

现在我们考虑加权积分

\[ I^w(f) = \int_{-1}^1 f(x)w(x)\,dx \]

其中的权函数 \(w(x)\) 在 \([-1,1]\) 满足非负性

\[ w(x) \ge 0 \]

因此 \(I(f)\) 可以视作权函数 \(w(x)=1\) 的特例。

带权的 Gauss 积分公式形式为

\[ I^w(f) = \int_{-1}^1 f(x)w(x)\,dx \approx I_n^w(f) = \sum_{i=0}^n A_i f(x_i) \]

带权积分对应的带权内积定义为

\[ (f,g)^w_{L^2} = \int_{-1}^1 f(x)g(x)w(x)\,dx \]

与前文相同的分析可知，记 \(q(x) = \Pi_{i=0}^n(x-x_i)\)，那么 \(q(x)\) 需要满足

\[ \begin{aligned} (p,q)^w_{L^2} &= \int_{-1}^1 p(x)q(x)w(x)\,dx = 0\,&\forall p(x)\in P^n \end{aligned} \]

即 n+1 次多项式 \(q(x)\) 需要在带权内积的意义下，与任意至多 n 次多项式都正交。

带权正交多项式

与 \(w=1\) 的正交多项式一样，我们把内积取为带权内积，就可以得到不同带权内积意义下的正交多项式。（正交多项式在不同的定义下可能会差一个系数，未必保证首项系数为 1）

例如权函数取

\[ w(x) = \frac{1}{\sqrt{1-x^2}} \]

对应的正交多项式称为（第一类） Chebyshev 多项式

\[ p_0(x)=1,p_{n+1}(x) = \frac{1}{2^n}\cos((n+1)\arccos x), n \ge 0 \]

前几项依次为

\[ \begin{aligned} p_0(x) &= 1\\ p_1(x) &= x\\ p_2(x) &= x^2 -\frac12 \end{aligned} \]

Chebyshev 多项式 \(p_{n+1}(x)\) 的零点有显式公式

\[ x_i = \cos(\frac{2i+1}{2n+2}\pi),i=0,\dots,n \]

并且积分系数也很有特点，全都是一样的值

\[ A_i = \frac{\pi}{n+1} \]

即可得到 n+1 个点对应的 Gauss-Chebyshev 积分公式

\[ \int_{-1}^1 f(x) \frac{1}{\sqrt{1-x^2}}\,dx \approx \frac{\pi}{n+1} \sum_{i=0}^n f\left(\cos(\frac{2i+1}{2n+2}\pi)\right) \]

这个公式很有特点，除了显然的三角函数背景之外，积分系数一样意味着乘法计算量显著地减少了。

在权重 \(w(x)=1\) 时，正交多项式为 Legendre 多项式，我们也称此时的 Gauss 积分公式为 Gauss-Legendre 积分公式。

Gauss 积分的误差

对于 \([-1,1]\) 的带权 Gauss 积分公式

\[ I^w(f) = \int_{-1}^1 f(x)w(x)\,dx \approx I_n^w(f) = \sum_{i=0}^n A_i f(x_i) \]

数值积分的误差 \(\widetilde{R}\) 为

\[ \widetilde{R} = \int_{-1}^1 f(x)w(x)\,dx - \sum_{i=0}^n A_i f(x_i)=\frac{f^{(2n+2)}(\eta)}{(2n+2)!} \int_{-1}^1 \Pi_{i=0}^n (x-x_i)^2 w(x)\,dx \]

证明：对 \(f(x)\) 在 \(\{x_i\}\) 进行 Hermite 插值， 2n+2 个条件，得到一个 2n+1 次的插值多项式 \(g(x)\)，有插值误差

\[ f(x)-g(x) = R(x) = \frac{f^{(2n+2)}(\xi_x)}{(2n+2)!}\Pi_{i=0}^n (x-x_i)^2 \]

注意到 Gauss 积分公式对于 \(g(x)\) 是精确成立的，再利用插值性质可以得到如下等式

\[ \int_{-1}^1 g(x)w(x)\,dx = \sum_{i=0}^n A_i g(x_i)= \sum_{i=0}^n A_i f(x_i) \]

因此数值积分误差

\[ \begin{aligned} \widetilde{R} =& \int_{-1}^1 f(x)w(x)\,dx - \sum_{i=0}^n A_i f(x_i)\\ =& \int_{-1}^1 f(x)w(x)\,dx - \int_{-1}^1 g(x)w(x)\,dx \\ =& \int_{-1}^1 \frac{f^{(2n+2)}(\xi_x)}{(2n+2)!}\Pi_{i=0}^n (x-x_i)^2\,dx\\ =& \frac{f^{(2n+2)}(\eta)}{(2n+2)!}\int_{-1}^1 \Pi_{i=0}^n (x-x_i)^2\,dx \end{aligned} \]

证毕。

任意区间的 Gauss 积分

已知标准区间 \([-1,1]\) 上的带权 Gauss 积分公式

\[ I^w(f) = \int_{-1}^1 f(x)w(x)\,dx \approx I_n^w(f) = \sum_{i=0}^n A_i f(x_i) \]

对于一般区间 \([a,b]\) 的数值积分，可以通过变量代换转化为 \([-1,1]\) 的积分

\[ t = \frac{a+b}2 + \frac{b-a}2 x \]

从而利用 Gauss 积分公式

\[ \begin{aligned} \int_a^b f(t)w(t)\,dt &= \int_{-1}^1 f(\frac{a+b}2 + \frac{b-a}2 x)\, w(\frac{a+b}2 + \frac{b-a}2 x)\frac{b-a}2 \,dx\\ &\approx \frac{b-a}2 \sum_{i=0}^n A_i f(\frac{a+b}2 + \frac{b-a}2 x_i) \end{aligned} \]

只需要注意多了一个变换系数，即两个积分区间之比。