第六章 数值积分

在科学计算中，常常需要计算定积分： $$I = \int_a^b f(x) dx$$

但存在以下困难： 1. $f(x)$ 的原函数无法用初等函数表示（如 $e^{-x^2}$，$\frac{\sin x}{x}$ 等） 2. $f(x)$ 的表达式未知，只知道离散数据点 3. 原函数表达式过于复杂

因此需要数值积分方法，通过被积函数的有限个函数值来近似计算积分值。

$$\int_a^b f(x) dx \approx \sum_{k=0}^{n} A_k f(x_k)$$

其中：

1. $x_k$：求积节点
2. $A_k$：求积系数（仅与节点有关，与被积函数无关）

这样的公式称为机械求积公式。

定义 6.1（代数精度） 若求积公式对所有次数不超过 $m$ 的多项式精确成立，但对 $m+1$ 次多项式不精确成立，则称该公式具有 $m$ 次代数精度。

用Lagrange插值多项式 $L_n(x)$ 近似 $f(x)$： $$\int_a^b f(x) dx \approx \int_a^b L_n(x) dx = \sum_{k=0}^{n} f(x_k) \int_a^b l_k(x) dx = \sum_{k=0}^{n} A_k f(x_k)$$

其中 $A_k = \int_a^b l_k(x) dx$。

定理 6.1 $n+1$ 个节点的插值型求积公式至少具有 $n$ 次代数精度。

取等距节点 $x_k = a + kh$，$h = \frac{b-a}{n}$，$k = 0, 1, \ldots, n$。

令 $x = a + th$，则： $$A_k = \frac{b-a}{n} \int_0^n \prod_{j=0, j \neq k}^{n} \frac{t-j}{k-j} dt = (b-a) C_k^{(n)}$$

其中 $C_k^{(n)}$ 称为 Cotes系数，满足 $\sum_{k=0}^{n} C_k^{(n)} = 1$。

Newton-Cotes公式： $$\int_a^b f(x) dx \approx (b-a) \sum_{k=0}^{n} C_k^{(n)} f(x_k)$$

梯形公式（$n=1$）： $$\int_a^b f(x) dx \approx \frac{b-a}{2}[f(a) + f(b)]$$

Cotes系数：$C_0^{(1)} = C_1^{(1)} = \frac{1}{2}$

代数精度：1次

Simpson公式（$n=2$，抛物线公式）： $$\int_a^b f(x) dx \approx \frac{b-a}{6}\left[f(a) + 4f\left(\frac{a+b}{2}\right) + f(b)\right]$$

Cotes系数：$C_0^{(2)} = \frac{1}{6}$，$C_1^{(2)} = \frac{4}{6}$，$C_2^{(2)} = \frac{1}{6}$

代数精度：3次（由于对称性，高于预期的2次）

Simpson 3/8公式（$n=3$）： $$\int_a^b f(x) dx \approx \frac{b-a}{8}\left[f(a) + 3f\left(\frac{2a+b}{3}\right) + 3f\left(\frac{a+2b}{3}\right) + f(b)\right]$$

代数精度：3次

Cotes公式（$n=4$）： $$\int_a^b f(x) dx \approx \frac{b-a}{90}\left[7f(a) + 32f(x_1) + 12f(x_2) + 32f(x_3) + 7f(b)\right]$$

代数精度：5次

Cotes系数表：

当 $n \geq 8$ 时，Cotes系数出现负值，公式数值不稳定。因此高阶Newton-Cotes公式很少使用。

定理 6.2 设 $f(x) \in C^2[a, b]$，则梯形公式的截断误差为： $$R_T = \int_a^b f(x) dx - \frac{b-a}{2}[f(a)+f(b)] = -\frac{(b-a)^3}{12}f(\xi), \quad \xi \in (a, b)$$ ==== 6.4.2 Simpson公式的误差 ==== 定理 6.3 设 $f(x) \in C^4[a, b]$，则Simpson公式的截断误差为： $$R_S = \int_a^b f(x) dx - \frac{b-a}{6}\left[f(a) + 4f\left(\frac{a+b}{2}\right) + f(b)\right] = -\frac{(b-a)^5}{2880}f^{(4)}(\xi)$$ ===== 6.5 复化求积公式 ===== ==== 6.5.1 复化梯形公式 ==== 将 $[a, b]$ 分成 $n$ 等份，$h = \frac{b-a}{n}$，$x_k = a + kh$，在每个小区间上用梯形公式： $$T_n = \sum_{k=0}^{n-1} \frac{h}{2}[f(x_k) + f(x_{k+1})] = \frac{h}{2}\left[f(a) + 2\sum_{k=1}^{n-1}f(x_k) + f(b)\right]$$ 误差： $$R_n = -\frac{b-a}{12}h^2 f(\xi) = O(h^2)$$

将 $[a, b]$ 分成 $2n$ 等份（$n$ 个小区间，每个小区间用Simpson公式），$h = \frac{b-a}{2n}$：

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

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

误差： $$R_n = -\frac{b-a}{180}\left(\frac{h}{2}\right)^4 f^{(4)}(\xi) = O(h^4)$$

$$C_n = \frac{h}{90}\left[7f(a) + 32\sum f_{\text{奇}} + 12\sum f_{\text{特殊}} + 32\sum f_{\text{奇}} + 7f(b)\right]$$

误差：$O(h^6)$

设 $I(h)$ 是步长为 $h$ 时某量的近似，有展开式： $$I(h) = I + c_1 h^p + c_2 h^{2p} + c_3 h^{3p} + \cdots$$

则： $$I\left(\frac{h}{2}\right) = I + c_1 \frac{h^p}{2^p} + c_2 \frac{h^{2p}}{4^p} + \cdots$$

