[科普中国]-数值微分

简介

数值微分(numerical differentiation)根据函数在一些离散点的函数值，推算它在某点的导数或高阶导数的近似值的方法。通常用差商代替微商，或者用一个能够近似代替该函数的较简单的可微函数（如多项式或样条函数等）的相应导数作为能求导数的近似值。例如一些常用的数值微分公式(如两点公式、三点公式等)就是在等距步长情形下用插值多项式的导数作为近似值的。此外，还可以采用待定系数法建立各阶导数的数值微分公式，并且用外推技术来提高所求近似值的精确度。当函数可微性不太好时，利用样条插值进行数值微分要比多项式插值更适宜。如果离散点上的数据有不容忽视的随机误差，应该用曲线拟合代替函数插值，然后用拟合曲线的导数作为所求导数的近似值，这种做法可以起到减少随机误差的作用。数值微分公式还是微分方程数值解法的重要依据。

有限差分法最简单的方式是使用有限差分近似。1

简单的二点估计法是计算经过(x,f(x))及邻近点(x+h,f(x+h))二点形成割线的斜率选择一个小的数值h，表示x的小变化，可以是正值或是负值。其斜率为

此表示法是牛顿的差商，也称为一阶均差。

割线斜率和切线斜率有些差异，差异大约和h成正比。若h近似于0，则割线斜率近似于切线斜率。因此，函数f真正在x处真正的斜率是割线趋近切线时的差商：

若直接将h用0取代会得到除以零的结果，因此计算导数需要一些较不直觉的的方式。

同様的，切线斜率也可以用(x - h)和x二点的割线斜率近似。

另外一种二点估计法是用经过(x-h,f(x-h))和(x+h,f(x+h))二点的割线，其斜率为

上述公式称为对称差分，其一次项误差相消，因此割线斜率和切线斜率的差和成正比。对于很小的h而言这个值比单边近似还要准确。特别的是公式虽计算x点的斜率，但不会用到函数在x点的数值。

估计误差为：

其中c为在x-h和x+h之间的某一点。此误差没有包括因为有限准确度而产生的舍入误差。

很多工程计算机都是用对称差分来计算导数，像德州仪器（TI）的TI-82、TI-83、TI-84及TI-85，其h=0.001。

虽然在实务十分常用，但上述二种方式的数值微分常被研究者批评，尤其是被一些鼓励使用自动微分的研究者批评，因为上述的数值微分其精确度不高，若计算器精准度是六位数，用对称差分计算导数只有三位数的精确度，而若是找到一计算斜率的函数，仍可以有几乎六位数的精确度。例如假设f(x) = x，用2x计算斜率有几乎完整的准确度，而用差分近似就会有上述的问题。

利用浮点数在浮点数运算下，不同的 h造成的舍入误差及公式误差，只有在特定值下误差才是最小值

若计算时使用浮点数，就需要考虑h要取到多小。若选的太小，相减之后会有大的舍入误差，事实上整个有限差分的公式都是病态的，若h够小，导数不为零的情形下，在相消后会得到数值微分为零的结果，若h大小，计算割线斜率的结果就会更加准确，但用割线斜率估算切线斜率的误差就更大了。

一种可以产生够小的h，但又不会产生舍入误差的方式是（不过x不能为0），其中最小浮点数ε大约是2.2×10数量级。。以下是一个一个可以平衡舍入误差和公式误差，有最佳精确度的h为

（不过f"(x) = 0时不成立），而且需要有关函数的资讯。

上述的最小浮点数是针对双精度（64-bit）变数，单精度变数在这类计算几乎不太实用。其计算结果在二进制中不太可能是“整数”。虽然x是可以用浮点数表示的数字，但x+h几乎不会也是可用浮点数表示（而且和x不同）的数字，因此x+h需调整为机器可读的数字，因此会出现(x+h) -x不等于h的情形，因此用二个函数计算值计算微分时，二个位置的差不会是h。几乎所有的十进制分数在二进制下都会是循环小数（都像1/3在十进制中的情形一様），例如h= 0.1在二进制下会是循环小数，是 0.000110011001100...。因此在浮点数下一个可能计算的方式是：

h:=sqrt(eps)*x; xph:=x + h; dx:=xph - x; slope:=(F(xph) - F(x))/dx;

先计算(x+h) -x的值，再用这个值作为微分算式的分母，不过若是用电脑计算，编译器最优化的机能可能会认为dx和h相同，因此让上述的方式失效。若是用C或其他类似的编程语言，可以让xph宣告成Volatile变量，以避免此一问题。

高阶方法也有用更高阶估计导数的方法，或是估计高阶导数的方法。

以下就是一阶导数的五点法（一维下的五点模版）

其中

微分求积微分求积（Differential quadrature）是用函数在特定位置数值的加权和来近似导数，2其名称类似数值积分中用的求积（quadrature），也就是像梯形法或是辛普森法中用的加权和，有许多方式可找出加权的系数，在求解偏微分方程时会用到微分求积。

复变的方法传统用有限差分近似数值微分的方式是病态的，不过若f是全纯函数，在实轴上的值都是实数，可以用复平面中靠近 x的位置来求值，此方式为数值稳定的方式，例如一阶导数可以用以下的复数导数公式计算：

上述公式只在计一阶导数时有效，若要拓展到任意阶导数，需要用到多重复数，结果也会是多重复数的导数。

而任意阶的导数可以用柯西积分公式计算：

其中积分会用数值积分计算。

Lyness和Moler在1967年提出用复变数来计算数值微分。Abate和Dubner提出一种用复数拉普拉斯转换的数值反演为基础的算法。3