Gauss-Legendre 和 Gauss-Lobatto 积分表计算
在数值计算的编程实践中经常需要获取 Gauss-Legendre 和 Gauss-Lobatto 积分点和权重， 除了直接打表，还可以利用 Legendre 多项式自身的特点和牛顿迭代法通过高效的数值计算获得相应的积分点和权重，这也是本文关注的内容。
本文的主要动机是对下面两份 MATLAB 代码的学习和解释：
基本概念
Gauss 型数值积分
首先需要介绍一些关于 Gauss 型数值积分的基本内容。 对于标准积分区间 \([-1,1]\) 的如下形式的数值积分公式： \[ I(f) = \int_{-1}^1 f(x)\,dx \approx I_n(f) = \sum_{i=1}^n w_i f(x_i) \]
其中节点 \(x_i \in [-1,1]\)。这些节点和权重的选取有一共 \(2n\) 个自由度，至多可以达到 \(2n-1\) 阶的代数精度。如果加上某些额外的约束条件，则可达到的最优代数精度也要相应降低。
对于一般的有限区间 \([a,b]\) 上的数值积分，通过变量代换可以转化为 \([-1,1]\) 上的数值积分（记 \(t_i =\frac{a+b}2 + \frac{b-a}2 x_i\)）
\[ \begin{aligned} \int_a^b f(t)w(t)\,dt &\approx \frac{b-a}2 \sum_{i=1}^n w_i f(t_i) \end{aligned} \]
我们关注的是如何快速计算出标准积分区间 \([-1,1]\) 上的节点 \(\{x_i\}\) 和权重 \(\{w_i\}\)。 这些值的计算通常需要涉及到一元多项式的求根，只有对较小的 \(n\) 才可以给出解析表达式， 对于足够大的 \(n\)，我们只能希望在浮点数误差意义下获得足够精确的数值结果。
Legendre 正交多项式
在积分区间 \([-1,1]\) 和标准内积定义下，如下递归定义的正交多项式序列称为 Legendre 多项式： \[ \begin{aligned} P_0(x) &= 1\\ P_1(x) &= x\\ P_2(x) &= \frac12(3x^2 - 1)\\ P_3(x) &= \frac12(5x^3 - 3x) \end{aligned} \]
后续的递推关系为 \[ P_{n+1}(x) = \frac{2n+1}{n+1}\, x\,P_n(x) - \frac{n}{n+1}\, P_{n-1}(x) \]
Legendre 多项式满足单位正交性
\[ \int_{-1}^1 P_m(x) P_n(x)\,dx = \frac{2}{2n+1} \delta_{m\,n} \]
关于 Legendre 多项式有很多不同的定义，有的定义与上述定义的区别是给每一个多项式都乘以一个常数以确保均为首一多项式，有的定义则是乘以一个常数确保多项式的模为1，这些没什么实质影响。
Legendre 多项式有非常多的特殊性质，例如：
对于奇数 \(n\)，\(P_n(x)\) 是一个奇函数；对于偶数 \(n\)，\(P_n(x)\) 是一个偶函数。
\(P_n\) 在两个端点处的值具有如下规律
\[ P_n(1) = 1,\quad P_n(-1) = (-1)^n \]
Legendre 多项式的导函数在 \((-1,1)\) 满足如下关系
\[ \frac{d}{dx}P_n(x) = \frac{n}{x^2-1} \Big[x P_n(x) - P_{n-1}(x)\Big] \]
这些性质可以直接用递推关系归纳证明。
理论可以保证 Legendre 多项式 \(P_n(x)\) 具有 \(n\) 个互异的实根，并且全部分布于 \((-1,1)\) 范围内，这个性质直接保证了 Gauss-Legendre 积分公式的存在性。
由于 Legendre 多项式具有特殊的对称性（偶函数和奇函数交错），这导致下面各种节点和权重的计算也保持了明显的对称性：节点序列关于原点对称，对称节点的权重对应相等。
考虑如下特征值问题 \[ \frac{d}{dx}\left((1-x^2) \frac{d}{dx}\right) u = - \lambda u \] 它的特征值为 \(n(n+1)\)，对应的特征向量为 Legendre 多项式 \(P_n(x)\)，\(n=0,1,2,\cdots\)。
在下文中，Gauss 型积分公式的节点和权重的计算其实都可以归结为求特殊多项式的所有零点。本文中的做法是直接猜测一组合适的初值，然后通过牛顿迭代求出所有根。 除此之外，还有一类数值算法（The Golub-Welsch algorithm），它的思路是利用多项式的递推关系构造一个对应的方阵，使得目标多项式的所有零点恰为这个方阵的全部特征值，进而可以利用求特征值的各种数值算法进行求解。
Gauss-Legendre 积分
如果通过调整所有节点和权重，使代数精度达到最优的 \(2n-1\) 阶，就可以得到 Gauss-Legendre 积分公式，例如 \(n=2\) 和 \(n=3\) 的具体公式如下：
\[ \begin{aligned} \int_{-1}^1 f(x)\,dx &\approx f(-\frac{\sqrt{3}}3) + f(\frac{\sqrt3}3) \\ \int_{-1}^1 f(x)\,dx &\approx \frac59 f(-\sqrt{\frac35}) + \frac89 f(0) + \frac59f(\sqrt{\frac35}) \end{aligned} \]
理论推导易得，对于使用 \(n\) 个节点的 Gauss-Legendre 积分， \(n\) 个节点恰为 Legendre 多项式 \(P_n(x)\) 的 \(n\) 个实根，对应的权重 \(w_i\) 可以使用如下公式计算
\[ w_i = \frac{2}{(1-x_i^2)[P_n'(x_i)]^2} \]
这个公式的推导参考 Hildebrand, F. B. Introduction to Numerical Analysis. New York: McGraw-Hill, pp. 323-325, 1956.
我们使用牛顿法求解 \(P_n(x)\) 的所有根，使用如下公式计算导函数值 \[ P_n'(x) = \frac{n}{x^2-1} \Big[x P_n(x) - P_{n-1}(x) \Big] \] 更新公式为 \[ x^{new} = x - \frac{P_n(x)}{P_n'(x)} = x - \left(\frac{x^2-1}{n}\right) \frac{P_n(x)}{x P_n(x) - P_{n-1}(x)} \]
采用的初值为 Chebyshev 节点加上一组小扰动，Chebyshev 节点为 \[ y_i = \cos(\frac{2i-1}{2n} \pi),\quad i=1,\cdots,n \] 加上的小扰动为 \[ \delta_i = \frac{0.27}n \sin(\frac{n-1}{n+1}\left(\frac{2i-n-1}{n-1} \right)),\quad i=1,\cdots,n \]
Gauss-Lobatto 积分
与 Gauss-Legendre 积分不同，我们强制要求积分节点包括两个端点，即 \(x_1 = -1\)，\(x_{n} = 1\)， 调整剩余的 \(2n-2\) 个自由度，可以达到最优的 \(2n-3\) 阶的代数精度，这样就得到了 Gauss-Lobatto 积分公式， 例如 \(n=2\)，\(n=3\) 和 \(n=4\) 的具体公式如下：
\[ \begin{aligned} \int_{-1}^1 f(x)\,dx &\approx f(-1) + f(1)\\ \int_{-1}^1 f(x)\,dx &\approx \frac13 f(-1) + \frac43 f(0) + \frac13 f(1)\\ \int_{-1}^1 f(x)\,dx &\approx \frac16 f(-1) + \frac56 f(-\sqrt{\frac15}) + \frac56 f(\sqrt{\frac15}) +\frac16 f(1) \end{aligned} \]
注：
- 不含内部节点的 Gauss-Lobatto 积分公式其实就是梯形公式；
- 只含一个内部节点的 Gauss-Lobatto 积分公式其实就是著名的 Simpson's rule。
理论推导可得，对于使用 \(n\) 个节点的 Gauss-Lobatto 积分， \(n\) 个节点为如下 \(n\) 次多项式的所有根 \[ Q_n(x) := (1-x^2) P_{n-1}'(x) \] 对应的权重表达式为 \[ w_i = \frac{2}{n(n-1) [P_{n-1}(x_i)]^2}, \quad i=1,\dots,n \]
Gauss-Lobatto 积分的内部节点为 \(P_{n-1}'(x)\) 的所有 \(n-2\) 个实数根。
使用牛顿法求解 \(Q_n(x)\) 的所有根时，我们需要计算 \(Q_n(x)\) 以及它的导数 \(Q_n'(x)\)，可以利用如下公式快速计算：
- 函数值：注意到导函数满足如下公式 \[ \frac{d}{dx}P_{n-1}(x) = \frac{n-1}{x^2-1} \Big[x P_{n-1}(x) - P_{n-2}(x)\Big] \] 因此 \[ Q_n(x) = -(n-1) \Big[x P_{n-1}(x) - P_{n-2}(x)\Big] \]
- 导数值：注意到根据前面的特征值问题，恰有 \[ \frac{d}{dx}Q_n(x) = - (n-1) n P_{n-1}(x) \]
因此牛顿迭代的更新公式可以整理为 \[ x^{new} = x - \frac{Q_n(x)}{Q_n'(x)} = x - \frac{x P_{n-1}(x) - P_{n-2}(x)}{n P_{n-1}(x)} \]
牛顿迭代采用的初值为 Chebyshev-Gauss-Lobatto 节点 \[ y_i = \cos(\frac{i-1}{n-1} \pi),\quad i=1,\cdots,n \]
虽然前面的导函数关系对于 \(x=\pm1\) 无意义，但是这里 \(x=\pm1\) 就是 \(Q_n(x)\) 的根，初值已经直接包含了 \(\pm1\)，上述牛顿迭代的公式对于 \(x=\pm1\) 也保持不变，因此没有问题。
补充
编程实现
计算 Gauss-Legendre 和 Gauss-Lobatto 数值积分表的 MATLAB 源文件放置在 Github 仓库：gaussquad。 主要是基于两份参考代码（Legendre-Gauss Quadrature Weights and Nodes，Legende-Gauss-Lobatto nodes and weights）进行了一些修改和完善，包括增加了一些注释，优化了代码结构。
顺手使用 C++
实现了这两个数值算法，并且提供了运行期运算和编译期运算的两种版本：
前者返回 std::vector，后者返回 std::array。
为了实现 constexpr，后者的算法中使用了支持
constexpr 的简化版三角函数。 完整的源文件放置在 Github
仓库：gaussquadcpp，使用
C++20 标准。
由于浮点数的输出格式和自动舍入策略的不同（以及简化版三角函数的影响），C++ 实现的两个版本和 MATLAB 版本得到的数值结果（均保留15位精度）并不是完全一致的，观察到存在如下差异：
1vs0.9999999999999990vs-1.34940133673351e-79
这些差异显然在浮点数运算中可以忽略不计。
数值积分表
Gauss-Legendre 数值积分表：（\(n=2,\dots,8\)）
1 | n = 2 |
Gauss-Lobatto 数值积分表：（\(n=2,\dots,8\)）
1 | n = 2 |
