辛普森数值积分公式的代数精度分析及复化误差估计

辛普森公式是数值积分中一个经典且高精度的近似计算方法，用于估算函数在某个区间上的定积分值。掌握它的精度特性与误差控制方法是高效应用的前提。

1. 理解单区间辛普森公式

核心公式与系数

对于一个连续函数 $f(x)$ 在区间 $[a, b]$ 上的定积分 $\int_a^b f(x) dx$，辛普森公式的近似值 $S(f)$ 由以下公式给出：

$$ S(f) = \frac{b-a}{6} \left[ f(a) + 4f\left(\frac{a+b}{2}\right) + f(b) \right] $$

其中，$a$ 和 $b$ 是积分下限和上限，$\frac{a+b}{2}$ 是区间的中点。公式中的系数 $1, 4, 1$ 是固定的权重。

代数精度是指一个数值积分公式能够精确积分所有次数不高于 $m$ 的多项式的最高次数 $m$。检验过程不需要复杂理论，按以下步骤操作即可。

确定测试函数序列：准备一组幂函数作为测试用例，依次是 $f_0(x)=1$, $f_1(x)=x$, $f_2(x)=x^2$, $f_3(x)=x^3$, $f_4(x)=x^4$。
计算精确积分值：在标准区间 $[0, 1]$ 上（这不会影响精度判断），计算每个函数的精确积分值 $I(f_k) = \int_0^1 f_k(x) dx$。
- 对于 $f_0(x)=1$，精确积分 $I(f_0) = 1$。
- 对于 $f_1(x)=x$，精确积分 $I(f_1) = 1/2$。
- 对于 $f_2(x)=x^2$，精确积分 $I(f_2) = 1/3$。
- 对于 $f_3(x)=x^3$，精确积分 $I(f_3) = 1/4$。
- 对于 $f_4(x)=x^4$，精确积分 $I(f_4) = 1/5$。
应用辛普森公式：将每个测试函数 $f_k(x)$ 和区间 $[0, 1]$ 代入辛普森公式 $S(f_k)$，其中 $a=0, b=1$。
- 计算 $S(f_0)$：$\frac{1}{6}[1 + 4*1 + 1] = 1$，与精确值 $1$ 相等。
- 计算 $S(f_1)$：$\frac{1}{6}[0 + 4*(1/2) + 1] = 1/2$，与精确值 $1/2$ 相等。
- 计算 $S(f_2)$：$\frac{1}{6}[0 + 4*(1/4) + 1] = 1/3$，与精确值 $1/3$ 相等。
- 计算 $S(f_3)$：$\frac{1}{6}[0 + 4*(1/8) + 1] = 1/4$，与精确值 $1/4$ 相等。
- 计算 $S(f_4)$：$\frac{1}{6}[0 + 4*(1/16) + 1] = 3/8$，精确值为 $1/5=0.2$，而 $3/8=0.375$，两者不相等。
判定结论：观察上述步骤。辛普森公式精确积分了所有次数不高于 $3$ 的多项式（即 $1, x, x^2, x^3$），但对 $x^4$ 失败。因此，辛普森公式的代数精度为 $3$。这意味着它能自动消除被积函数中 $1, x, x^2, x^3$ 这些部分产生的误差。

为了提高计算大范围积分时的精度，防止单个区间过大导致误差过大，我们采用“复化”或“分段”策略。

划分积分区间：将原始积分区间 $[a, b]$ 均匀地分成 $n$ 个子区间，$n$ 必须是一个偶数（这是辛普森公式的要求）。每个子区间的长度称为步长 $h$，计算公式为 $h = \frac{b-a}{n}$。那么 $n+1$ 个等距节点为 $x_k = a + k \cdot h$，其中 $k = 0, 1, 2, ..., n$。
应用基本公式到子区间：在每两个相邻节点 $x_{2k-2}, x_{2k}, x_{2k}$ 构成的子区间 $[x_{2k-2}, x_{2k}]$ 上应用基本的辛普森公式。注意，每次应用会覆盖两个子区间（$x_{2k-2}$ 到 $x_{2k}$），因此共需要 $n/2$ 次应用。
求和得到复化公式：将所有子区间上的辛普森近似值相加。整理系数后，得到复化辛普森公式的最终形式：

