Numerical Differentiation

1、Chapter 6 Numerical Differentiation数值微分数值导数的公式对于开发求解常微分方程（ODE）和偏微分方程（PDE）边值问题的算法很重要（参见第九章和第十章）。数值微分中经常采用已知函数为例，以便数值逼近的结果与准确值相比较。为了说明问题，考虑贝塞尔函数（the Bessel function），它的值列表可以在标准参考书中查到。闭区间上八个等距点和。它所基于的原理是插值多项式的微分。考虑寻求。插值多项式经过三个结点和，并由此得到。这个二次多项式和它在点处的切线如图所示。如果使用五个插值点，可以得到更好的逼近。经过点和的多项式，并由此得到。这个四次多项式和它在点处

2、的切线如图所示。这个导数的准确值是，在使用和时的误差分别为和。在这一章，我们研究数值微分精确度的相关理论。此处所举例子是可查值列表的贝塞尔函数，以便进行逼近值与真值的比较。采用的方法是用与函数相应的插值多项式的导数值来近似其导数值。分别采用了二次插值多项式和四次插值多项式来近似，并计算了相应误差，显然后者近似效果优于前者。事实上，这里提供了一种近似函数导数值的方法，即用与函数相应的插值多项式的导数值来近似其导数值。还要进行的思考是，这种方法的近似效果如何？怎么提高近似精度？再就是，有没有其他的近似方法？在前面的插值理论中已经讨论过，若函数的导数在定义闭区间上一致有界，且步长，则选择较大的将得到

3、较小的（即有误差界），从而高次逼近多项式将产生较小的误差。也讨论过高次逼近时在区间上发生的振荡（节点数增加，振荡更剧烈），即龙格现象，并提出了应对策略分段低次插值。所以说，通过插值多项式的导数来近似函数的导数值是切实可行的，并且由于多项式求导简单易行，这种方法就显得更加方便。考虑近似精度的话，相对而言，用高次插值多项式近似误差要比低次的小，当然这是在不致引起振荡（次数不太高）的前提下。毕竟插值多项式的导出有时是很复杂的，比如三次样条插值函数，需要数次求导积分并利用函数值、导数连续以及边界条件求待定系数来推倒得到，计算复杂，况且我们所需的只是逼近导数值，没有必要把插值多项式（函数的逼近）求出来，

4、所以本节探索仅用函数值的组合去寻求导数值的近似值，这是很自然的。6.1 Approximating The Derivative导数的近似值6.1.1 The Limit of the Difference Quotient差商的极限 考虑逼近函数的导数的数值过程：这个方法看起来很直接：选取序列，使得，计算下面序列的极限：读者可能注意到我们只需要计算序列中的有限项，将作为答案。这样的问题经常被提出：为什么计算？同样会问：如何取使得是的最佳逼近？为回答这个问题，我们再看一个例子，会发现为什么这不是件简单的事。 例如，考虑函数，分别使用步长和去构造点和点之间的割线。当变小时，割线逼近于切线，如图6

5、.2所示。尽管图6.2给出了处理过程的可视化表示，我们必须计算时的值去获得可以接受的数值解，而当取这个值时割线与切线的图像已不可分辨。上面处理为的过程，目的就是一个将连续问题离散化，把问题转化为计算机可以处理的离散数据输入。下面问题是计算不能无休止地进行，只要达到一定的精度要求计算就可以停下来，也就是文中提到的有限项。从下面的例子我们还将看到，这不仅仅是考虑计算必须停止所要求的，也是因为并非越小越好，太小甚至会扩大误差。例6.1 设且。使用步长，其中，计算差商。精度为小数点后九位。 在计算中所用到的和的值如表6.1所示。 最大值不能得到好的近似，因为步长太大，使得两点分割太远，差商是经过这两点

6、的割线的斜率，不能很好地近似切线。当以小数点后9位的精度计算公式时，根据得到近似值，根据得到近似值。如果太小，那么算得的与就非常接近。作差时就会导致精度损失，因为两个量几乎相等了。值太小了以至于存储的值与是相等的，因此计算的差商是0。在例6.1中，极限值是。可以看到，取值可以给出最佳逼近，。从例6.1看出，找到方程的数值极限并不容易。这个差商序列一开始向收敛，是最接近的，然后开始远离。程序6.1中表明，当时序列中的项停止计算。这是为了在项开始远离极限值之前，决定最佳逼近。将这个判定条件应用到例6.1中，有，因此是我们选择的答案。现在我们开始讨论如何根据较大的值得到合理精度的近似值公式。也就是说

