第七章 数值微分

数值微分是用离散方法近似计算函数导数或微分的方法。给定函数 $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)$，二阶精度

二阶中心差分： $$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)$$

整理即得。

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

三点向前差分（二阶精度）： $$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)$

过 $(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)$$（向后）

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

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

对于中心差分公式： $$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}$。

$$\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}$$

$$\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.1 已知 $f(x) = e^x$ 在以下点的值：

求 $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$。

基础练习

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.