高斯积分

上一节中我们介绍了 Romberg 积分，它通过 Richardson 外展显著提高了梯形法则的积分精度。它的内核仍然固定采样步长 $h$ 的Newton-Cotes 积分公式。从积分的几何意义上来考虑这个固定的采样步长，显然是不太精确的。如果采样点数目一样，那么采样点的分布对于最后的计算结果会有很大影响。高斯研究出了一系列的正交的多项式 $\{P_n(x)\}$，并且发现对于 $n$ 采样点，如果参考 $P_n(x) = 0$ 的根选择采样点，就会得到 $2n-1$ 的代数精度。而 Newton-Cotes 积分公式的代数精度只有 $n+1$。后来人们发现高斯的这组多项式，正是物理界及其重要的勒让德多项式。

本文，我们先简要介绍勒让德多项式的性质，再证明高斯积分的 $2n-1$ 的代数精度，最后给出可以自由指定阶数 $n$ 的 Gaussian-Legendre 积分实现。

1. 勒让德多项式

最初勒让德为了计算天体的引力势，引入了一组多项式 $\{P_0(x), P_1(x), \cdots, P_n(x), \cdots \}$，其递推公式如下: $$ \begin{equation}\label{f1} (n+1)P_{n+1}(x) = (2n+1)xP_n(x) - nP_{n-1}(x) \end{equation} $$ 前两项为 $P_0(x) = 1, P_1(x) = x$，参见例程。后来人们把这组多项式命名为Legendre 多项式，它们有一种被称为正交性的特性，使得任何次数 $m < n$ 的多项式 $Q_m(x)$，与 $P_n(x)$ 的乘积在区间 $[-1, 1]$ 上的积分为 0，即:

$$ \begin{equation}\label{f2} \int_{-1}^1 Q_m(x)P_n(x) dx = 0 \end{equation} $$

此外，$n$ 阶勒让德多项式 $P_n(x)$ 具有对称性 $P_n(-x) = (-1)^nP_n(x)$。在开区间 $(-1, 1)$ 内恰好有 $n$ 个互不相同的实零点，这些零点关于 $x = 0$ 对称分布。根据这些特性，我们可以写出区间 $[-1, 1]$ 上的高斯积分。

2. 区间 $[-1,1]$ 上的高斯积分

若 $x_1, x_2, \cdots, x_n$ 是 $n$ 阶勒让德多项式 $P_n(x)$ 的零点，那么以这些零点作为采样点，在区间 $[-1, 1]$ 上对函数 $f(x)$ 积分，根据Lagrange多项式数值积分有:

$$ \begin{equation}\label{f3} \int_{-1}^1 f(x) dx = \sum_{i=1}^n c_i f(x_i) + \frac{1}{n!}\int_{-1}^1 \left(\prod_{i=1}^n(x-x_i)f^{(n)}(\varepsilon(x)) \right) dx \end{equation} $$

其中 $c_i, i = \{1,2,\cdots, n\}$ 为 Lagrange 基函数在 $[-1,1]$ 上的定积分:

$$ \begin{equation}\label{f4} c_i = \int_{-1}^1 \sum_{j=1,j\neq i}^n \frac{x - x_j}{x_i - x_j} dx \end{equation} $$

如果 $f(x)$ 是一个次数小于 $2n$ 的多项式，那么($\ref{f3}$)的截断误差项为 0。也就是说，以 $n$ 阶勒让德多项式 $P_n(x)$ 的零点作为采样点的数值积分的代数精度能达到 $2n-1$。 首先，如果 $f(x)$ 是一个次数小于 n 的多项式，那么其 $n$ 阶导数 $f^(n)(x) = 0$，显然其截断误差为 0。

其次，如果 $f(x)$ 的次数大于等于 $n$ 但是小于 $2n$，那么根据多项式的除法，有:

$$ \begin{equation}\label{f5} f(x) = Q(x)P_n(x) + R(x) \end{equation} $$

上式中 $Q(x), R(x)$ 是两个次数小于 $n$ 的多项式。那么根据勒让德多项式的特性 ($\ref{f2}$) 在区间 $[-1, 1]$ 上的积分有:

$$ \int_{-1}^{1} f(x) dx = \int_{-1}^{1} Q(x)P_n(x) dx + \int_{-1}^{1} R(x) dx = \int_{-1}^{1} R(x) dx $$

上式表明在区间 $[-1, 1]$ 上，$f(x)$ 与 $R(x)$ 的积分相等。由于 $x_i$ 是 $P_n(x)$ 的零点，所以式 ($\ref{f5}$) 可推出 $f(x_i) = R(x_i)$，有:

