好的，作为数值积分与微分领域的专家，我将为你讲解一个新颖且基础的题目。
数值微分：Richardson 外推法在中心差分公式中的应用与精度提升
题目描述
数值微分的核心任务是用函数在离散点的值来近似计算其导数值。最常用的方法之一是中心差分公式。例如，计算函数 \(f(x)\) 在点 \(x_0\) 处的一阶导数，常用三点中心差分公式：
\[f'(x_0) \approx D(h) = \frac{f(x_0 + h) - f(x_0 - h)}{2h} \]
其中，\(h\) 为步长。这个公式的截断误差（由泰勒展开忽略的高阶项引起）是 \(O(h^2)\) 阶的，即误差大致与 \(h^2\) 成正比。
我们的问题是：如何在不减小步长 \(h\) 的前提下（因为过小的 \(h\) 会放大舍入误差），显著提高上述中心差分公式的精度？
解题思路是运用 Richardson 外推法，通过组合不同步长的计算结果，消去误差展开式中的低阶误差项，从而得到更高阶精度的近似。
解题过程循序渐进讲解
步骤1：分析中心差分公式的误差结构
首先，我们对 \(f(x_0 + h)\) 和 \(f(x_0 - h)\) 进行泰勒展开：
\[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) \]
两边同除以 \(2h\)，得到我们使用的近似公式 \(D(h)\)：
\[D(h) = \frac{f(x_0 + h) - f(x_0 - h)}{2h} = f'(x_0) + \frac{h^2}{6} f'''(x_0) + O(h^4) \]
关键结论：近似值 \(D(h)\) 与真值 \(f'(x_0)\) 之间的误差可以写成关于 \(h^2\) 的幂级数：
\[D(h) = f'(x_0) + C_2 h^2 + C_4 h^4 + C_6 h^6 + \dots \]
这里 \(C_2 = \frac{f'''(x_0)}{6}\)，\(C_4, C_6, \dots\) 是由更高阶导数决定的常数。没有 \(h\) 的奇数次幂项（因为中心差分的对称性消除了它们）。
步骤2：构造 Richardson 外推
Richardson 外推的精髓在于：如果我们用两个不同的步长（例如 \(h\) 和 \(h/2\)）分别计算 \(D(h)\) 和 \(D(h/2)\)，那么它们的误差展开式具有相同的结构，但系数不同。
- 对于步长 \(h\)：
\[ D(h) = f'(x_0) + C_2 h^2 + C_4 h^4 + C_6 h^6 + \dots \quad \text{(1)} \]
- 对于步长 \(h/2\)：
\[ D(h/2) = f'(x_0) + C_2 \left(\frac{h}{2}\right)^2 + C_4 \left(\frac{h}{2}\right)^4 + C_6 \left(\frac{h}{2}\right)^6 + \dots \]
\[ = f'(x_0) + \frac{C_2}{4} h^2 + \frac{C_4}{16} h^4 + \frac{C_6}{64} h^6 + \dots \quad \text{(2)} \]
观察(1)和(2)式，我们发现它们都含有我们想要求解的 \(f'(x_0)\)，并且误差项都从 \(h^2\) 项开始。为了消去误差中占主导地位的 \(h^2\) 项，我们可以将(1)式乘以某个系数减去(2)式。
具体地，我们希望找到一个组合 \(a \cdot D(h) + b \cdot D(h/2)\)，使得 \(h^2\) 项的系数为零：
\[a \cdot (C_2 h^2) + b \cdot (\frac{C_2}{4} h^2) = 0 \]
同时，为了保证组合结果仍然近似于 \(f'(x_0)\)，我们需要：
\[a + b = 1 \]
这是一个简单的二元一次方程组：
\[\begin{cases} a + b = 1 \\ a + \frac{b}{4} = 0 \end{cases} \]
解得：\(a = -\frac{1}{3}, \quad b = \frac{4}{3}\)。
步骤3：得到高阶精度公式
将解得的系数代入组合：
\[f'(x_0) \approx D_1(h) = \frac{4}{3} D(h/2) - \frac{1}{3} D(h) \]
这里的下标 \(_1\) 表示进行了一次外推。
现在我们来分析新公式 \(D_1(h)\) 的误差。将(1)式和(2)式代入：
\[D_1(h) = \frac{4}{3} \left[ f'(x_0) + \frac{C_2}{4} h^2 + \frac{C_4}{16} h^4 + \dots \right] - \frac{1}{3} \left[ f'(x_0) + C_2 h^2 + C_4 h^4 + \dots \right] \]
合并同类项：
- \(f'(x_0)\) 项：\(\frac{4}{3} - \frac{1}{3} = 1\)。
- \(h^2\) 项：\(\frac{4}{3} \cdot \frac{C_2}{4} - \frac{1}{3} \cdot C_2 = \frac{C_2}{3} - \frac{C_2}{3} = 0\)。
- \(h^4\) 项：\(\frac{4}{3} \cdot \frac{C_4}{16} - \frac{1}{3} \cdot C_4 = \frac{C_4}{12} - \frac{C_4}{3} = -\frac{C_4}{4}\)。
因此，
\[D_1(h) = f'(x_0) - \frac{C_4}{4} h^4 + O(h^6) \]
奇迹发生了！ 新公式 \(D_1(h)\) 的误差从原来的 \(O(h^2)\) 阶提升到了 \(O(h^4)\) 阶，精度实现了飞跃。
步骤4：进一步外推与一般化
这个过程可以继续迭代。如果我们再用步长 \(h/4\) 计算 \(D(h/4)\)，然后对 \(D(h/2)\) 和 \(D(h/4)\) 应用同样的外推逻辑（但目标是消去 \(h^4\) 项），可以得到误差为 \(O(h^6)\) 阶的公式 \(D_2(h)\)。
在实际编程中，这通常表现为构建一个外推表格（类似于龙贝格积分表）：
步长 \ 外推阶 | 0阶(O(h^2)) 1阶(O(h^4)) 2阶(O(h^6))
-----------------------------------------------------
h | D(h)
h/2 | D(h/2) D_1(h)
h/4 | D(h/4) D_1(h/2) D_2(h)
... | ... ... ...
表中元素的递推公式为：
\[D_{k}(h) = \frac{4^k D_{k-1}(h/2) - D_{k-1}(h)}{4^k - 1} \]
其中 \(D_0(h) = D(h)\) 就是最初的中心差分公式。每向右一列，精度提高两阶。
总结
通过以上步骤，我们系统地解决了如何提升中心差分公式精度的问题：
- 分析误差：利用泰勒展开，明确了中心差分公式 \(D(h)\) 的误差是 \(h^2\) 的偶次幂级数。
- 设计外推：通过组合两个不同步长（\(h\) 和 \(h/2\)）的 \(D(h)\) 值，建立方程组，消去了误差的主项 \(h^2\) 项。
- 获得高阶公式：得到了新的近似公式 \(D_1(h) = \frac{4D(h/2) - D(h)}{3}\)，其精度从 \(O(h^2)\) 提升至 \(O(h^4)\)。
- 迭代推广：此过程可以系统化、迭代进行，构建外推表格，理论上可以获得任意高阶精度的数值微分公式，且无需使用更小的步长来牺牲数值稳定性。
这种方法的核心价值在于，它用计算量（多算几个步长的函数值并进行组合）巧妙地换取了精度阶数的显著提升，是数值计算中“用智慧替代蛮力”的典范。