数值微分三点公式的误差分析与步长选择
问题描述
数值微分的核心问题在于，当函数的解析导数不可用或难以计算时，如何利用函数在一些离散点上的函数值来近似计算其导数值。三点公式是数值微分中最常用且基础的方法之一，它利用一个点及其左右相邻两个点的函数值来构造导数的近似。然而，这个近似会引入误差，且误差与所选取的节点步长 \(h\) 密切相关。本题目旨在详细分析三点公式的误差来源，推导其截断误差表达式，并探讨如何选择最优的步长 \(h\) 以平衡截断误差和舍入误差，从而获得最精确的导数值近似。
背景与动机
假设我们有一个光滑函数 \(f(x)\)，但我们没有其导函数 \(f'(x)\) 的解析表达式，只知道它在某些离散点 \(x_0, x_1, x_2, ...\) 上的函数值。数值微分的目标就是通过这些离散的函数值来估计 \(f'(x)\) 在某个点（比如中点）的值。三点公式是构建这种近似的一个简单而有效的方法。但使用不恰当的步长可能导致结果极不准确，因此理解其误差构成并优化步长至关重要。
解题过程
步骤1：三点公式的推导
我们考虑三个等距节点：\(x_0 - h, x_0, x_0 + h\)，对应的函数值为 \(f(x_0 - h), f(x_0), f(x_0 + h)\)。
- 中心差分公式（近似 \(f'(x_0)\)）：
这是最常用的三点公式。我们构造通过这三点的二次插值多项式 \(P_2(x)\)，然后对 \(P_2(x)\) 在 \(x_0\) 处求导，作为 \(f'(x_0)\) 的近似。
利用 Lagrange 插值或 Taylor 展开，可以得到：
\[ f'(x_0) \approx \frac{f(x_0 + h) - f(x_0 - h)}{2h} \]
这个公式称为**中心差分公式**。
- 前向差分公式（近似 \(f'(x_0)\)）：
我们也可以利用点 \(x_0, x_0 + h, x_0 + 2h\) 来构造。但更常见的是，用点 \(x_0, x_0+h, x_0+2h\) 构造二次多项式，并求其在 \(x_0\) 处的导数，得到三点前向差分公式：
\[ f'(x_0) \approx \frac{-3f(x_0) + 4f(x_0 + h) - f(x_0 + 2h)}{2h} \]
- 后向差分公式（近似 \(f'(x_0)\)）：
同理，利用点 \(x_0, x_0-h, x_0-2h\) 可得三点后向差分公式：
\[ f'(x_0) \approx \frac{3f(x_0) - 4f(x_0 - h) + f(x_0 - 2h)}{2h} \]
本题目我们将以最常用的**中心差分公式**为主要分析对象。
步骤2：截断误差分析
截断误差来源于用有限项 Taylor 展开的多项式来近似无穷可微的函数。
- 对中心差分公式进行 Taylor 展开：
\[ f(x_0 + h) = f(x_0) + h f'(x_0) + \frac{h^2}{2!} f''(x_0) + \frac{h^3}{3!} f'''(x_0) + \frac{h^4}{4!} f^{(4)}(x_0) + O(h^5) \]
\[ f(x_0 - h) = f(x_0) - h f'(x_0) + \frac{h^2}{2!} f''(x_0) - \frac{h^3}{3!} f'''(x_0) + \frac{h^4}{4!} f^{(4)}(x_0) + O(h^5) \]
- 将两式相减：
\[ f(x_0 + h) - f(x_0 - h) = 2h f'(x_0) + \frac{2h^3}{3!} f'''(x_0) + O(h^5) \]
\[ f(x_0 + h) - f(x_0 - h) = 2h f'(x_0) + \frac{h^3}{3} f'''(x_0) + O(h^5) \]
- 推导中心差分公式及其误差：
将上式两边同除以 \(2h\)：
\[ \frac{f(x_0 + h) - f(x_0 - h)}{2h} = f'(x_0) + \frac{h^2}{6} f'''(x_0) + O(h^4) \]
因此，中心差分公式的截断误差 $ E_T $ 为：
\[ E_T = f'(x_0) - \frac{f(x_0 + h) - f(x_0 - h)}{2h} = -\frac{h^2}{6} f'''(x_0) + O(h^4) \]
**关键结论**：中心差分公式的截断误差主项与 $ h^2 $ 成正比。我们说这个公式具有**二阶精度**。当 $ h $ 减小时，截断误差以 $ h^2 $ 的速度减小。
步骤3：舍入误差分析
在实际计算中，函数值 \(f(x)\) 是带有舍入误差的。设机器精度为 \(\epsilon\)（例如，双精度浮点数 \(\epsilon \approx 2.22 \times 10^{-16}\)），则计算得到的函数值 \(\tilde{f}(x)\) 满足 \(|\tilde{f}(x) - f(x)| \lessapprox \epsilon |f(x)|\)。为简化分析，通常假设函数值的绝对误差界为 \(\delta\)，即 \(|\tilde{f}(x) - f(x)| \leq \delta\)，其中 \(\delta\) 与函数值和机器精度相关。
- 中心差分公式的舍入误差传播：
我们计算的是 \(\tilde{D} = \frac{\tilde{f}(x_0 + h) - \tilde{f}(x_0 - h)}{2h}\)。
每个函数值的误差最大为 \(\delta\)，则分子的绝对误差最大约为 \(2\delta\)。
因此，最终结果的舍入误差 \(E_R\) 满足：
\[ |E_R| = |\tilde{D} - D| \lessapprox \frac{2\delta}{2h} = \frac{\delta}{h} \]
其中 $ D $ 是使用精确函数值计算的中心差分结果。
- 关键观察：舍入误差与 \(1/h\) 成正比。当 \(h\) 非常小时，分母 \(h\) 很小，会放大函数值中的微小误差（包括测量误差和计算机的舍入误差），导致计算结果极不稳定。
步骤4：总误差与最优步长选择
总误差 \(E_{total}\) 是截断误差 \(E_T\) 和舍入误差 \(E_R\) 的绝对值之和的一个上界估计。
- 总误差模型：
\[ |E_{total}| \lessapprox |E_T| + |E_R| \approx \frac{M h^2}{6} + \frac{\delta}{h} \]
其中 $ M = \max_{\xi \in [x_0-h, x_0+h]} |f'''(\xi)| $，我们用它来近似 $ |f'''(x_0)| $ 的上界。$ \delta $ 是函数值的绝对误差水平。
- 寻找最优步长 \(h_{opt}\)：
总误差是 \(h\) 的函数。为了最小化总误差，我们对 \(h\) 求导（将总误差上界视为函数 \(g(h) = A h^2 + B / h\)，其中 \(A = M/6, B = \delta\)）：
\[ \frac{d}{dh} (A h^2 + B h^{-1}) = 2A h - B h^{-2} = 0 \]
\[ 2A h^3 = B \]
\[ h^3 = \frac{B}{2A} = \frac{\delta}{2(M/6)} = \frac{3\delta}{M} \]
\[ h_{opt} = \sqrt[3]{\frac{3\delta}{M}} \]
- 对最优步长的理解：
- \(h_{opt}\) 与 \(\delta^{1/3}\) 成正比，与 \(M^{1/3}\) 成反比。
- 如果函数值非常精确（\(\delta\) 很小），我们可以使用较小的 \(h\) 来减小截断误差。
- 如果函数的三阶导数值很大（函数变化剧烈），截断误差项 \(Mh^2\) 增长快，为了控制它，最优步长 \(h_{opt}\) 应该小一些。
- 如果函数值误差 \(\delta\) 很大，为了避免舍入误差被过度放大，最优步长 \(h_{opt}\) 应该大一些。
- 在最优点，截断误差和舍入误差的量级大致平衡：\(E_T \sim h_{opt}^2\)， \(E_R \sim \delta / h_{opt}\)，代入 \(h_{opt}\) 可得两者均约为 \(\delta^{2/3} M^{1/3}\) 量级。
步骤5：实例与验证
假设我们想计算 \(f(x) = e^x\) 在 \(x_0 = 0\) 处的导数，真实值 \(f'(0) = 1\)。假设 \(\delta = 10^{-16}\)（双精度舍入误差级别），并估计 \(M = |f'''(0)| = e^0 = 1\)。
- 计算最优步长：
\[ h_{opt} = \sqrt[3]{\frac{3 \times 10^{-16}}{1}} \approx \sqrt[3]{3 \times 10^{-16}} \approx 6.69 \times 10^{-6} \]
我们可以取 $ h \approx 10^{-5} $ 量级。
-
数值实验：
- 取 \(h = 10^{-2}\)：截断误差主导，约为 \(h^2/6 \approx 1.67\times 10^{-5}\)。
- 取 \(h = 10^{-5}\)：接近最优，总误差很小。
- 取 \(h = 10^{-8}\)：舍入误差开始主导，约为 \(\delta / h \approx 10^{-8}\)，精度反而不如 \(h = 10^{-5}\) 时。
- 取 \(h = 10^{-12}\)：舍入误差约为 \(10^{-4}\)，结果已严重失真。
这个简单的实验验证了最优步长的存在性：步长并非越小越好。
总结
通过本题目，我们深入剖析了数值微分三点公式（以中心差分公式为代表）的误差构成。截断误差（近似误差）与 \(h^2\) 成正比，而舍入误差（数据/计算误差）与 \(1/h\) 成正比。总误差是两者的叠加，存在一个使总误差最小的最优步长 \(h_{opt} \approx \sqrt[3]{3\delta / M}\)。在实际应用中，我们需要根据函数的光滑性（通过 \(M\) 体现）和数据的精度（\(\delta\) ）来合理选择步长，而非盲目地追求极小的 \(h\)。这一分析框架同样适用于其他数值微分公式（如前向/后向差分）的误差分析与步长优化。