数值分析讲义（二）：数值积分
建议先阅读 数值分析讲义（一）：插值方法，因为本章的 Newton-Cotes 求积会直接用到插值多项式和 Lagrange 基函数。本文接着讨论数值积分，也就是如何用有限个函数值近似计算定积分。
数值积分
本章讨论如何近似计算定积分
\[\int_a^b f(x)\,dx\]的常用方法。
积分任务
给定可积函数 $f:[a,b]\to\mathbb{R}$，计算
即使函数形式并不复杂，也未必能解析地写出原函数，例如
\[\frac{\sin x}{x} \qquad\text{和}\qquad e^{-x^2}.\]这时就需要数值积分。许多数值积分方法都基于同一个想法：先在若干支撑点
\[(x_i,f(x_i)),\qquad x_i\in[a,b],\]上用多项式 $p_n$ 近似 $f$，再对这个多项式积分。这样得到的近似为
\[I_n(f)=\int_a^b p_n(x)\,dx,\]这类方法称为插值型求积。
2.1 Newton-Cotes 求积
2.1.1 闭型 Newton-Cotes 求积
对 $n\in\mathbb{N}$，选择等距支撑点
\[x_i=a+ih,\qquad i=0,\ldots,n,\qquad h=\frac{b-a}{n}.\]于是插值多项式 $p_n$ 的 Lagrange 表示为
\[p_n(x)=\sum_{i=0}^{n} f(x_i)L_{i,n}(x), \qquad L_{i,n}(x)= \prod_{\substack{j=0\\ j\ne i}}^{n} \frac{x-x_j}{x_i-x_j}.\]由此得到数值求积公式
\[I_n(f) =\int_a^b p_n(x)\,dx =\sum_{i=0}^{n} f(x_i)\int_a^b L_{i,n}(x)\,dx.\]作代换 $x=a+sh$，其中 $s\in[0,n]$，可得：
闭型 Newton-Cotes 公式
\[I_n(f) =h\sum_{i=0}^{n}\alpha_{i,n}f(x_i), \qquad \alpha_{i,n} =\int_0^n \prod_{\substack{j=0\\ j\ne i}}^{n} \frac{s-j}{i-j}\,ds. \tag{2.1}\]数 $\alpha_{0,n},\ldots,\alpha_{n,n}$ 称为权重。它们与 $f$ 以及区间 $[a,b]$ 无关，因此可以预先制表。并且总有
\[\sum_{i=0}^{n}h\alpha_{i,n}=b-a.\]定义 2.1.1
积分公式
如果它至少能精确积分所有次数不超过 $n$ 的多项式, 称为具有 $n$ 阶精确度。
闭型 Newton-Cotes 公式 $I_n(f)$ 按构造具有 $n$ 阶精确度。
重要的是估计误差
\[E_n(f):=I(f)-I_n(f).\]由推论 1.1.4 可知
\[|f(x)-p_n(x)| \le \frac{|f^{(n+1)}(\xi)|}{(n+1)!}(b-a)^{n+1},\]其中 $\xi\in[a,b]$。因此，对某个（可能不同的）$\xi\in[a,b]$ 有
\[\left| \int_a^b f(x)\,dx-\int_a^b p_n(x)\,dx \right| \le \int_a^b |f(x)-p_n(x)|\,dx \le \frac{|f^{(n+1)}(\xi)|}{(n+1)!}(b-a)^{n+2}.\]通过 Taylor 展开可以得到更精细的余项估计。结果如下表所示。
| $n$ | $\alpha_{i,n}$ | 最大误差 $E_n(f)$ | 名称 |
|---|---|---|---|
| 1 | $\frac12,\frac12$ | $-\frac{f^{(2)}(\xi)}{12}h^3$ | 梯形规则 |
| 2 | $\frac13,\frac43,\frac13$ | $-\frac{f^{(4)}(\xi)}{90}h^5$ | Simpson 规则 |
| 3 | $\frac38,\frac98,\frac98,\frac38$ | $-\frac{3f^{(4)}(\xi)}{80}h^5$ | 3/8 规则 |
| 4 | $\frac{14}{45},\frac{64}{45},\frac{24}{45},\frac{64}{45},\frac{14}{45}$ | $-\frac{8f^{(6)}(\xi)}{945}h^7$ | Milne 规则 |
需要注意的是，当 $n\ge 7$ 时，Newton-Cotes 权重会出现负值。负权重会带来相消，使公式在数值上越来越不稳定，所以实际计算中通常不会单纯通过增大 $n$ 来提高精度。
示意图 2.1： 闭型 Newton-Cotes 会使用端点。梯形规则用一条直线连接端点，Simpson 规则则用通过三个节点的二次曲线近似函数。
闭型公式把端点 $a,b$ 也作为节点。二次插值通常能更好贴合曲率，但高阶 Newton-Cotes 不一定更稳定。
例：考虑 $f(x)=\log_2(x)$ 在区间 $[a,b]=[2,4]$ 上的积分。精确积分为
\[I(f) =\int_2^4 f(x)\,dx =\int_2^4 \log_2(x)\,dx =\left.\frac{1}{\ln(2)}(x\ln(x)-x)\right|_2^4\] \[=\frac{4\ln(4)-4-2\ln(2)+2}{\ln(2)} \approx 3.11461.\]梯形规则给出：$h=\frac{2}{1}=2$，因此
\[I_1(f) =2\left(\frac12 f(2)+\frac12 f(4)\right) =2\left(\frac12\log_2(2)+\frac12\log_2(4)\right) =2\left(\frac12+1\right)=3.\]对 Simpson 规则，有 $h=\frac{2}{2}=1$，并且
\[I_2(f) =1\left(\frac13 f(2)+\frac43 f(3)+\frac13 f(4)\right) \approx 3.11328.\]使用 $\xi\in[2,4]$ 时，梯形规则的误差估计为
\[E_1(f) =-\frac{f''(\xi)h^3}{12} =\frac{1}{\xi^2\ln(2)}\frac{8}{12} \le \frac{1}{2^2\ln(2)}\frac{2}{3} \approx 0.24045.\]这与实际误差 $0.11461$ 可以接受地吻合。Simpson 规则的误差为
\[E_2(f) =-\frac{f^{(4)}(\xi)h^5}{90} =\frac{6}{\xi^4\ln(2)}\frac{1}{90} \le \frac{6}{2^4\ln(2)}\frac{1}{90} \approx 0.00601,\]这与实际误差 $0.00133$ 吻合得很好。
2.1.2 开型 Newton-Cotes 求积
这里对 $n\in\mathbb{N}\cup{0}$，选择位于开区间 $(a,b)$ 内的等距支撑点
\[x_i=a+ih,\qquad i=1,\ldots,n+1,\qquad h=\frac{b-a}{n+2}.\]用完全类似的方式处理，可以得到插值型积分公式：
开型 Newton-Cotes 公式
\[\tilde I_n(f) =h\sum_{i=1}^{n+1}\tilde\alpha_{i,n}f(x_i), \qquad \tilde\alpha_{i,n} = \int_0^{n+2} \prod_{\substack{j=1\\ j\ne i}}^{n+1} \frac{s-j}{i-j}\,ds.\]求积误差也有与闭型情形类似的公式：
| $n$ | $\tilde\alpha_{i,n}$ | 最大误差 $\tilde E_n(f)$ | 名称 |
|---|---|---|---|
| 0 | $2$ | $\frac{f^{(2)}(\xi)}{3}h^3$ | 矩形规则 |
| 1 | $\frac32,\frac32$ | $\frac{3f^{(2)}(\xi)}{4}h^3$ | |
| 2 | $\frac83,-\frac43,\frac83$ | $\frac{28f^{(4)}(\xi)}{90}h^5$ |
示意图 2.2： 开型 Newton-Cotes 不使用端点。矩形规则只取区间内部的一个点，常见做法是用中点函数值乘以区间长度。
当函数在端点不可取值、端点值不可靠，或者算法本身要求避开端点时，开型公式会更自然。
例：再次考虑 $f(x)=\log_2(x)$ 在 $[2,4]$ 上的积分，其精确值为 $I(f)\approx 3.11461$。矩形规则给出：$h=\frac22=1$，于是
\[\tilde I_0(f)=2f(3)\approx 3.16992.\]误差估计为
\[\tilde E_0(f) =\frac{f^{(2)}(\xi)h^3}{3} =-\frac{1}{\xi^2\ln(2)}\frac13 \ge -\frac{1}{2^2\ln(2)}\frac13 \approx -0.12022,\]而精确误差为 $-0.05529$。
2.2 复化 Newton-Cotes 公式
Newton-Cotes 公式只在积分区间较小、节点数不过大时比较可靠。因此，更常用的做法是把区间 $[a,b]$ 分解成多个小区间，并在每个小区间上分别使用 Newton-Cotes 公式，最后把这些子积分相加。
把区间 $[a,b]$ 分解为 $m$ 个长度为
\[H=\frac{b-a}{m}\]的子区间，在这些子区间上分别应用 $n$ 次 Newton-Cotes 公式，并把结果相加。令
\[N=m\cdot n,\qquad H=\frac{b-a}{m},\qquad h=\frac{H}{n}=\frac{b-a}{N},\] \[x_i=a+ih,\qquad i=0,\ldots,N,\] \[y_j=a+jH,\qquad j=0,\ldots,m.\]由于
\[I(f)=\sum_{j=0}^{m-1}\int_{y_j}^{y_{j+1}}f(x)\,dx,\]得到：复化闭型 Newton-Cotes 公式
\[S_N^{(n)}(f) =h\sum_{j=0}^{m-1}\sum_{i=0}^{n} \alpha_{i,n}f(x_{jn+i}).\]权重 $\alpha_{i,n}$ 仍由 (2.1) 给出。求积误差
\[R_N^{(n)}(f)=I(f)-S_N^{(n)}(f)\]可以通过累加各个子区间上的误差得到。
定理 2.2.1：设 $f\in C^{n+2}([a,b])$。则存在中间点 $\xi\in(a,b)$，使得
\[R_N^{(n)}(f)= \begin{cases} C(n)f^{(n+2)}(\xi)(b-a)h^{n+2}, & n \text{ 为偶数},\\ C(n)f^{(n+1)}(\xi)(b-a)h^{n+1}, & n \text{ 为奇数}. \end{cases}\]这里 $C(n)$ 是只依赖于 $n$ 的常数。
下面列出最常用的几个复化公式及其求积误差。
复化梯形规则
闭型，$n=1$，$h=\frac{b-a}{m}$：
\[S_N^{(1)}(f) =\frac{h}{2} \sum_{j=0}^{m-1} \left(f(x_j)+f(x_{j+1})\right), \qquad x_j=a+jh.\]误差：
\[R_N^{(1)}(f) =-\frac{f''(\xi)}{12}(b-a)h^2.\]复化 Simpson 规则
闭型，$n=2$，$h=\frac{b-a}{2m}$：
\[S_N^{(2)}(f) =\frac{h}{3} \sum_{j=0}^{m-1} \left(f(x_{2j})+4f(x_{2j+1})+f(x_{2j+2})\right), \qquad x_j=a+jh.\]误差：
\[R_N^{(2)}(f) =-\frac{f^{(4)}(\xi)}{180}(b-a)h^4.\]复化矩形规则
开型，$n=0$，$2m=N$，$h=\frac{b-a}{N}$：
\[\tilde S_N^{(0)}(f) =2h\sum_{j=1}^{m} f(x_{2j-1}), \qquad x_j=a+jh.\]误差：
\[\tilde R_N^{(0)}(f) =\frac{f''(\xi)}{6}(b-a)h^2.\]示意图 2.3： 复化公式通过缩小步长 $h$ 来降低误差。梯形规则的误差主阶通常随 $h^2$ 下降，Simpson 规则则随 $h^4$ 下降。
复化思想的重点不是把单个区间上的插值次数无限提高，而是把区间切细，让每一段上的低阶公式更可靠。
例：再次考虑 $f(x)=\log_2(x)$ 在 $[2,4]$ 上的积分，精确值为 $I(f)\approx 3.11461$，并取 $m=2$。
-
对复化梯形规则，有 $n=1$、$N=m=2$、$h=\frac22=1$，得到
\[S_N^{(1)}(f) =\frac12\left(f(2)+f(3)+f(3)+f(4)\right) \approx \frac12(1+1.5849+1.5849+2) =3.08496.\]误差估计给出
\[R_N^{(1)}(f) =-\frac{f''(\xi)}{12}(b-a)h^2 = \frac{1}{\xi^2\ln(2)}\frac{2}{12} \le \frac{1}{2^2\ln(2)}\frac{2}{12} \approx 0.06011.\]这明显减小了简单梯形规则的误差；真实误差为 $0.02965$。
-
对复化 Simpson 规则，有 $n=2$、$N=4$、$h=\frac24=\frac12$，并且
\[S_N^{(2)}(f) = \frac{1}{2\cdot 3} \left( (f(2)+4f(2.5)+f(3)) +(f(3)+4f(3.5)+f(4)) \right) \approx 3.11450.\]误差估计为 $R_N^{(2)}(f)\le 0.00037$，真实误差为 $0.00011$。
-
对复化矩形规则，得到 $h=\frac24=\frac12$，并且
\[\tilde S_N^{(0)}(f) = \frac22\bigl(f(2.5)+f(3.5)\bigr) \approx 3.12928.\]误差估计为
\[\tilde R_N^{(0)}(f) = \frac{f''(\xi)}{6}(b-a)h^2 = -\frac{1}{\xi^2\ln(2)}\frac16\frac24 \ge -\frac{1}{2^2\ln(2)}\frac16\frac12 \approx -0.03005,\]真实误差为 $-0.01467$。
来源、版权与使用说明
本文主要参考 TU Darmstadt 信息学专业公开仓库中的数值分析基础课讲义： mathe3-script-2011-SoSe.pdf 原仓库包含 The Unlicense 授权说明。本文作为个人学习、翻译与知识整理用途发布，文中的中文表述、补充解释和图表重制不代表原作者或官方立场。
本文中的个人整理、中文表述、补充解释以及我重新制作的图表，可在注明作者与原文链接的前提下，用于非商业学习、交流和引用。由于本文部分内容基于 TU Darmstadt 公开讲义的翻译与整理，原始讲义及其中可能包含的材料仍应以其原作者、原仓库及相关授权说明为准。若需进行商业使用、系统转载、出版，或大规模改编，建议同时确认原始材料的授权状态。
如文中存在翻译、公式、术语或理解上的疏漏，或相关权利方认为内容使用不当，欢迎联系我指出，我会及时处理或删除。
AI feedback
AnonymousLoading AI feedback…