Lagrange多项式数值积分
工程中，有时需要我们把一个不断变化的量，在一个时间段或空间范围内累积起来求和，就需要用到积分。比如对速度的积分得到位移量。 根据牛顿莱布尼茨公式 \(\int_{a}^b f(x) dx = F(b) - F(a)\)，只要找到被积函数 \(f(x)\) 的原函数 \(F(x)\)，就可以轻松写出解析积分值。 但是很不幸，我们会面临数值微分相同的问题。 要么拿到的是一些离散的采样点，连被积函数 \(f(x)\) 的表达式都不知道，更不用说原函数 \(F(x)\) 了。要么 \(f(x)\) 的形式过于复杂，不知道怎么求原函数。 数值积分就是要解决这种离散数据或黑盒函数的累积问题。
数值积分通过对离散数据点的加权求和 \(\sum_{i=0}^n a_i f(x_i)\) 来近似计算连续的定积分。插值多项式是最经典的构造数值积分的方法。 本节我们借助拉格朗日多项式来初步探讨一下数值积分。 通过对它的截断误差分析，我们会发现相比于数值微分而言，数值积分是一种稳定的数值计算问题。
1. Lagrange多项式数值积分
从几何角度来看，定积分 \(\int_a^b f(x) dx\) 的几何意义是曲线与 \(x\) 轴围成的曲边梯形的净面积。 区间 \([a,b]\) 上，\(x\) 轴上方的面积为正，\(x\) 轴下方的面积为负，这些正负面积相加的结果，就是所谓的净面积。 尽管我们不知道被积函数 \(f(x)\) 的表达式，但是我们可以找一个多项式 \(P(x)\) 来逼近它。 只要 \(P(x)\) 能足够逼近 \(f(x)\) 那么它们与 \(x\) 轴围成的曲边梯形的净面积就会很接近。
我们可以选择 Larange 多项式来逼近 \(f(x)\)。如果点 \(\{x_0, \cdots, x_n\}\) 是闭区间 \([a, b]\) 上连续函数 \(f(x) \in C^{n+1}[a,b]\) 的 \(n+1\) 个不同的采样点。 记拉格朗日基函数 \(L_i(x) = \prod_{i \neq k}\frac{x - x_i}{x_k - x_i}\)，根据逼近误差公式，有:
$$ f(x) = \sum_{i = 0}^n f(x_i)L_i(x) + \frac{f^{(n+1)}(\varepsilon(x))}{(n+1)!} \prod_{i=0}^n(x - x_i) $$其中 \(\varepsilon(x) \in [a,b]\) 是一个关于 \(x\) 的数值。对 \(f(x)\) 在闭区间 \([a,b]\) 上积分有:
$$ \begin{equation}\label{f1} \int_a^b f(x) dx = \underbrace{\int_a^b \left(\sum_{i=0}^n f(x_i)L_i(x) \right) dx}_{Int(x)} + \underbrace{\frac{1}{(n+1)!}\int_a^b\left(\prod_{i=0}^n(x-x_i)f^{(n+1)}(\varepsilon(x)) \right) dx}_{E(f)} \end{equation} $$我们用上式中的第一项 \(Int(x)\) 来近似 \(f(x)\) 在闭区间 \([a,b]\) 上的积分，有:
$$ \begin{equation}\label{f2} \begin{aligned} \int_a^b f(x) dx & \approx \int_a^b \left(\sum_{i=0}^n f(x_i)L_i(x) \right) dx \\ & = \sum_{i=0}^n \left[ f(x_i) \cdot\int_a^b L_i(x) dx \right] \end{aligned} \end{equation} $$给定积分区间 \([a,b]\)，拉格朗日基函数的积分 \(\int_a^b L_i(x) dx\) 就是一个固定值，记为 \(a_i\)。那么上式可简写为下式。我们完全按照该公式实现了 LagrangePolynomialIntegral。
$$ \begin{equation}\label{f3} \int_a^b f(x) dx \approx \sum_{i=0}^n a_i f(x_i) \end{equation} $$式(\(\ref{f1}\))中的第二项，就是截断误差:
$$ \begin{equation}\label{f4} E(f) = \frac{1}{(n+1)!}\int_a^b\left(\prod_{i=0}^n(x-x_i)f^{(n+1)}(\varepsilon(x)) \right) dx \end{equation} $$我们将连续的积分运算，改成离散的加权求和。如果我们对区间 \([a,b]\) 做等距采样，就可以得到 Newton-Cotes 积分公式。
2. Newton-Cotes 积分公式
根据采样点是否包含区间端点\(a,b\)，Newton-Cotes 积分公式有闭式(Closed)和开式(open)两种，如下图1所示。一般情况下我们采用 Closed Newton-Cotes 公式， 但遇到被积函数在端点处没有定义或者存在奇异点的时候，会选择 Open Newton-Cotes 公式。
 |
 |