7、，太大，割线不能很好地近似切线，但是也不能太小，否则由于精度损失，反而使误差增大。上面的例子已经显示出，求的极限的数值近似解时，选择适当的的重要性。在表格6.1中最后一列看到，在时最接近真实值，然后逐渐远离。其实，将相邻的两两作差，可以发现，这个差值是由大变小，又由小变大的，所以必定存在，使得，这个就是所求。在实际的程序操作中，便可把上面的不等式，作为迭代运算停止的判定条件。说完了存在性，那么为什么要以作为判定条件呢？这不禁让人想到极限收敛的Cauchy判别方法，它的要义是要求当自变量值在一定范围内时，函数值无限接近，也可以说函数值的差要任意小。我们在计算精度的限制下，不可能做到函数值之差任意

8、小，但起码可以保证找到函数值差别最小的，这就是上面的，这样这个就是最佳逼近所需的。6.1.2 The Central-difference Formulas中心差分公式如果函数在点的左边和右边的值可计算，则最佳二点公式包含两边的两个对称的横坐标。定理6.1（精度为的中心差分公式）设，且，则进一步，存在数，满足其中项称为截断误差（the truncation error）。证明 从关于的二阶泰勒展开式开始，写出与：用式减去式，结果是因为是连续的，根据中值定理可以找到值，使得将此式代入，得上式中右边的第一项是中心差值公式，第二项是截断误差。定理证毕。假设三阶导数的值变化不快，则式中的截断误差以与同

9、样的方式趋近于0，表示为。当用计算机进行计算是，不宜将选得太小。因此，求解近似值的公式具有精度为的截断误差项是有用的。定理6.2（精度为的中心差分公式）设，且，则而且存在数，满足其中。证明 一种导出公式式的方法如下。从关于四阶泰勒展开式开始，作和的差：然后，使用步长代替，写出以下近似：接下来用方程的8倍减去。包含的项就被消去了，得到如果的符号或是正或是负，而且它的值变化不快的话，在区间上可以找到一个值，使得把式代入，得到计算的结果式中右端的第一项是中心差分公式，第二项是截断误差。定理证毕。 假设在上有界，中的截断误差以与同样的方式趋近于0，用符号表示。现在来比较一下公式和。假设有五阶连续导数，

10、与大约相等。这样四阶公式的截断误差是的，比二阶公式中的截断误差更快地趋近于0，这样可以使用更大的步长。 回顾一下前面介绍的两个定理以及得到的差分公式。在定理6.1中，用式减去式，消去了。同样的，在定理6.2中，用方程的8倍减去。包含的项就被消去了。剩下的部分，根据中值定理并作必要的假设，合并整理出仅含的部分和误差项部分。这样得到的的表达式中就不包含导数信息，实现了仅用函数值来近似。定理6.2得到的四阶公式的截断误差是的，比定理6.1得到的二阶公式中的截断误差更快地趋近于0，这样就可以使用较大的步长值以得到合理精度下的近似值。下面的例子，就在不同的步长下，分别采用公式和近似导数值，很好地说明了这

11、一点。例6.2 令步长取，使用公式和，计算的近似值。在所有的计算中均采用小数点后九位的精度。与准确值比较。 取使用公式，得到取使用公式，得到由公式和公式得到的逼近值的误差分别是-0.000011941和0.000000017。此例中，当取时，公式比公式近似的效果要好。误差分析将解释这个例子并阐明为什么会这样。其他的计算结果在表6.2中列出。再来分析一下由公式和公式得到的数值微分的结果。首先，用式进行逼近时，在时，误差最小。而用式进行逼近时，在时，误差就可以达到最小。换言之，使用式逼近导数值时，可以取较大步长实现要求精度。其次，可以看到，用式得到的最佳近似值的误差即最小误差是-0.0000000

12、91，而利用式得到的最佳近似值的误差仅为0.000000017，所以说，即便是在不考虑步长的优势下，公式也比公式的近似效果要好。下面一个小节的工作就是将舍入误差考虑进去对误差进行分析，并通过考察总误差取最小以找到步长的最优值，讨论的只是纵向最优，没有进行公式和公式的横向比较。6.1.3 Error Analysis and Optimum Step Size误差分析和优化步长关于数值微分的一个重要课题是研究计算机的舍入误差。下面将对此进行更深入的分析。假设使用计算机进行数值计算，而且有和其中和是数值和的近似值，和分别是相关的舍入误差。下面的结论说明了数值微分中误差分析的复杂性。简单的进行一下符