$$ \int_{-1}^{1} R(x) dx = \sum_{i=1}^n c_i R(x_i) = \sum_{i=1}^n c_i f(x_i) $$

所以，如果 $f(x)$ 是一个次数小于 $2n$ 的多项式，那么一定有关系:

$$ \begin{equation}\label{f6} \int_{-1}^{1} f(x) dx = \sum_{i=1}^n c_i f(x_i) \end{equation} $$

对于任意的积分区间 $[a,b]$ 我们可以定义如下的关系，将积分区间映射到 $[-1, 1]$ 上:

$$ \begin{equation}\label{f7} t = \frac{2x - a - b}{b - a} \Longleftrightarrow x = \frac{1}{2}[(b-a)t + a + b] \end{equation} $$

如此，任意区间的积分都可以转换成 $[-1, 1]$ 上的积分:

$$ \begin{equation}\label{f8} \int_a^b f(x)dx = \int_{-1}^{1} f\left(\frac{(b - a)t + a + b}{2}\right) \frac{(b-a)}{2}dt \end{equation} $$

3. GaussianIntegration 实现

对于式 ($\ref{f4}$) 而言，确定了勒让德多项式的阶数 $n$，就确定了积分权重 $c_i$。所以在实际应用中，我们完全可以先计算一个采样点和积分权重的表。这样给定任意一个被积函数 $f(x)$，我们都可以直接根据式 ($\ref{f6}$)($\ref{f8}$) 通过加权求和完成数值积分。我们定义了类 GaussianIntegration，在其构造函数中完成积分表的计算。

展开 GaussianIntegration 的构造函数 >>>

        GaussianIntegration(int n) : _n(n)
        {
            assert(_n > 1);

            _roots.resize(_n);
            _weights.resize(_n);
            int m = (_n + 1) / 2;

            auto P = LegendrePolynomial<Scalar>(_n);
            for (int i = 0; i < m; ++i) {
                // 契比雪夫初值估计
                Scalar z = std::cos(static_cast<Scalar>(M_PI) * (i + 0.75) / (_n + 0.5));
                Scalar r = P.NewtonRaphson(z);

                Scalar pxi, dpxi;
                P.Horner(r, pxi, dpxi);

                _roots[i] = -r;
                _roots[_n - 1 - i] = r;

                // 权重公式
                Scalar w = 2.0 / ((1.0 - r * r) * dpxi * dpxi);
                _weights[i] = w;
                _weights[_n - 1 - i] = w;
            }
        }

构造函数的实现，如上面代码片段所示。其中，关于勒让德多项式求根和积分权重 $c_i$ 的计算，两个技术细节需要注意一下。这里我们直接按照教科书中说法给出结论，以后专门研究勒让德多项式时会再来看这两点。

高斯积分以勒让德多项式的零点为采样点，这涉及到了多项式求根的问题。虽然通过Müller算法可以解出多项式的所有根，但该方法涉及到综合除法，计算复杂而且数值稳定性并不太好。实际上如果我们能给出一个比较好的初始估计，就可以直接使用 NewtonRaphson 算法了。上面代码片段中，采用的契比雪夫初值估计公式 $$ \cos \left(\pi \frac{i + 0.75}{n + 0.5}\right),\qquad i = 0, 1, \cdots, \left\lfloor \frac{n-1}{2} \right\rfloor $$ 利用三角函数，在 $[-1, 1]$ 区间内估算所有零点的大致位置。这个估计值和真正的勒让德零点极其接近，据说其估计误差是 $O(1/n^2)$ 级别的，可以保证牛顿迭代的收敛性。计算出采样点 $x_{i}, i \in \{1, 2, \cdots, n\}$ 之后，可以直接通过采样点导数值得到积分权重 $c_i$:

$$ c_i = \frac{2}{(1-x_{i}^{2})[P_{n}'(x_{i})]^{2}} $$

完成采样点和积分权重的计算之后，给定任意函数都可以直接通过成员函数 Evaluate 按照式 ($\ref{f8}$) 写出任意区间 $[a,b]$ 上的积分。

4. 完

最后需要强调一点。虽然理论上提高了勒让德多项式的阶数，就可以增加采样点，进而获得更高精度的高斯积分。但实际工程中，人们并不会使用极高阶的高斯积分。仍然采用分段的方式，在每个小区间上使用一个低阶的高斯积分，然后求和。

高斯积分

1. 勒让德多项式

2. 区间 \([-1,1]\) 上的高斯积分

3. GaussianIntegration 实现

4. 完