| 图1.a Closed Newton-Cotes 公式 | 图1.b Open Newton-Cotes 公式 |
如上图1.a 所示，\(n+1\) 个采样点可以写作 \(x_i = x_0 + ih, i \in \{0,1,\cdots,n\}\)，并且 \(x_0 = a, x_n = b, h = (b-a)/n\)。 那么按照式 (\(\ref{f3}\)) 写出近似积分就是闭式 Newton-Cotes 公式。如果 \(n\) 是偶数并且 \(f \in C^{n+2}[a,b]\)，那么它的截断误差为:
$$ \begin{equation}\label{f5} E(f) = \frac{h^{n+3}f^{(n+2)(\xi)}}{(n+2)!}\int_0^n t^2(t-1)\cdots(t-n) dt \end{equation} $$如果 \(n\) 是奇数，并且 \(f \in C^{n+1}[a,b]\)，那么它的截断误差为:
$$ \begin{equation}\label{f6} E(f) = \frac{h^{n+2}f^{(n+1)(\xi)}}{(n+1)!}\int_0^n t(t-1)\cdots(t-n) dt \end{equation} $$要证明上式(\(\ref{f5}\))(\(\ref{f6}\))，我们需要引入几个引理。首先令 \(x = x_0 + th\) 将截断误差(\(\ref{f4}\))写成关于 \(h\) 的参数形式，可以将积分区间从 \([a,b]\) 映射到 \([0, n]\) 上:
$$ \begin{aligned} E(f) & = \frac{1}{(n+1)!}\int_a^b\left(\prod_{i=0}^n((x_0+th)-(x_0+ih))f^{(n+1)}(\varepsilon(t)) \right) d(th) \\ & = \frac{h^{n+2}}{(n+1)!}\int_0^n\left(f^{(n+1)}(\varepsilon(t)) \prod_{i=0}^n (t-i) \right) d(t) \\ \end{aligned} $$引理1: 记 \(\psi_n(t) = \prod_{i = 0}^n(t-i)\)，取区间 \([0,n]\) 的中点 \(n/2\) 有:
$$ \psi_n(\frac{n}{2} + \alpha) = (-1)^{n+1}\psi_n(\frac{n}{2} - \alpha) $$证: 令正整数 \(j = n - i\)，有 \(\prod_{i = 0}^n(\frac{n}{2} + \alpha - i) = \prod_{j = 0}^n(\frac{n}{2} + \alpha - (n - j)) = \prod_{j = 0}^n-(\frac{n}{2} - \alpha - j) = (-1)^{n+1} \prod_{j = 0}^n(\frac{n}{2} - \alpha - j)\) □
引理2: 将区间 \([0, n]\) 用正数 \(i\) 分割成 \(n\) 个子区间，相邻的两个子区间中，函数 \(\psi_n(t)\) 的符号相反，并且:
- \(0 < t + 1 < \frac{n}{2}\) 时，\(|\psi_n(t+1)| < |\psi_n(t)|\)
- \(\frac{n}{2} < t < n\) 时，\(|\psi_n(t+1)| > |\psi_n(t)|\)
证: 根据 \(\psi_n(t)\) 的定义，\(t\) 在相邻两个子区间 \([i, i+1], [i+1, i+2]\)切换时，\(\prod_{i = 0}^n (t-i)\) 中只有第 \(i+1\) 项符号发生改变，所以 \(\psi_n(t)\) 的符号翻转。 此外，我们有:
$$ \psi_n(t) = t(t-1)(t-2)\cdots(t-n) $$ $$ \psi_n(t+1) = (t+1)t(t-1)\cdots(t-n+1) $$上述两式相除后整理得:
$$ \left| \psi_n(t+1) \right| = \left| \psi_n(t) \right| \cdot \left| \frac{t+1}{t-n} \right| $$当 \(0 < t+1 < \frac{n}{2}\) 时，有 \(n - t > n - (t + 1) > \frac{n}{2}\)，所以 \(\frac{|t+1|}{|t-n|} < 1\)，因此 \(| \psi_n(t+1) | < | \psi_n(t) |\)。
当 \(\frac{n}{2} < t < n\) 时，有 \(n - t < \frac{n}{2}\)，所以 \(\frac{|t+1|}{|t-n|} > 1\)，因此 \(| \psi_n(t+1) | < | \psi_n(t) |\)。□
引理3: 记函数 \(\Psi_n(x) = \int_0^x \psi_n(t) dt\)。那么:
- \(n\) 为偶数时 \(\Psi_n(0) = \Psi_n(n) = 0\)，并且 \(0 < x < n\) 时 \(\Psi_n(x) > 0\)。
- \(n\) 为奇数时 \(\Psi_n(0) = 0, \Psi_n(n) = 2\Psi_n(n/2)\)，并且 \(0 < x < n\) 时 \(\Psi_n(x) < 0\)。
证: 当 \(n\) 为偶数时，根据引理1，\(\psi_n(t)\) 关于 \(n/2\) 反对称，即 \(\psi_n(\frac{n}{2} + \alpha) = -\psi_n(\frac{n}{2} - \alpha)\)。 按照积分的几何意义，有 \(\Psi_n(0) = \Psi_n(n) = 0\)。当 \(0 < t < 1\) 时，有 \(t-0 > 0\)，以及 \(t-i < 0, i \in \{1,\cdots,n\}\)。 此时 \(\psi_n(t)\) 中一共有偶数项负数，一项正数。所以 \(\psi_n(t) > 0, \Psi_n(x) > 0\)。当 \(1 < t < 2\) 时，根据引理2知 \(\psi_n(t)<0\)， 并且 \(|\psi_n(t)| < |\psi_n(t-1)|\)，因此根据积分几何意义 \(\Psi_n(x)\) 中负面积一定小于正面积，所以此时 \(\Psi_n(x) > 0\)。 重复类似的推理，我们可以断言 \(0 < x < n/2\) 时 \(\Psi_n(x) > 0\)。又因为 \(n\) 为偶数时，\(\psi_n(t)\) 关于 \(n/2\) 反对称，可以推知 \(x > n/2\) 时也成立。
当 \(n\) 为奇数时，\(\psi_n(t)\) 关于 \(n/2\) 对称，即 \(\psi_n(\frac{n}{2} + \alpha) = -\psi_n(\frac{n}{2} - \alpha)\)。所以有 \(\Psi_n(0) = 0, \Psi_n(n) = 2\Psi_n(n/2)\)。 当 \(0 < t < 1\) 时，\(\psi_n(t)\) 中共有奇数项负数和一项正数，所以 \(\psi_n(t) < 0\)。根据引理2我们也可以断言 \(0 < x < n/2\) 时 \(\Psi_n(x) < 0\)， 根据此时 \(\psi_n(t)\) 关于 \(n/2\) 对称，所以 \(x > n/2\) 时，也有 \(\Psi_n(x) < 0\)。
当 \(n\) 为偶数时，根据分步积分法则有:
$$ \begin{aligned} E(f) & = \frac{h^{n+2}}{(n+1)!}\int_0^n f^{(n+1)}(\varepsilon(t)) \frac{d \Psi_n(t)}{dt} dt \\ & = \frac{h^{n+2}}{(n+1)!}\left(\left. f^{(n+1)}(\varepsilon(t)) \Psi_n(t) \right|_0^n - \int_0^n \Psi_n(t) \frac{d f^{(n+1)}(\varepsilon(t))}{dt} dt\right) \end{aligned} $$根据引理3，此时 \(\Psi_n(0) = \Psi_n(n) = 0\) 所以 \(\left. f^{(n+1)}(\varepsilon(t)) \Psi_n(t) \right|_0^n = 0\)，有:
$$ E(f) = \frac{h^{n+2}}{(n+1)!}\int_0^n \Psi_n(t) \frac{d f^{(n+1)}(\varepsilon(t))}{dt} dt = \frac{h^{n+2}}{(n+1)!}\int_0^n \Psi_n(t) f^{(n+2)}(\varepsilon(t))\varepsilon'(t) dt $$实际上证明到这里我就卡住了，我不知道系数 \(h\) 和 \(n+2\) 是怎么搞出来的，\(\varepsilon'(t)\) 又是怎么消失的。据说参考 Hermite 插值增加一条约束就可以搞定。 根据引理3还有 \(\Psi_n(x) ≥ 0\)，根据加权积分中值定理: 如果 \(f(x)\) 连续，\(g(x)\) 在\([a,b]\)上不变号，那么存在 \(\xi \in [a,b]\) 使得 \(\int_a^b f(x)g(x) dx = f(\xi)\int_a^b g(x) dx\)。有
$$ \int_0^n \Psi_n(t) f^{(n+2)}(\varepsilon(t)) dt = f^{(n+2)}(\xi)\int_0^n \Psi_n(t) dt, \qquad \xi \in [0, n] $$再次使用分步积分法则有:
$$ \begin{aligned} \int_0^n \Psi_n(t) dt & = \left. t \Psi_n(t) \right|_0^n - \int_0^n t \frac{\Psi_n(t)}{dt} dt \\ & = - \int_0^n x \psi_n(t) dt \\ & = - \int_0^n t^2(t-1) \cdots (t-n) dt \end{aligned} $$开式 Newton-Cotes 公式以 \(x_0 = a + h\) 作为第一个采样点，最后一个采样点为 \(x_n = b - h\)，\(n+1\) 个采样点仍然写作 \(x_i = x_0 + ih, i \in \{0,\cdots,n\}\)。 我们把区间的两端标记为 \(x_{-1} = a, x_{n+1} = b\)，如上图1.b 所示，采样步长 \(h = (b-a)/(n+2)\)。如果 \(n\) 是偶数并且 \(f \in C^{n+2}[a,b]\)，那么它的截断误差为:
$$ \begin{equation}\label{f7} E(f) = \frac{h^{n+3}f^{(n+2)(\xi)}}{(n+2)!}\int_{-1}^{n+1} t^2(t-1)\cdots(t-n) dt \end{equation} $$如果 \(n\) 是奇数，并且 \(f \in C^{n+1}[a,b]\)，那么它的截断误差为:
$$ \begin{equation}\label{f8} E(f) = \frac{h^{n+2}f^{(n+1)(\xi)}}{(n+1)!}\int_{-1}^{n+1} t(t-1)\cdots(t-n) dt \end{equation} $$从上述截断误差的表达式上看，当 \(n\) 为偶数时算法的代数精度可以达到 \(n + 1\)。奇数时则只有 \(n\)。代数精度(Degree of Accuracy / Precision) 是一种描述数值积分公式精度的方法。 如果一个数值积分公式能够对所有次数不超过 \(n\) 的多项式都能做到 0 误差的完美计算，而对 \(n + 1\) 次多项式存在误差，那么我们说该公式具有 \(n\) 阶代数精度。
从代数精度上来看，理论上采样点越多积分公式就越精准。但是在讨论Lagrange多项式的时候，我们就已经提过一种龙格现象。 这种等距采样的高阶的 Lagrange 多项式会在区间的两端剧烈震荡，并不能很好的拟合目标函数 \(f(x)\)，所以实际的数值积分结果也很差。为了克服这一问题， 人们会用低阶的分段多项式来作插值，或者放弃等距采样。我们将在本部分的其它文章中深入探讨。
3. 梯形法则和 Simpson 法则
根据式(\(\ref{f1}\))我们最常见的就是梯形法则(Trapezoidal)和 Simpson 法则。它们都是闭式的 Newton-Cotes，其中梯形法则对区间的两个端点采样，Simpson 则对两个端点和区间的中点采样。 它们按照式(\(\ref{f1}\))展开化简之后就是:
$$ \begin{aligned} \text{梯形法则}: & \int_a^b f(x)dx = \frac{h}{2} \left[f(x_0) + f(x_1)\right] - \frac{h^3}{12}f''(\xi) \\ \text{Simpson法则}: & \int_a^b f(x)dx = \frac{h}{3} \left[f(x_0) + 4f(x_1) + f(x_2)\right] - \frac{h^5}{90}f^{(4)}(\xi) \end{aligned} $$从这两个简化之后的表达是中可以看到，这里是在乘以采样步长 \(h\)。 相比于数值微分的除以 \(h\) 而言，这是一种很稳定的形式。 这样缩小采样步长降低截断误差的同时，不会因为除以一个小数而把舍入误差放大。
 |
 |
| 图2.a 梯形积分法则 | 图2.b Simpson积分法则 |
显然，梯形法则和 Simpson 法则的截断误差是比较大的，在实际应用中通常采用增加采样点，以分段的方式进行数值积分。 我们为此实现了CompositeTrapeRule, CompositeSimpsonRule。