张叶安的小站

本页面只读。您可以查看源文件，但不能更改它。如果您觉得这是系统错误，请联系管理员。
====== 第六章 数值积分 ======

===== 6.1 引言 =====

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

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

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

===== 6.2 数值积分的基本概念 =====

==== 6.2.1 求积公式的一般形式 ====

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

其中：
  - $x_k$：求积节点
  - $A_k$：求积系数（仅与节点有关，与被积函数无关）

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

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

===== 6.3 Newton-Cotes公式 =====

==== 6.3.1 插值型求积公式 ====

用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$ 次代数精度。

==== 6.3.2 Newton-Cotes公式 ====

取等距节点 $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)$$

==== 6.3.3 常用Newton-Cotes公式 ====

**梯形公式**（$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$ | $C_k^{(n)}$ |
|-----|-------------|
| 1 | $\frac{1}{2}, \frac{1}{2}$ |
| 2 | $\frac{1}{6}, \frac{4}{6}, \frac{1}{6}$ |
| 3 | $\frac{1}{8}, \frac{3}{8}, \frac{3}{8}, \frac{1}{8}$ |
| 4 | $\frac{7}{90}, \frac{32}{90}, \frac{12}{90}, \frac{32}{90}, \frac{7}{90}$ |
| 5 | $\frac{19}{288}, \frac{75}{288}, \frac{50}{288}, \frac{50}{288}, \frac{75}{288}, \frac{19}{288}$ |
| 6 | $\frac{41}{840}, \frac{216}{840}, \frac{27}{840}, \frac{272}{840}, \frac{27}{840}, \frac{216}{840}, \frac{41}{840}$ |

==== 6.3.4 Newton-Cotes公式的稳定性 ====

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

===== 6.4 求积公式的误差估计 =====

==== 6.4.1 梯形公式的误差 ====

**定理 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)$$

==== 6.5.2 复化Simpson公式 ====

将 $[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)$$

==== 6.5.3 复化Cotes公式 ====

$$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)$

===== 6.6 Romberg积分 =====

==== 6.6.1 Richardson外推 ====

设 $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})$$

==== 6.6.2 Romberg算法 ====

对梯形公式（$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表**：
| $k$ | $T_0^{(k)}$ | $T_1^{(k)}$ | $T_2^{(k)}$ | $T_3^{(k)}$ | ... |
|-----|-------------|-------------|-------------|-------------|-----|
| 0 | $T_0^{(0)}$ | | | | |
| 1 | $T_0^{(1)}$ | $T_1^{(0)}$ | | | |
| 2 | $T_0^{(2)}$ | $T_1^{(1)}$ | $T_2^{(0)}$ | | |
| 3 | $T_0^{(3)}$ | $T_1^{(2)}$ | $T_2^{(1)}$ | $T_3^{(0)}$ | |
| ... | ... | ... | ... | ... | ... |

其中：
  - $T_0^{(k)}$：复化梯形值
  - $T_1^{(k)}$：复化Simpson值
  - $T_2^{(k)}$：复化Cotes值
  - $T_3^{(k)}$：Romberg值

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

===== 6.7 Gauss求积公式 =====

==== 6.7.1 代数精度的最高限度 ====

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

==== 6.7.2 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$ 次正交多项式的零点。

==== 6.7.3 Gauss-Legendre求积 ====

区间：$[-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点）：
  - $x_0 = -\frac{1}{\sqrt{3}} \approx -0.5774$，$A_0 = 1$
  - $x_1 = \frac{1}{\sqrt{3}} \approx 0.5774$，$A_1 = 1$

$n=2$（3点）：
  - $x_0 = -\sqrt{\frac{3}{5}} \approx -0.7746$，$A_0 = \frac{5}{9}$
  - $x_1 = 0$，$A_1 = \frac{8}{9}$
  - $x_2 = \sqrt{\frac{3}{5}} \approx 0.7746$，$A_2 = \frac{5}{9}$

$n=3$（4点）：
  - $x_0 = -0.8611$，$A_0 = 0.3479$
  - $x_1 = -0.3400$，$A_1 = 0.6521$
  - $x_2 = 0.3400$，$A_2 = 0.6521$
  - $x_3 = 0.8611$，$A_3 = 0.3479$

==== 6.7.4 其他Gauss型求积 ====

**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.8 典型例题 =====

**例题 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$。

===== 6.9 习题 =====

**基础练习**

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.