两式组合消去 $h^p$ 项： $$I_1(h) = \frac{2^p I(h/2) - I(h)}{2^p - 1} = I + O(h^{2p})$$

对梯形公式（$p=2$）： $$T_0^{(k)} = T_{2^k}$$（$k$ 层二分后的梯形值）

递推公式： $$T_m^{(k)} = \frac{4^m T_{m-1}^{(k+1)} - T_{m-1}^{(k)}}{4^m - 1}$$

Romberg表：

其中：

1. $T_0^{(k)}$：复化梯形值
2. $T_1^{(k)}$：复化Simpson值
3. $T_2^{(k)}$：复化Cotes值
4. $T_3^{(k)}$：Romberg值

收敛速度：$T_m^{(k)} = I + O(h^{2(m+1)})$

$n+1$ 个节点的求积公式最高可达到 $2n+1$ 次代数精度。达到此精度的公式称为Gauss求积公式，相应的节点称为Gauss点。

定理 6.4 节点 $x_0, x_1, \ldots, x_n$ 是Gauss点的充分必要条件是：多项式 $\omega_{n+1}(x) = \prod_{k=0}^{n}(x-x_k)$ 与任意次数不超过 $n$ 的多项式 $P(x)$ 正交，即： $$\int_a^b \rho(x) \omega_{n+1}(x) P(x) dx = 0$$

推论：Gauss点是 $[a, b]$ 上关于权函数 $\rho(x)$ 的 $n+1$ 次正交多项式的零点。

区间：$[-1, 1]$，权函数：$\rho(x) = 1$

Gauss点为Legendre多项式 $P_{n+1}(x)$ 的零点。

求积系数： $$A_k = \frac{2}{(1-x_k^2)[P'_{n+1}(x_k)]^2}$$

Gauss-Legendre求积节点和系数（部分）：

$n=1$（2点）：

1. $x_0 = -\frac{1}{\sqrt{3}} \approx -0.5774$，$A_0 = 1$
2. $x_1 = \frac{1}{\sqrt{3}} \approx 0.5774$，$A_1 = 1$

$n=2$（3点）：

1. $x_0 = -\sqrt{\frac{3}{5}} \approx -0.7746$，$A_0 = \frac{5}{9}$
2. $x_1 = 0$，$A_1 = \frac{8}{9}$
3. $x_2 = \sqrt{\frac{3}{5}} \approx 0.7746$，$A_2 = \frac{5}{9}$

$n=3$（4点）：

1. $x_0 = -0.8611$，$A_0 = 0.3479$
2. $x_1 = -0.3400$，$A_1 = 0.6521$
3. $x_2 = 0.3400$，$A_2 = 0.6521$
4. $x_3 = 0.8611$，$A_3 = 0.3479$

Gauss-Chebyshev求积： 区间 $[-1, 1]$，权函数 $\rho(x) = \frac{1}{\sqrt{1-x^2}}$

节点：$x_k = \cos\left(\frac{2k+1}{2n+2}\pi\right)$，$A_k = \frac{\pi}{n+1}$

Gauss-Laguerre求积： 区间 $[0, +\infty)$，权函数 $\rho(x) = e^{-x}$

Gauss-Hermite求积： 区间 $(-\infty, +\infty)$，权函数 $\rho(x) = e^{-x^2}$

例题 6.1 用梯形公式和Simpson公式计算 $I = \int_0^1 e^x dx$。

解：

梯形公式： $$T = \frac{1}{2}(e^0 + e^1) = \frac{1}{2}(1 + 2.71828) = 1.85914$$

（精确值 $e - 1 \approx 1.71828$，误差约 $0.14$）

Simpson公式： $$S = \frac{1}{6}\left[e^0 + 4e^{0.5} + e^1\right] = \frac{1}{6}[1 + 4 \times 1.64872 + 2.71828] = 1.71886$$

（误差约 $0.0006$）

例题 6.2 用Romberg算法计算 $\int_0^1 \frac{4}{1+x^2} dx$（精确值为 $\pi$）。

解：

$T_0^{(0)} = \frac{1}{2}[f(0) + f(1)] = \frac{1}{2}[4 + 2] = 3$

$T_0^{(1)} = \frac{1}{4}[f(0) + 2f(0.5) + f(1)] = \frac{1}{4}[4 + 2 \times 3.2 + 2] = 3.1$

$T_1^{(0)} = \frac{4T_0^{(1)} - T_0^{(0)}}{3} = \frac{12.4 - 3}{3} = 3.13333$

继续二分计算Romberg表，可快速收敛到 $\pi \approx 3.14159$。

基础练习

1. 确定下列求积公式的代数精度：

 $$\int_{-1}^1 f(x)dx \approx f\left(-\frac{1}{\sqrt{3}}\right) + f\left(\frac{1}{\sqrt{3}}\right)$$

2. 推导梯形公式的误差公式。

3. 用复化Simpson公式（$n=4$）计算 $\int_0^1 e^{-x^2} dx$。

进阶练习

4. 证明Simpson公式对三次多项式精确成立。

5. 设 $f \in C^4[a, b]$，证明复化Simpson公式的误差：

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

6. 推导Romberg递推公式 $T_m^{(k)} = \frac{4^m T_{m-1}^{(k+1)} - T_{m-1}^{(k)}}{4^m - 1}$。

编程实践

7. 编写复化Simpson公式的自适应版本，根据精度要求自动确定 $n$。

8. 实现Romberg算法，计算 $\int_0^1 \frac{\sin x}{x} dx$。

1. 李庆扬等. 数值分析（第5版）. 清华大学出版社, 2008. 2. Davis P J, Rabinowitz P. Methods of Numerical Integration. Academic Press, 1984.