====== 第七章 数值微分 ======

===== 7.1 引言 =====

数值微分是用离散方法近似计算函数导数或微分的方法。给定函数 $f(x)$ 在离散点上的值，估计其导数 $f'(x)$、二阶导数 $f''(x)$ 等。

数值微分的重要性：
1. 函数以表格形式给出时
2. 微分方程数值解中需要近似导数
3. 优化算法中需要计算梯度

===== 7.2 差分公式 =====

==== 7.2.1 基于Taylor展开的差分公式 ====

由Taylor展开：
$$f(x+h) = f(x) + hf'(x) + \frac{h^2}{2}f''(x) + \frac{h^3}{6}f'''(x) + O(h^4)$$

$$f(x-h) = f(x) - hf'(x) + \frac{h^2}{2}f''(x) - \frac{h^3}{6}f'''(x) + O(h^4)$$

==== 7.2.2 一阶差分公式 ====

**向前差分**：
$$f'(x) \approx \frac{f(x+h) - f(x)}{h}$$

误差：$R = -\frac{h}{2}f''(\xi) = O(h)$，一阶精度

**向后差分**：
$$f'(x) \approx \frac{f(x) - f(x-h)}{h}$$

误差：$R = \frac{h}{2}f''(\xi) = O(h)$，一阶精度

**中心差分**：
$$f'(x) \approx \frac{f(x+h) - f(x-h)}{2h}$$

误差：$R = -\frac{h^2}{6}f'''(\xi) = O(h^2)$，二阶精度

==== 7.2.3 二阶差分公式 ====

**二阶中心差分**：
$$f''(x) \approx \frac{f(x+h) - 2f(x) + f(x-h)}{h^2}$$

误差：$R = -\frac{h^2}{12}f^{(4)}(\xi) = O(h^2)$

推导：将 $f(x+h)$ 和 $f(x-h)$ 的Taylor展开相加：
$$f(x+h) + f(x-h) = 2f(x) + h^2f''(x) + \frac{h^4}{12}f^{(4)}(\xi)$$

整理即得。

==== 7.2.4 高阶差分公式 ====

利用更多点可得到更高精度的公式。

**三点向前差分（二阶精度）**：
$$f'(x) \approx \frac{-3f(x) + 4f(x+h) - f(x+2h)}{2h}$$

**三点向后差分（二阶精度）**：
$$f'(x) \approx \frac{3f(x) - 4f(x-h) + f(x-2h)}{2h}$$

**五点中心差分（四阶精度）**：
$$f'(x) \approx \frac{f(x-2h) - 8f(x-h) + 8f(x+h) - f(x+2h)}{12h}$$

误差：$O(h^4)$

===== 7.3 插值型数值微分 =====

==== 7.3.1 两点公式 ====

过 $(x_0, f(x_0))$ 和 $(x_1, f(x_1))$ 的线性插值：
$$L_1(x) = f(x_0) + f[x_0, x_1](x - x_0)$$

则：
$$f'(x) \approx L'_1(x) = f[x_0, x_1] = \frac{f(x_1) - f(x_0)}{x_1 - x_0}$$

**端点公式**：
$$f'(x_0) = \frac{f(x_1) - f(x_0)}{h} - \frac{h}{2}f''(\xi)$$（向前）
$$f'(x_1) = \frac{f(x_1) - f(x_0)}{h} + \frac{h}{2}f''(\xi)$$（向后）

==== 7.3.2 三点公式 ====

过 $(x_0, f(x_0))$、$(x_1, f(x_1))$、$(x_2, f(x_2))$ 的抛物线插值，$x_k = x_0 + kh$：

**中点公式**：
$$f'(x_1) = \frac{f(x_2) - f(x_0)}{2h} - \frac{h^2}{6}f'''(\xi)$$

**端点公式**：
$$f'(x_0) = \frac{-3f(x_0) + 4f(x_1) - f(x_2)}{2h} + \frac{h^2}{3}f'''(\xi)$$
$$f'(x_2) = \frac{f(x_0) - 4f(x_1) + 3f(x_2)}{2h} + \frac{h^2}{3}f'''(\xi)$$

**二阶导数公式**：
$$f''(x_1) = \frac{f(x_0) - 2f(x_1) + f(x_2)}{h^2} - \frac{h^2}{12}f^{(4)}(\xi)$$

===== 7.4 外推法 =====

==== 7.4.1 Richardson外推原理 ====

设 $I(h)$ 是 $I$ 的近似，具有渐近展开：
$$I(h) = I + c_1 h^p + c_2 h^{2p} + c_3 h^{3p} + \cdots$$

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

组合得更高阶近似：
$$G_1(h) = \frac{2^p I(h/2) - I(h)}{2^p - 1} = I + O(h^{2p})$$

可继续外推：
$$G_k(h) = \frac{2^{kp} G_{k-1}(h/2) - G_{k-1}(h)}{2^{kp} - 1}$$

==== 7.4.2 应用于数值微分 ====

对于中心差分公式：
$$f'(x) = \frac{f(x+h) - f(x-h)}{2h} - \frac{h^2}{6}f'''(x) - \frac{h^4}{120}f^{(5)}(x) - \cdots$$

令 $G_0(h) = \frac{f(x+h) - f(x-h)}{2h}$，则：

**一次外推**：
$$G_1(h) = \frac{4G_0(h/2) - G_0(h)}{3} = \frac{-f(x+h) + 8f(x+h/2) - 8f(x-h/2) + f(x-h)}{6h}$$

误差：$O(h^4)$

