23 数值微分

在没有函数微分的解析表达式时, 可以用数值方法由函数值近似计算函数的微分值。 按照微分定义, 当\(f(x)\)在点\(x\)处可微时, 对很小的\(h\)\[\begin{align} f'(x) =& \frac{f(x+h) - f(x)}{h} + O(h) \tag{23.1} \\ =& \frac{f(x) - f(x-h)}{h} + O(h) \tag{23.2} \\ =& \frac{f(x+h) - f(x-h)}{2h} + O(h^2) \tag{23.3} \end{align}\] 这三种近似分别称为数值微分的向前差商向后差商中心差商公式, 其中中心差商公式有更好的精度。 \(h\)选取应该很小, 使得继续减小\(h\)时按差商公式计算的近似导数值不再变化或变化量小于给定误差界限。 但是,\(h\)也不能过分小, 因为计算\(f(x)\)值只能到一定精度, 当\(h\)太小时\(f(x+h)\)\(f(x)\)的差别已经不是由\(h\)引起的而是由\(f(x)\)计算误差引起的。 如果\(f(x)\)的计算精度可以达到机器单位\(U\)量级, 取\(h\)的一个经验法则是取\(h/|x|\)\(U^{1/2}\)左右, 对中心差商公式取\(h/|x|\)\(U^{1/3}\)左右 (见 Monahan (2001) §8.6)。 多元函数的梯度可以类似地计算。

可以类似地计算高阶微分,如 \[\begin{align} f''(x) = \frac{[f(x+h) - f(x)] - [f(x) - f(x-h)]}{h^2} + O(h), \tag{23.4} \end{align}\] 公式分子有意地写成了两个差分的差分, 以避免有效位数损失。

公式(23.1)(23.3)是用割线斜率近似切线斜率, 相当于用两点的函数值作线性插值后用插值函数的微分近似原始函数微分。 一般地,对函数\(f(x)\)可以计算插值多项式\(P_n(X)\), 用\(P_n'(x)\)近似\(f'(x)\)。这种计算微分的方法称为插值型求导公式。 设多项式\(P_{n-1}(x)\)满足 \[\begin{aligned} P_{n-1}(x_i) = f(x_i), \ i=1,2,\dots,n \end{aligned}\] 由定理21.2, 插值近似的余项为 \[\begin{aligned} f(x) - P_{n-1}(x) = \frac{f^{(n)}(\xi)}{n!} g(x) \end{aligned}\] 其中\(g(x)=\prod_{i=1}^n (x - x_i)\), \(\xi\)依赖于\(x\)的位置。 插值型求导公式的余项为 \[\begin{align} f'(x) - P_{n-1}'(x) = \frac{f^{(n)}(\xi)}{n!} \cdot g'(x) + \frac{g(x)}{n!} \cdot \frac{d}{dx} \frac{f^{(n)}(\xi)}{n!} \tag{23.5} \end{align}\] 因为(23.5)中的第二项中\(\xi\)\(x\)的函数关系未知, 对任意\(x\)处用\(P_{n-1}'(x)\)近似\(f'(x)\)可能有很大误差。 如果仅在节点\(x_i, i=1,2,\dots,n\)处计算\(P_{n-1}'(x)\), 则有 \[\begin{align} f'(x_i) = P_n'(x_i) + \frac{f^{(n)}(\xi)}{n!} g'(x_i), \ i=1,2,\dots,n. \tag{23.6} \end{align}\]

根据(23.6), 我们导出对等距节点\((x_i, z_i)\)( \(i=1,2,\dots,n\))计算数值微分的公式。 设\(z_i = f(x_i)\)\(x_i - x_{i-1}=h\)

用线性插值近似求导的公式为 \[\begin{aligned} f'(x_i) = \frac{z_{i+1} - z_i}{h} - \frac{h}{2} f''(\xi_1), \ i=1,2,\dots,n-1 \ (\xi_1 \in (x_i, x_{i+1}) \end{aligned}\]\[\begin{aligned} f'(x_i) = \frac{z_i - z_{i-1}}{h} + \frac{h}{2} f''(\xi_2), \ i=2,3,\dots,n \ (\xi_2 \in (x_{i-1}, x_i) \end{aligned}\]

用二次多项式插值近似求导的公式为 \[\begin{aligned} f'(x_i) =& \frac{-3z_i + 4 z_{i+1} - z_{i+2}}{2h} + \frac{h^2}{3} f^{(3)}(\xi_1), \\ &\ i=1,2,\dots,n-2, \ \xi_1 \in (x_i, x_{i+2}) \nonumber \end{aligned}\]\[\begin{aligned} f'(x_i) =& \frac{-z_{i-1} + z_{i+1}}{2h} - \frac{h^2}{6} f^{(3)}(\xi_2), \\ & \ i=2,3,\dots,n-1, \ \xi_2 \in (x_{i-1}, x_{i+1}) \nonumber \end{aligned}\]\[\begin{aligned} f'(x_i) =& \frac{z_{i-2} - 4 z_{i-1} + 3z_i}{2h} + \frac{h^2}{3} f^{(3)}(\xi_3), \\ &\ i=3,4,\dots,n, \ \xi_3 \in (x_{i-2}, x_i). \end{aligned}\]

用二次多项式插值近似求二阶导数的公式为 \[\begin{aligned} f''(x_i) =& \frac{z_{i-1} - 2 z_i + z_{i+1}}{h^2} - \frac{h^2}{12} f^{(4)}(\xi_1), \\ & \ i=2,3,\dots, n-1, \ \xi_1 \in (x_{i-1}, x_{i+1}). \end{aligned}\]

用多项式插值方法近似计算函数微分, 在一般的\(x\)值处计算的误差可能很大。 可以采用样条函数来进行插值, 用样条函数的导数来计算被插值函数的导数。

习题

习题1

\(\boldsymbol x = (x_1, x_2, x_3, x_4)^T\), \(f(\boldsymbol x) = [\log(x_1 x_2 )/(x_3 x_4)]^2\), \(x_i>0, i=1,2,3,4\)。 求\(f(\boldsymbol x)\)的梯度\(\nabla f(\boldsymbol x)\)的表达式, 并编写R程序用数值微分方法估计\(\nabla f(\boldsymbol x)\), 比较两者的结果,并比较不同的数值微分方法的结果。

References

Monahan, John F. 2001. Numerical Methods of Statistics. Cambridge University Press.