数值微分中的五点中心差分公式：精度阶数与步长权衡分析
题目描述
数值微分是近似计算函数导数值的方法，其中差分公式是最基础的途径。我们考虑一个在区间上充分光滑的函数 \(f(x)\)。已知函数在某些离散点上的取值，要求用五点中心差分公式来近似计算其某一点的一阶导数值，并分析该公式的精度阶数（截断误差与步长的关系）以及在实际计算中如何选择合适的步长以避免舍入误差过大。
具体任务：
- 推导五点中心差分公式（即利用 \(x_{i-2}, x_{i-1}, x_i, x_{i+1}, x_{i+2}\) 五个等距节点处的函数值来估计 \(f'(x_i)\)）。
- 分析该公式的局部截断误差主项及其精度阶数。
- 讨论舍入误差的来源及其与步长的关系，给出选择步长的实用指导原则。
解题过程循序渐进讲解
1. 五点中心差分公式的推导
我们假设节点等距，步长为 \(h\)，即：
\[x_k = x_i + k h \quad (k = -2, -1, 0, 1, 2) \]
对应函数值记为 \(f_{i+k} = f(x_i + k h)\)。
思路：用这些点构造一个四次多项式 \(P_4(x)\) 来逼近 \(f(x)\)，然后对多项式求导，在 \(x_i\) 处得到导数的近似值。
更直接的方法是利用泰勒展开联立消去高阶项。将各点函数值在 \(x_i\) 处展开：
\[\begin{aligned} f_{i-2} &= f_i - 2h f'_i + \frac{(2h)^2}{2!} f''_i - \frac{(2h)^3}{3!} f'''_i + \frac{(2h)^4}{4!} f^{(4)}_i - \frac{(2h)^5}{5!} f^{(5)}_i + O(h^6) \\ f_{i-1} &= f_i - h f'_i + \frac{h^2}{2!} f''_i - \frac{h^3}{3!} f'''_i + \frac{h^4}{4!} f^{(4)}_i - \frac{h^5}{5!} f^{(5)}_i + O(h^6) \\ f_{i+1} &= f_i + h f'_i + \frac{h^2}{2!} f''_i + \frac{h^3}{3!} f'''_i + \frac{h^4}{4!} f^{(4)}_i + \frac{h^5}{5!} f^{(5)}_i + O(h^6) \\ f_{i+2} &= f_i + 2h f'_i + \frac{(2h)^2}{2!} f''_i + \frac{(2h)^3}{3!} f'''_i + \frac{(2h)^4}{4!} f^{(4)}_i + \frac{(2h)^5}{5!} f^{(5)}_i + O(h^6) \end{aligned} \]
其中 \(f_i = f(x_i)\)，\(f'_i = f'(x_i)\)，\(f^{(k)}_i = f^{(k)}(x_i)\)。
我们希望用这四个相邻点（和 \(f_i\) 本身）的线性组合来近似 \(f'_i\)，设：
\[f'_i \approx \frac{A f_{i-2} + B f_{i-1} + C f_i + D f_{i+1} + E f_{i+2}}{h} \]
代入泰勒展开并匹配 \(f_i, f'_i, f''_i, \dots\) 的系数。
系数匹配方程：
- \(f_i\) 的系数：\(A + B + C + D + E = 0\)（因为左边无常数项）。
- \(f'_i\) 的系数：\((-2A - B + D + 2E) = 1\)（因为分母有 \(h\)，右边应得 \(h \cdot f'_i\) 的系数为 1）。
- \(f''_i\) 的系数：\((4A + B + D + 4E)/2 = 0\)（我们希望消去二阶项）。
- \(f'''_i\) 的系数：\((-8A - B + D + 8E)/6 = 0\)（消去三阶项）。
- \(f^{(4)}_i\) 的系数：\((16A + B + D + 16E)/24 = 0\)（消去四阶项）。
解此线性方程组（五个未知数五个方程）：
由对称性可设 \(A = E\)，\(B = D\)，则方程简化为：
\[\begin{cases} 2A + 2B + C = 0 \\ (-2A - B + B + 2A) = 0 \quad\text{(自动满足)} \\ \text{实际上需仔细重列：} \end{cases} \]
重新按完整方程求解：
\[\begin{cases} A + B + C + B + A = 0 \implies 2A + 2B + C = 0 \\ -2A - B + B + 2A = 0 \quad\text{(自动成立，需用正确方程)} \end{cases} \]
正确列出 \(f'_i\) 系数方程（注意系数是 \(h\) 的一次项）：
左端 \(h f'_i\) 的系数来自泰勒展开中 \(h f'_i\) 项：
对 \(f_{i-2}\)：系数 \(-2A\)
对 \(f_{i-1}\)：系数 \(-B\)
对 \(f_{i+1}\)：系数 \(D\)
对 \(f_{i+2}\)：系数 \(2E\)
因此方程：\(-2A - B + D + 2E = 1\)。
利用对称性 \(A = E, B = D\) 代入：
\[-2A - B + B + 2A = 0 = 1 \]
矛盾！说明不能简单假设对称性为 \(A = E, B = D\)。实际上，因为要消去偶阶导数项（二阶、四阶），系数应呈奇对称：\(A = -E, B = -D\)。但这样会导致 \(f'_i\) 的系数为零。所以这里需要的是五点中心差分，系数应对称：\(A = -E, B = -D\) 不是我们要的。
正确对称性：系数应满足 \(A = -E, B = -D, C = 0\) 才可能是中心差分吗？不一定。实际上五点中心差分公式是对称的：系数满足 \(A = -E, B = -D, C=0\) 吗？我们来验证常见公式。
已知经典五点中心差分公式（一阶导数）为：
\[f'(x_i) \approx \frac{-f_{i+2} + 8f_{i+1} - 8f_{i-1} + f_{i-2}}{12h} \]
即：
\[A = 1/(12h),\quad B = -8/(12h) = -2/(3h),\quad C = 0,\quad D = 8/(12h) = 2/(3h),\quad E = -1/(12h)。 \]
记 \(A' = A h\) 等，则 \(A' = 1/12, B' = -8/12, C' = 0, D' = 8/12, E' = -1/12\)。
验证上述泰勒展开系数匹配：
\(A' + B' + C' + D' + E' = (1 -8 +0 +8 -1)/12 = 0\) ✅
\(f'_i\) 系数：\((-2A' -B' + D' + 2E') = (-2/12 + 8/12 + 8/12 - 2/12) = 12/12 = 1\) ✅
\(f''_i\) 系数：\((4A' + B' + D' + 4E')/2 = (4/12 -8/12 +8/12 -4/12)/2 = 0\) ✅
\(f'''_i\) 系数：\((-8A' - B' + D' + 8E')/6 = (-8/12 +8/12 +8/12 -8/12)/6 = 0\) ✅
\(f^{(4)}_i\) 系数：\((16A' + B' + D' + 16E')/24 = (16/12 -8/12 +8/12 -16/12)/24 = 0\) ✅
\(f^{(5)}_i\) 系数：\((-32A' - B' + D' + 32E')/120 = (-32/12 +8/12 +8/12 -32/12)/120 = (-48/12)/120 = -4/120 = -1/30\)。
因此局部截断误差的主项为：
\[R = \frac{h^5}{5!} f^{(5)}_i \times (-1/30) \times (120) \quad\text{？} \]
更准确：把 \(f_{i-2}, f_{i-1}, f_{i+1}, f_{i+2}\) 的泰勒展开代入公式，合并后 \(f^{(5)}_i\) 的系数是：
\[\frac{(-32A' - B' + D' + 32E')}{120} h^5 f^{(5)}_i = \frac{-48/12}{120} h^5 f^{(5)}_i = \frac{-4}{120} h^5 f^{(5)}_i = -\frac{1}{30} h^5 f^{(5)}_i。 \]
2. 局部截断误差与精度阶数
局部截断误差 \(T_i\) 为：
\[T_i = f'(x_i) - \frac{-f_{i+2} + 8f_{i+1} - 8f_{i-1} + f_{i-2}}{12h} \]
由泰勒展开推导可知，所有低于五阶的导数项都被精确消去，第一个非零误差项来自五阶导数：
\[T_i = -\frac{h^4}{30} f^{(5)}(x_i) + O(h^6) \]
因此该公式是 四阶精度 的（误差与 \(h^4\) 成正比）。
3. 舍入误差与步长选择
实际计算中，函数值 \(f_k\) 有舍入误差，记为 \(\epsilon_k\)，通常 \(|\epsilon_k| \le \epsilon_{\text{mach}} \cdot |f_k|\)，其中 \(\epsilon_{\text{mach}}\) 是机器精度。
五点公式的舍入误差放大系数（即公式中系数的绝对值之和除以 \(h\)）为：
\[\text{系数绝对和} = \frac{1 + 8 + 8 + 1}{12h} = \frac{18}{12h} = \frac{1.5}{h} \]
因此舍入误差引起的导数误差大约为：
\[E_{\text{round}} \approx \frac{1.5 \epsilon_{\text{mach}} \cdot M}{h} \]
其中 \(M\) 是 \(|f|\) 的量级。
总误差 ≈ 截断误差 + 舍入误差：
\[E_{\text{total}} \approx \frac{h^4}{30} |f^{(5)}| + \frac{1.5 \epsilon_{\text{mach}} M}{h} \]
对 \(h\) 求导找极小点：
\[\frac{dE_{\text{total}}}{dh} \approx \frac{4h^3}{30} |f^{(5)}| - \frac{1.5 \epsilon_{\text{mach}} M}{h^2} = 0 \]
\[ \frac{2}{15} h^5 |f^{(5)}| \approx 1.5 \epsilon_{\text{mach}} M \]
\[ h_{\text{opt}} \approx \left( \frac{11.25 \epsilon_{\text{mach}} M}{|f^{(5)}|} \right)^{1/5} \]
若 \(M \sim 1\)，\(|f^{(5)}| \sim 1\)，\(\epsilon_{\text{mach}} \approx 10^{-16}\)（双精度），则：
\[h_{\text{opt}} \approx (11.25 \times 10^{-16})^{1/5} \approx (1.125\times 10^{-15})^{0.2} \approx 10^{-3} \]
数量级在 \(10^{-3}\) 左右。
4. 实用指导
- 精度阶数：五点中心差分公式是四阶精度的，比三点中心差分（二阶精度）更精确。
- 步长选择：步长不宜太小，否则舍入误差占主导；不宜太大，否则截断误差大。可用上述 \(h_{\text{opt}}\) 公式估计，或通过试验取误差最小的 \(h\)。
- 适用条件：函数需充分光滑（至少五阶可导），且节点等距。
通过以上步骤，我们完成了五点中心差分公式的推导、误差分析和步长选择讨论。