**二次外推**：
$$G_2(h) = \frac{16G_1(h/2) - G_1(h)}{15}$$

误差：$O(h^6)$

===== 7.5 差分公式的稳定性分析 =====

==== 7.5.1 舍入误差的影响 ====

设 $f(x \pm h)$ 的舍入误差为 $\varepsilon_{\pm}$，$|\varepsilon_{\pm}| \leq \varepsilon$。

对中心差分：
$$f'(x) \approx \frac{\tilde{f}(x+h) - \tilde{f}(x-h)}{2h} = \frac{f(x+h) - f(x-h)}{2h} + \frac{\varepsilon_+ - \varepsilon_-}{2h}$$

舍入误差上界：$\frac{\varepsilon}{h}$。

==== 7.5.2 最优步长 ====

总误差 = 截断误差 + 舍入误差

对中心差分：
$$\text{总误差} \approx \frac{h^2}{6}M_3 + \frac{\varepsilon}{h}$$

其中 $M_3 = \max |f'''(x)|$。

对 $h$ 求导并令为0，得**最优步长**：
$$h_{opt} = \sqrt[3]{\frac{3\varepsilon}{M_3}}$$

若 $\varepsilon \approx 10^{-16}$（双精度），$M_3 \approx 1$，则 $h_{opt} \approx 10^{-5}$。

===== 7.6 偏导数的数值计算 =====

==== 7.6.1 一阶偏导数 ====

$$\frac{\partial f}{\partial x}(x, y) \approx \frac{f(x+h, y) - f(x-h, y)}{2h}$$

$$\frac{\partial f}{\partial y}(x, y) \approx \frac{f(x, y+k) - f(x, y-k)}{2k}$$

==== 7.6.2 二阶偏导数 ====

$$\frac{\partial^2 f}{\partial x^2}(x, y) \approx \frac{f(x+h, y) - 2f(x, y) + f(x-h, y)}{h^2}$$

$$\frac{\partial^2 f}{\partial y^2}(x, y) \approx \frac{f(x, y+k) - 2f(x, y) + f(x, y-k)}{k^2}$$

**混合偏导数**：
$$\frac{\partial^2 f}{\partial x \partial y} \approx \frac{f(x+h, y+k) - f(x+h, y-k) - f(x-h, y+k) + f(x-h, y-k)}{4hk}$$

===== 7.7 典型例题 =====

**例题 7.1** 已知 $f(x) = e^x$ 在以下点的值：
| $x$ | 1.8 | 1.9 | 2.0 | 2.1 | 2.2 |
|-----|-----|-----|-----|-----|-----|
| $f(x)$ | 6.0496 | 6.6859 | 7.3891 | 8.1662 | 9.0250 |

求 $f'(2.0)$ 的近似值。

**解**：

用中心差分（$h=0.1$）：
$$f'(2.0) \approx \frac{f(2.1) - f(1.9)}{0.2} = \frac{8.1662 - 6.6859}{0.2} = 7.4015$$

用五点公式：
$$f'(2.0) \approx \frac{f(1.8) - 8f(1.9) + 8f(2.1) - f(2.2)}{1.2}$$
$$= \frac{6.0496 - 8 \times 6.6859 + 8 \times 8.1662 - 9.0250}{1.2} = 7.3891$$

（精确值 $f'(2.0) = e^{2.0} = 7.3891$，五点公式更精确）

**例题 7.2** 用Richardson外推求 $f(x) = \sin x$ 在 $x = \frac{\pi}{3}$ 处的导数。

**解**：$f'(x) = \cos x$，$f'(\frac{\pi}{3}) = 0.5$。

$G_0(h) = \frac{\sin(\frac{\pi}{3}+h) - \sin(\frac{\pi}{3}-h)}{2h}$

$G_0(\frac{\pi}{6}) = \frac{\sin(\frac{\pi}{2}) - \sin(\frac{\pi}{6})}{\frac{\pi}{3}} = \frac{1 - 0.5}{1.047} = 0.4775$

$G_0(\frac{\pi}{12}) = \frac{\sin(\frac{5\pi}{12}) - \sin(\frac{\pi}{4})}{\frac{\pi}{6}} = 0.4970$

一次外推：
$$G_1(\frac{\pi}{6}) = \frac{4 \times 0.4970 - 0.4775}{3} = 0.5032$$

已接近精确值 $0.5$。

===== 7.8 习题 =====

**基础练习**

1. 推导二阶中心差分公式并给出截断误差。

2. 已知 $f(0) = 1$，$f(0.1) = 1.1052$，$f(0.2) = 1.2214$，用三点公式求 $f'(0.1)$。

3. 证明五点中心差分公式：
   $$f'(x) = \frac{-f(x+2h) + 8f(x+h) - 8f(x-h) + f(x-2h)}{12h} + O(h^4)$$

**进阶练习**

4. 分析向前差分公式的舍入误差，并求最优步长。

5. 设 $f \in C^5[a, b]$，推导用Richardson外推求 $f'(x)$ 的三阶近似公式。

6. 证明：若 $f$ 是 $n$ 次多项式，则 $n$ 阶差商 $f[x_0, x_1, \ldots, x_n]$ 为常数（与节点位置无关）。

**编程实践**

7. 编写程序比较不同差分公式的精度随步长的变化。

8. 实现Richardson外推算法，计算函数在某点的各阶导数。

===== 参考文献 =====

1. 李庆扬等. 数值分析（第5版）. 清华大学出版社, 2008.
2. Fornberg B. Generation of Finite Difference Formulas on Arbitrarily Spaced Grids. Mathematics of Computation, 1988.