$$ S_n(f) = \frac{h}{3} \left[ f(a) + f(b) + 4 \sum_{k=1}^{n/2} f(x_{2k-1}) + 2 \sum_{k=1}^{n/2-1} f(x_{2k}) \right] $$

其中，$f(x_{2k-1})$ 是所有奇数编号节点的函数值，系数为 $4$；$f(x_{2k})$ 是所有内部偶数编号节点的函数值，系数为 $2$。端点 $a$ 和 $b$ 的系数为 $1$。

复化公式的误差与被积函数的光滑程度（即其高阶导数的性质）及步长 $h$ 有关。

$$ E_n(f) = I(f) - S_n(f) = -\frac{(b-a)}{180} h^4 f^{(4)}(\xi) $$

其中，$\xi$ 是区间 $[a, b]$ 内的某个未知点，$f^{(4)}(\xi)$ 是函数在该点的四阶导数。

解读误差公式：
- 误差与 $h^4$ 成正比。这意味着步长 $h$ 减半（即节点数 $n$ 翻倍），误差大约减小为原来的 $1/16$。这是四阶收敛的特性。
- 误差与被积函数在积分区间内的四阶导数的绝对值最大值有关。如果函数本身是三次多项式，其四阶导数 $f^{(4)}(x)=0$，则误差 $E_n(f)=0$，这与代数精度为 $3$ 的结论完全一致。
- 系数 $-\frac{b-a}{180}$ 是一个常数因子。
实用误差估计（事后误差估计）：在实际计算中，$\xi$ 和 $f^{(4)}(\xi)$ 未知。一个常用的方法是理查森外推法的思想。我们可以用不同步长（如 $h$ 和 $h/2$）计算两个近似值 $S_n(f)$ 和 $S_{2n}(f)$，然后利用误差公式的形式来估计误差和更高精度的近似值。
- 由于误差与 $h^4$ 成正比，我们有近似关系：$I(f) - S_n \approx C h^4$ 和 $I(f) - S_{2n} \approx C (h/2)^4 = C h^4 / 16$。
- 将两个方程相减消去 $I(f)$，可以解出 $I(f)$ 的一个更优估计，并估计误差。更简单的做法是，比较两次计算结果：
  $$ \frac{|S_{2n} - S_n|}{|S_n - S_{2n}|} \approx 15 $$
  这个比值可以帮助判断计算是否已经达到足够精度。如果两次计算结果非常接近（即该比值远大于 $15$），则说明 $S_{2n}$ 是一个相当精确的结果。

在实际编程中，可以编写一个函数来计算复化辛普森积分并控制误差。

编写核心计算函数：创建一个函数，输入参数包括被积函数 func，积分下限 a，上限 b，以及子区间数量 n（偶数）。函数内部根据复化辛普森公式进行求和计算并返回结果 $S_n$。
实现误差控制循环：在主程序中，设置一个初始偶数 $n$（如 $n=4$）和一个期望的误差容限 tol（如 1e-6）。进入循环：
- 调用核心函数计算当前 $n$ 下的积分值 $S_{old}$。
- 将 $n$ 加倍（n = 2 * n），再次计算积分值 $S_{new}$。
- 计算两次结果的绝对误差估计 error = abs(S_new - S_old)。
- 判断误差 error 是否小于容限 tol。如果是，则认为精度满足要求，跳出循环，以 $S_{new}$ 作为最终结果。否则，继续循环加倍 $n$。
注意事项：此方法假定了误差随 $h^4$ 单调递减。对于非常“粗糙”或高阶导数变化剧烈的函数，可能需要调整策略或接受更粗糙的精度控制。节点数量 $n$ 不能无限增加，需考虑计算资源和舍入误差的限制。