13、号的说明，与中的的表示源于自变量，即，这里只是一种符号表示，对于往更高阶推广，使用更方便一些。推论6.1设函数满足定理6.1中的假设，并利用计算公式误差分析可通过如下的方程进行解释其中这里的总误差项（the total error term）是舍入误差（round-off error）与截断误差（truncation error）之和。推论6.1设函数满足定理6.1中的假设，且进行数值计算。如果，且，则式右边最小时的值为 在式右边关于求导，令导数为0，得到上式中的最优值。当较小时，式中包含的部分相对较大。在例6.2中，当时，就会碰到这种情况。舍入误差为截断误差为式中的误差项为实际上，当时的导数

14、近似值可用下式进行计算显然少了4位有效数字。误差是-0.000003909，接近预计的误差-0.000004001。在这一部分，使用例6.2，说明了上面两个推论的结果，并表明当时，舍入误差相对较大。上面和 对应的舍入误差和，数量级都在，的数量级是，以致总的舍入误差的数量级为，而截断误差项的数量级为，这样舍入误差相对于截断误差较大。而当时，和 对应的舍入误差和，数量级仍是，的数量级是，这样总的舍入误差与截断误差项的数量级均为，在同一数量级上！此时的最优，不只是减少舍入误差或者截断误差，实际上保证了舍入误差与截断误差和的最小，是一种博弈的结果。当将式用于例6.2时，可用界和值计算舍入误差。的优化值

15、为0.001144714。步长时，最接近于优化值0.001144714，而且在包含式的4个选项中，它给出了的最佳近似值（如表6.2和图6.3所示）。 对于式的误差分析跟上面的类似。设用计算机进行数值计算，并有函数。推论6.2设函数满足定理6.2的假设，使用计算公式误差分析可用下列方程进行解释： 其中这里的总误差项是舍入误差与截断误差的和。推论6.2设函数满足定理6.2的假设，且进行数值计算。如果且，则式最小时的值为当将式用于例6.2时，可用界和值计算舍入误差。的优化值为0.022388475，步长时最接近优化值0.022388475，而且在包含式的4个选项中，给出了的最佳近似值（如表6.2和图

16、6.4所示）。 这里可以将式和的结果再次横向比较。图6.3和6.4直观地展示了式不仅可以去较大的步长，而且精度确实比式要高，即优于。反映在图像上就是图线的最低点更低，离轴更近，而的标度可以取更大。下面再次提到了篇首的做法，用插值多项式的微分近似函数的微分。通过另一种推导也可得到数值微分公式，即对插值多项式进行微分。例如，经过点，和的2次多项式的拉格朗日表达式为将它展开可得通过类似的计算可得到经过点和4次多项式对这些多项式进行微分，可得和，与表6.2中下面的值相符。和以及它们在点的切线分别如图6.5和图6.5所示。6.1.4 Richardsons Extrapolation Richardso

17、n外推法 这一节将重点研究式和式之间的关系。设,且用和分别表示以和为步长根据式得到的的近似值，表示为和如果对式乘以4，并减去式，则可消去包含的项，结果为对式进一步分解，可得式中最后一个表达式就是中心差分公式。例6.3 设，并利用式和式，取，说明式中的线性组合如何用来求出式给出的的近似值。精度为小数点后9位。利用式和式，取，可得和式中的线性组合为这与例6.2中直接用式所得的的近似值相同。这种从低阶公式中推导出求解高阶导数的方法称为外推法。（这句话表述不够清楚，意思是将误差阶由提高到的方法，称为外推法）这里已经看到了如何使用步长和消去包含的项。为了说明如何消去，用和分别表示以和为步长根据式得到的精度为的的近似值。则近似值可表示为和设只为正或负值，而且变化不快，则可用假设来消去式和式中的，结果为下面的结论描述了提高计算精度的一般形式。定理6.3（Richardson外推法）设的两个精度为的近似值分别为和，而且它们满足：和这样可得到改进的近似值表达式 本节由用插值多项式的微分来近似函数微分引入，探讨了如何