数值微分的五点中心差分公式：高精度构造及其截断误差分析
数值微分是计算函数导数近似值的重要方法。当函数没有解析导数表达式，或导数难以直接计算时，数值微分提供了一种基于离散函数值的实用工具。在众多数值微分公式中，五点中心差分公式因其较高的精度和对称性而被广泛使用。今天，我将为你详细讲解五点中心差分公式的推导、构造过程以及其截断误差的精确分析。
问题描述
给定一个足够光滑的函数 \(f(x)\)，我们希望计算它在某一点 \(x_0\) 的一阶导数 \(f'(x_0)\) 的近似值。我们能够获取的信息是函数在一些离散点上的值，这些点通常以等间距分布，步长为 \(h\)。具体地，已知 \(f(x_0 - 2h), f(x_0 - h), f(x_0), f(x_0 + h), f(x_0 + 2h)\) 这五个点的函数值。我们的目标是：仅利用这五个点的函数值，构造一个高精度的数值微分公式来近似 \(f'(x_0)\)，并分析该公式的截断误差（即理论误差）。
解题过程
步骤1：核心思想与泰勒展开
数值微分的基本思想是利用函数在某点附近的泰勒展开式，将导数与函数值的线性组合联系起来。对于五点中心差分公式，我们利用点 \(x_0\) 及其左右对称的四个点进行展开，以获得更高的精度。
我们定义：
- \(x_0\) 为求导点。
- 步长 \(h > 0\)。
- 已知点：\(x_{-2} = x_0 - 2h\), \(x_{-1} = x_0 - h\), \(x_0\), \(x_1 = x_0 + h\), \(x_2 = x_0 + 2h\)。
假设 \(f(x)\) 在包含这些点的区间内无限次可微（至少五阶连续可导），我们在 \(x_0\) 处进行泰勒展开：
- \(f(x_0 + h) = f(x_0) + h f'(x_0) + \frac{h^2}{2!} f''(x_0) + \frac{h^3}{3!} f^{(3)}(x_0) + \frac{h^4}{4!} f^{(4)}(x_0) + \frac{h^5}{5!} f^{(5)}(x_0) + O(h^6)\)
- \(f(x_0 - h) = f(x_0) - h f'(x_0) + \frac{h^2}{2!} f''(x_0) - \frac{h^3}{3!} f^{(3)}(x_0) + \frac{h^4}{4!} f^{(4)}(x_0) - \frac{h^5}{5!} f^{(5)}(x_0) + O(h^6)\)
- \(f(x_0 + 2h) = f(x_0) + 2h f'(x_0) + \frac{(2h)^2}{2!} f''(x_0) + \frac{(2h)^3}{3!} f^{(3)}(x_0) + \frac{(2h)^4}{4!} f^{(4)}(x_0) + \frac{(2h)^5}{5!} f^{(5)}(x_0) + O(h^6)\)
- \(f(x_0 - 2h) = f(x_0) - 2h f'(x_0) + \frac{(2h)^2}{2!} f''(x_0) - \frac{(2h)^3}{3!} f^{(3)}(x_0) + \frac{(2h)^4}{4!} f^{(4)}(x_0) - \frac{(2h)^5}{5!} f^{(5)}(x_0) + O(h^6)\)
其中，\(O(h^6)\) 表示六次及更高次幂的项。
步骤2：构造线性组合以消去无关项
我们的目标是找到一个形如：
\[f'(x_0) \approx \frac{A f(x_0 - 2h) + B f(x_0 - h) + C f(x_0) + D f(x_0 + h) + E f(x_0 + 2h)}{h} \]
的公式，其中 \(A, B, C, D, E\) 是待定系数。
为了确定这些系数，我们将四个泰勒展开式代入上述近似公式的右边，然后与左边的 \(f'(x_0)\) 进行比较，并要求两边的对应项系数匹配。具体来说：
将 \(f(x_0 \pm h)\) 和 \(f(x_0 \pm 2h)\) 的泰勒展开代入右边，并按 \(f(x_0), f'(x_0), f''(x_0), \dots\) 的幂次整理。右边表达式为：
\[\frac{1}{h} \left[ A f(x_0 - 2h) + B f(x_0 - h) + C f(x_0) + D f(x_0 + h) + E f(x_0 + 2h) \right] \]
我们要求这个表达式在 \(h \to 0\) 时精确等于 \(f'(x_0)\)。为了使公式精确，我们希望：
- \(f(x_0)\) 的系数之和为零。
- \(f'(x_0)\) 的系数之和为 \(h\)（因为右边除以了 \(h\)，所以其系数和应为 \(h\) 才能与左边的 \(f'(x_0)\) 匹配）。
- 尽可能多地消去高阶项（如 \(f''(x_0), f^{(3)}(x_0)\) 等），以提高精度。
将泰勒展开式代入并整理后，得到各阶导数的系数方程：
- 常数项（\(f(x_0)\) 的系数）：
\[ A + B + C + D + E = 0 \]
- 一阶导项（\(f'(x_0)\) 的系数）：
\[ (-2A - B + D + 2E) h = h \quad \Rightarrow \quad -2A - B + D + 2E = 1 \]
（注意：等式右边是 \(f'(x_0)\)，其系数为1，但我们的表达式整体除以 \(h\)，所以这里 \(f'(x_0)\) 的系数应为 \(h\) 才能经 \(1/h\) 后变成1。）
- 二阶导项（\(f''(x_0)\) 的系数）：
\[ \left( \frac{(2h)^2}{2!} A + \frac{h^2}{2!} B + \frac{h^2}{2!} D + \frac{(2h)^2}{2!} E \right) \frac{1}{h} = \frac{h}{2} (4A + B + D + 4E) \]
我们希望这项为零以提高精度，所以：
\[ 4A + B + D + 4E = 0 \]
- 三阶导项（\(f^{(3)}(x_0)\) 的系数）：
\[ \left( -\frac{(2h)^3}{3!} A - \frac{h^3}{3!} B + \frac{h^3}{3!} D + \frac{(2h)^3}{3!} E \right) \frac{1}{h} = \frac{h^2}{6} (-8A - B + D + 8E) \]
希望为零：
\[ -8A - B + D + 8E = 0 \]
- 四阶导项（\(f^{(4)}(x_0)\) 的系数）：
\[ \left( \frac{(2h)^4}{4!} A + \frac{h^4}{4!} B + \frac{h^4}{4!} D + \frac{(2h)^4}{4!} E \right) \frac{1}{h} = \frac{h^3}{24} (16A + B + D + 16E) \]
希望为零：
\[ 16A + B + D + 16E = 0 \]
我们得到了五个方程，但理论上我们可以保留四阶项来构造一个更高精度的公式。然而，我们只有五个未知数 \(A, B, C, D, E\)。为了获得尽可能高的精度，我们通常要求消去直到四阶导的项（这样截断误差将是 \(O(h^4)\) 量级）。让我们检查方程数：
- 我们有5个未知数。
- 我们列出了5个方程：常数项、一阶、二阶、三阶、四阶项为零（实际上常数项、二阶、三阶、四阶设为零，一阶设为1）。
现在解这个线性方程组。
步骤3：求解系数
方程组为：
\[\begin{cases} A + B + C + D + E = 0 \quad &(1)\\ -2A - B + D + 2E = 1 \quad &(2)\\ 4A + B + D + 4E = 0 \quad &(3)\\ -8A - B + D + 8E = 0 \quad &(4)\\ 16A + B + D + 16E = 0 \quad &(5) \end{cases} \]
首先，观察(3)、(4)、(5)方程。它们都是关于 \((A+E)\) 和 \((B+D)\) 的线性组合。实际上，将(3)和(5)相减，可以简化。
由(3): \(4(A+E) + (B+D) = 0\)
由(5): \(16(A+E) + (B+D) = 0\)
将两式相减：\(12(A+E) = 0 \Rightarrow A+E = 0\)
代入(3)：\(0 + (B+D) = 0 \Rightarrow B+D = 0\)
所以 \(E = -A\)，\(D = -B\)。
现在利用(4)检查一致性：将 \(D=-B, E=-A\) 代入(4)：
\[-8A - B + (-B) + 8(-A) = -8A - B - B - 8A = -16A - 2B = 0 \Rightarrow 8A + B = 0 \quad (4') \]
利用(2)：将 \(D=-B, E=-A\) 代入(2)：
\[-2A - B + (-B) + 2(-A) = -2A - B - B - 2A = -4A - 2B = 1 \Rightarrow 2A + B = -\frac{1}{2} \quad (2') \]
现在我们有两个方程：
\[8A + B = 0 \quad \text{和} \quad 2A + B = -\frac{1}{2} \]
相减：\((8A+B) - (2A+B) = 0 - (-\frac{1}{2}) \Rightarrow 6A = \frac{1}{2} \Rightarrow A = \frac{1}{12}\)
代入 \(8A+B=0\)：\(8 \times \frac{1}{12} + B = 0 \Rightarrow \frac{2}{3} + B = 0 \Rightarrow B = -\frac{2}{3}\)
于是：
\[A = \frac{1}{12}, \quad B = -\frac{2}{3}, \quad E = -A = -\frac{1}{12}, \quad D = -B = \frac{2}{3} \]
最后利用(1)求 \(C\)：
\[A + B + C + D + E = \frac{1}{12} - \frac{2}{3} + C + \frac{2}{3} - \frac{1}{12} = C = 0 \]
因此，所有系数为：
\[A = \frac{1}{12}, \quad B = -\frac{2}{3}, \quad C = 0, \quad D = \frac{2}{3}, \quad E = -\frac{1}{12} \]
步骤4：得到五点中心差分公式
将系数代入近似式：
\[f'(x_0) \approx \frac{ \frac{1}{12} f(x_0 - 2h) - \frac{2}{3} f(x_0 - h) + 0 \cdot f(x_0) + \frac{2}{3} f(x_0 + h) - \frac{1}{12} f(x_0 + 2h) }{h} \]
通常写作更简洁的形式：
\[\boxed{ f'(x_0) \approx \frac{ f(x_0 - 2h) - 8f(x_0 - h) + 8f(x_0 + h) - f(x_0 + 2h) }{12h} } \]
这个公式就是五点中心差分公式（实际上用了四个点，因为中心点 \(f(x_0)\) 的系数为零）。
步骤5：截断误差分析
为了分析截断误差，我们需要检查被我们忽略的高阶项。回顾泰勒展开，当我们代入系数后，所有直到四阶导的项都被精确消去了。接下来看五阶导项。
计算五阶导项 \(f^{(5)}(x_0)\) 的系数：
- 来自 \(f(x_0 - 2h)\)：\(- \frac{(2h)^5}{5!} \times A = - \frac{32h^5}{120} \times \frac{1}{12} = -\frac{32}{1440} h^5\)
- 来自 \(f(x_0 - h)\)：\(- \frac{h^5}{120} \times B = - \frac{h^5}{120} \times (-\frac{2}{3}) = +\frac{2}{360} h^5\)
- 来自 \(f(x_0 + h)\)：\(+ \frac{h^5}{120} \times D = + \frac{h^5}{120} \times \frac{2}{3} = +\frac{2}{360} h^5\)
- 来自 \(f(x_0 + 2h)\)：\(+ \frac{(2h)^5}{120} \times E = + \frac{32h^5}{120} \times (-\frac{1}{12}) = -\frac{32}{1440} h^5\)
将这些系数相加（并除以 \(h\) 因为公式右边有 \(1/h\)）：
\[\frac{1}{h} \times h^5 \left( -\frac{32}{1440} + \frac{2}{360} + \frac{2}{360} - \frac{32}{1440} \right) \]
简化：\(\frac{2}{360} = \frac{1}{180} = \frac{8}{1440}\)，所以括号内为：
\[\left( -\frac{32}{1440} + \frac{8}{1440} + \frac{8}{1440} - \frac{32}{1440} \right) = -\frac{48}{1440} = -\frac{1}{30} \]
因此五阶导项的总贡献为：
\[\frac{1}{h} \times h^5 \times \left(-\frac{1}{30}\right) = -\frac{h^4}{30} f^{(5)}(x_0) \]
这里我们假设各点的五阶导在 \(x_0\) 处相同（因为步长很小，且函数光滑）。实际上，更严格的误差分析需要利用泰勒余项。但通常我们取主导项作为截断误差的估计。
于是，完整的表达式为：
\[f'(x_0) = \frac{ f(x_0 - 2h) - 8f(x_0 - h) + 8f(x_0 + h) - f(x_0 + 2h) }{12h} - \frac{h^4}{30} f^{(5)}(\xi) \]
其中 \(\xi\) 是位于 \([x_0-2h, x_0+2h]\) 之间的某一点（根据泰勒展开的拉格朗日余项形式）。因此，截断误差为 \(-\frac{h^4}{30} f^{(5)}(\xi)\)，这表明五点中心差分公式具有四阶精度（误差与 \(h^4\) 成正比）。
总结
- 五点中心差分公式（实际使用四个对称点）：
\[ f'(x_0) \approx \frac{ f(x_0 - 2h) - 8f(x_0 - h) + 8f(x_0 + h) - f(x_0 + 2h) }{12h} \]
- 精度：四阶精度，截断误差为 \(O(h^4)\)，具体主项为 \(-\frac{h^4}{30} f^{(5)}(\xi)\)。
- 优点：相比两点或三点公式，精度显著提高，特别适用于对导数精度要求较高的场合。
- 注意事项：当 \(h\) 过小时，舍入误差可能增大（因函数值相减造成有效数字损失），因此需权衡步长选择。
通过以上循序渐进的推导和分析，你应该对五点中心差分公式的构造原理和误差特性有了清晰的理解。