数值微分与Richardson外展
在数学分析中，我们用导数来描述函数在某一点处的变化率。对于一个显式给出的连续函数 \(f(x)\)，我们可以通过求导法则写出它的导函数 \(f'(x)\)。这种逻辑推导的方法，被称为解析微分。 但是在工程实践中，我们通常无法获得函数的解析表达式，更无法直接得到它的导函数。大多数情况下，我们面对仅仅是一系列离散的采样数据。在一些高级的数学模型中， 函数 \(f(x)\) 的表达式可能及其复杂，手撕一个解析的导函数出来很困难。面对这些无法直接求导的离散数据或黑盒函数，我们该如何分析它的变化率？这正是数值微分要解决的问题。
数值微分研究的是如何利用有限的、离散的函数值信息，去近似计算函数在某一点处的导数或微商。基于拉格朗日多项式以及泰勒展开， 我们可以写出十分简洁的数值微分计算式。但是在这里我认为，计算表达式不应该是我们关注的重点。我们应该认真体会一下，在数值计算中的截断误差和舍入误差。
1. 前向微分与后向微分
函数 \(f(x)\) 在 \(x_0\) 处的导数定义为
$$ f'(x_0) = \lim_{h \to 0} \frac{f(x_0 + h) - f(x_0)}{h} $$这个定义在数学分析中很严谨。但很不幸，在计算机的离散数字世界里， 我们无法真正执行一个趋近于 0 的极限过程。我们能做的，只是给定一个足够小的步长 \(h\)，用有限的代数运算\(\frac{f(x_0 + h) - f(x_0)}{h}\) 来近似表示导数 \(f'(x_0)\)。
实际上这是一个很粗糙的近似。假设函数 \(f(x)\) 在闭区间 \([a,b]\) 上是 \(C^2\) 连续的，\(x_1 = x_0 + h\) 并且 \(x_0,x_1 \in (a,b)\)。 根据拉格朗日多项式的逼近误差公式，可以把 \(f(x)\) 写成一个多项式 \(P_1(x)\)一个余项:
$$ \begin{aligned} f(x) & = \frac{f(x_0)(x-x_1)}{x_0-x_1} + \frac{f(x_1)(x - x_0)}{x_1-x_0} + \frac{(x-x_0)(x-x_1)}{2!}f''(\varepsilon(x)) \\ & = \frac{f(x_0)(x - x_1)}{-h} + \frac{f(x_0 + h)(x - x_0)}{h} + \frac{(x-x_0)(x-x_1)}{2!}f''(\varepsilon(x)) \end{aligned} $$对上式在 \(x_0\) 处求导有:
$$ \begin{aligned} f'(x_0) & = \frac{f(x_0+h) - f(x_0)}{h} + \frac{(x_0 - x_1) + (x_0 - x_0)}{2}f''(\varepsilon(x_0)) + \frac{(x_0-x_0)(x_0-x_1)}{2}\frac{d {f''(\varepsilon(x_0))}}{dx} \\ & = \frac{f(x_0+h) - f(x_0)}{h} - \frac{h}{2}f''(\xi) \end{aligned} $$上式中的第一项就是我们微分近似公式，如果 \(h > 0\) 就是前向微分，\(h > 0\) 则是后向微分。 被我们抛弃掉的余项 \(- \frac{h}{2}f''(\xi)\)，就是数值计算中的截断误差(Truncation Error)。 它与 \(h\) 呈一阶的线性关系，我们可以用大\(O\)记号 \(O(h)\) 笼统的替代之:
$$ \begin{equation} f'(x_0) = \frac{f(x_0+h) - f(x_0)}{h} + O(h) \end{equation} $$上述的推导过程也可以通过泰勒展开的形式得到。如果函数 \(f(x)\) 在 \(x_{0}\) 附近足够光滑，那么它在 \(x_0 + h\) 处的函数值可以展开为:
$$ f(x_{0}+h)=f(x_{0})+hf^{\prime }(x_{0})+\frac{h^{2}}{2!}f^{\prime \prime }(x_{0})+\frac{h^{3}}{3!}f^{\prime \prime \prime }(x_{0})+\dots $$将 \(f'(x_0)\) 移到等号左侧有:
$$ \begin{equation}\label{f2} \begin{array}{rrcl} & hf^{\prime }(x_{0}) & = & f(x_{0}+h)-f(x_{0})-\frac{h^{2}}{2!}f^{\prime \prime }(x_{0})-\frac{h^{3}}{3!}f^{\prime \prime \prime }(x_{0})-\dots\\ \Longrightarrow & f^{\prime }(x_{0}) & = & \frac{f(x_{0}+h)-f(x_{0})}{h} - \left[\frac{h}{2}f''(x_0) + \frac{h^2}{6}f'''(x_0) + \cdots \right] \\ & & = & \frac{f(x_{0}+h)-f(x_{0})}{h} + O(h) \end{array} \end{equation} $$我们通常说 \(O(h)\) 为一阶精度，可以预期，如果缩小步长 \(h\)，那么截断误差也会等比例的缩小。但是在计算机的离散世界里面，我们只能以固定的位数来表示浮点数，这就带来了舍入误差。 近似项 \(\frac{f(x_{0}+h)-f(x_{0})}{h}\) 还要除以 \(h\)，如果步长过小，\(h\) 的舍入误差还会被除法放大。所以数值微分是不稳定的。
2. \(n+1\) 点法
截断误差与舍入误差的这种矛盾，让我们陷入了两难。在浮点数位数确定的情况下，我们没有什么好的手段改进舍入误差。要想进一步提高算法的精度，我们只能想办法减少截断误差的影响。 这可以通过更多的数据点来做到。
如果点 \(\{x_0, \cdots, x_n\}\) 是闭区间 \([a, b]\) 上连续函数 \(f(x) \in C^{n+1}[a,b]\) 的 \(n+1\) 个不同的采样点。记拉格朗日基函数 \(L_i(x) = \prod_{i \neq k}\frac{x - x_i}{x_k - x_i}\)， 那么函数 \(f(x)\) 可以写为:
$$ f(x) = \sum_{i = 0}^n f(x_i)L_i(x) + \frac{f^{(n+1)}(\varepsilon(x))}{(n+1)!} \prod_{i=0}^n(x - x_i) $$其中 \(\varepsilon(x) \in [a,b]\) 是一个关于 \(x\) 的数值。对上式求导有
$$ f'(x) = \sum_{i=0}^n f(x_i)L_i'(x) + \frac{1}{(n+1)!}\left[\frac{d {\prod_{i=0}^n(x - x_i)}}{dx}f^{(n+1)}(\varepsilon(x)) + \frac{d f^{(n+1)}(\varepsilon(x))}{dx}\prod_{i=0}^n(x - x_i)\right] $$对于 \(x_j \in \{x_0, \cdots, x_n\}\)，有
$$ \begin{aligned} f'(x_j) & = \sum_{i=0}^n f(x_i)L_i'(x_j) + \frac{1}{(n+1)!}\left[\frac{d {\prod_{i=0}^n(x - x_i)}}{dx}f^{(n+1)}(\varepsilon(x))\right] \\ & = \sum_{i=0}^n f(x_i)L_i'(x_j) + \frac{f^{(n+1)}(\varepsilon(x_j))}{(n+1)!}\prod_{i = 0, i \neq j}^n(x_j - x_i) \end{aligned} $$上式就是 \(f'(x_j)\) 的 \(n+1\) 点的近似公式。其中 \(\frac{f^{(n+1)}(\varepsilon(x_j))}{(n+1)!}\prod_{i = 0, i \neq j}^n(x_j - x_i)\) 是截断余项。 如果 \(n+1\) 个点 \(\{x_0, \cdots, x_n\}\) 是等距采样的，采样步长是 \(h\)，那么截断余项将于 \(h^n\) 呈正比关系。 所以我们用 \(\sum_{i=0}^n f(x_i)L_i'(x_j)\) 来近似表示 \(f'(x_j)\) 时，截断误差是 \(O(h^n)\)，即:
$$ \begin{equation}\label{f3} f'(x_j) = \sum_{i=0}^n f(x_i)L_i'(x_j) + O(h^n) \end{equation} $$我们说 \(n+1\) 点的数值微分具有n 阶精度。采样点数的增加虽然可以提高计算精度，但同时它也增加了计算量。 而且表达式中基本都是累乘、累加的运算，点数的增加，可能导致计算误差的累积，最终的计算精度可能并不会像我们预期的那样总是得到提升。 所以，一般我们都是使用三点法或者五点法来近似计算。套用上式(\(\ref{f3}\))，经过一通化简之后有:
$$ \begin{equation}\label{f4} \begin{cases} f'(x_0) = \frac{1}{2h}\left[f(x_0 + h) - f(x_0 - h)\right] - \frac{h^2}{6}f^{(3)}(\xi_0), & 三点中心公式 \\ f'(x_0) = \frac{1}{2h}\left[-3f(x_0) + 4f(x_0 + h) - f(x_0 + 2h)\right] - \frac{h^2}{3}f^{(3)}(\xi_1), & 三点端点公式 \\ \end{cases} \end{equation} $$上式中 \(\xi_0\) 位于 \(x_0 + h\) 和 \(x_0 - h\) 之间，\(\xi_1\) 位于 \(x_0\) 和 \(x_0 + 2h\) 之间。
$$ \begin{equation}\label{f5} \begin{cases} f'(x_0) = \frac{1}{12h}\left[f(x_0 - 2h) - 8f(x_0 - h) + 8f(x_0 + h) - f(x_0 + 2h)\right] - \frac{h^4}{30}f^{(5)}(\xi_2), & 五点中心公式 \\ f'(x_0) = \frac{1}{12h}\left[-25f(x_0) + 48f(x_0 + h) - 36f(x_0 + 2h) + 16f(x_0 + 3h) - 3f(x_0 + 4h)\right] + \frac{h^4}{5}f^{(5)}(\xi_3), & 五点端点公式 \\ \end{cases} \end{equation} $$上式中 \(\xi_2\) 位于 \(x_0 + 2h\) 和 \(x_0 - 2h\) 之间，\(\xi_3\) 位于 \(x_0\) 和 \(x_0 + 4h\) 之间。
3. Richardson外展
Richardson 外展是数值分析中非常经典且强大的思想，它通过巧妙组合低精度的计算结果，来获得更高的精度，无需推导复杂的公式。 我们再来看一下通过泰勒展开得到的前向微分公式 (\(\ref{f2}\)):
$$ \begin{equation}\label{f6} f^{\prime }(x_{0}) = \underbrace{\frac{f(x_{0}+h)-f(x_{0})}{h}}_{N(h)} - \left[\frac{h}{2}f''(x_0) + \frac{h^2}{6}f'''(x_0) + \cdots \right] \end{equation} $$如果我们取 \(h / 2\) 的步长，代入上式有:
$$ \begin{equation}\label{f7} f^{\prime }(x_{0}) = \underbrace{\frac{f(x_{0}+\frac{h}{2})-f(x_{0})}{\frac{h}{2}}}_{N(\frac{h}{2})} - \left[\frac{h}{4}f''(x_0) + \frac{h^2}{24}f'''(x_0) + \cdots \right] \end{equation} $$我们将式(\(\ref{f7}\))乘以 2 减去式(\(\ref{f6}\)) 有:
$$ \begin{equation}\label{f8} f^{\prime }(x_{0}) = 2N(\frac{h}{2}) - N(h) + \left[\frac{h^2}{12}f'''(x_0) + \cdots \right] \end{equation} $$显然，式(\(\ref{f8}\))的截断误差就是 \(O(h^2)\) 的了，也就达到了二阶精度。整个过程十分简单明了，先计算 \(N(h)\)，再减半步长计算 \(N(h/2)\)，最后加权消元。 我们还可以以类似的方法，消去 \(\frac{h^2}{12}f'''(x_0)\) 得到三阶精度。这个操作理论上可以一直持续下去，得到更高阶的精度。 但计算机的离散世界里面，如果盲目地一直外推下去，会发现不仅精度不再提升，计算结果反而会彻底崩溃。这仍然是舍入误差与截断误差的矛盾。 因为舍入误差的存在，步长就不能一直减半，所以也不能一直消除截断误差。
由于数值微分本身就是一个不稳定的问题，所以 Richardson 外展所能发挥的作用有限。但在数值积分问题中，Richardson 外展将大放异彩。