数值计算方法之数值微分与外推方法.ppt

第11章：数值微分与外推方法,设f(x)是x0-h,x0+h上连续可微的实函数，数值微分就是直接利用计算f(x)的程序来计算f(x)在x0处的导数值。对于工程应用来说，数值微分还是非常重要的：如果函数是用列表的方式给出的，如果函数不是初等函数，这时我们只能用数值微分方法求导数值。为了争夺市场，现在的软件商更喜欢开发所谓傻瓜软件，这种软件不能要求用户给出一个函数的导函数。直接采用数值微分方法可使相关的软件开发简单一些。现在的计算工具也为我们求数值微分提供了极大的方便，主要是计算时间和数值稳定性都得到明显改善，从而提升了的数值微分的实用价值。,11.1 利用差商外推加速,求数值微分的方法还是比较多的，当我们强烈推荐的方法是外推加速方法，基本步骤是：构造一个与步长有关的近似公式；构造一个步长收敛于零的变步长序列；利用变步长序列外推，得到一个加速收敛序列。利用外推方法加速的优点是：尽管理论上有一定的高度，但学习起来并不困难，而编程方法特别简单；方法具有通用性，在后面求数值积分和常微分方程数值解都用到了这种方法，也收到显著效果；外推方法的效果特别好，在几乎不增加计算量的同时，可以大幅度提升计算结果的精度。,1.利用差商替代微商,计算公式假设f(x)是x0-h,x0+h上连续可微的实函数，当h的值充分小时，我们用f(x)在x0,x0+h处的一阶差商fx0,x0+h作为f(x0)的近似值，从而得到一个与步长h有关的近似公式,1.利用差商替代微商,估计截断误差记fx0,x0+h为f(x)在x0,x0+h处的差商，在上面的(1)式中，把f(x0+h)在x0处一阶泰勒展开，简单处理后后即可得到|f(x0)-fx0,x0+h|=O(h)定义：如果一个关于微小增量h的近似公式的截断误差与hk成正比，则称该公式具有k阶精度。结论 上面的(1)式具有1阶精度。,2.利用变步长方法提高精度,构造一个单调收敛于零的步长序列进行计算对于给定的初始步长h00，我们可以令 hk=hk-1/2,k=1,2,即可简单地得到一个步长序列，相应地得到一个导数的近似值序列 Zk=fx0,x0+hk,k=0,1,基本结论从理论上讲，当k时，我们有hk0，从而zk收敛于f(x0)。一般的实际情况是，首先当k增大时，截断误差会显著减小，到一定程度后，舍入误差又会显著增大，所以存在临界的k值。,2.利用变步长方法提高精度,确定停机规则 利用一个循环结构计算序列（hk,zk）|k=0,1,当然是一件很容易的事情，但确定停机规则需要慎重。根据实际精度要求来决定是否停机。事先确定充分小的正数EPS，只要某个|zk+1-zk|EPS即停止计算。如果EPS太小，有可能永远达不到精度要求。根据最好的结果停机。如果永远达不到精度要求，我们也可以力争得到最好的结果，也就是发现|zk+1-zk|开始随k的增大而增大时，停止进一步的计算。也可以综合上面两条规格而形成一个综合性的规则。,3利用外推方法加速收敛,问题的提法对于前面变步长方法得到的序列（hk,zk）|k=0,1,我们可以把他们看成是自变量为h的函数 z=fx0,x0+h 的函数值列表，从而可以从中取出若干对数据作插值多项式p(h)，进而得到当h=0时多项式的值。虽然可以用Lagrange插值公式完成相应的计算，由于多项式但自变量为零时的值就是常数项，所以能够找到更好的计算方法。由于我们是利用若干个hk0构造插值多项式来推算fx0,x0+h当h趋近于零时的极限值，而插值多项式中的h可以取零，所以是一种外推式的插值方法，亦即插值点在所有插值基点所在的最小区间的外边。,3利用外推方法加速收敛,以抛物线插值微利说明处理方法假设利用(hk,zk)、(hk+1,zk+1)、(hk+2,zk+2)作插值抛物线z=a0+a1h+a2h2，利用hk+1=hk/2,hk+2=hk/4,不难写出由此解得,3利用外推方法加速收敛,建立递推格式在上面的(2)式中，可以把a0记为Ak+2,从而得到递推格式,3利用外推方法加速收敛,编程与案例计算利用上面的递推格式编程是很容易的，源代码可参看教材第305页程序11.01。利用程序11.01求平方根函数在x0=4处的导数值得计算结果又下表给出，不难看到外推法的效果。,图1.1,11.2 利用中心差商外推加速,利用中心差商外推加速与利用差商外推加速的基本思路完全相同，也采用三步走的方法：构造计算导数的近似公式，其误差也是与步长h有关；构造变步长序列，从而得到导数的近似值序列；利用所得到的近似值序列构造插值多项式，从而得到多项式的常数项。与上一节相比，差别只是计算导数的近似公式不同，所以我们关注的重点是这种计算公式的差别所导致的方法上的细微差别以及性能上的差异。,1.利用中心差商替代微商,计算公式假设f(x)是x0-h,x0+h上连续可微的实函数，当h的值充分小时，我们用f(x)在x0,x0+h处的中心差商fx0-h,x0+h 作为f(x0)的近似值，即,1 利用差商替代微商,估计截断误差把f(x0+h)和f(x0-h)分别在x0处2阶泰勒展开，得由此不难得到所以，中心差商是一个2阶公式。,2.利用变步长方法提高精度,构造一个单调收敛于零的步长序列进行计算与前面的处理方式完全相同，对于给定的初始步长h00，我们可以利用地退的方式得到一个序列基本结论对于这里得到的序列，也具有与前面相同的结论。由于中心差商具有2阶精度，所以这里的序列会收敛的更快一些，所能得到的最好结果也会更精确一些。,3利用外推方法加速收敛,偶函数的性质不难验证，对于固定的x0来说，中心差商fx0-h,x0+h作为h的函数是偶函数。不难证明，偶函数在原点处的所有奇数次阶的导均为零。偶函数在原点附近的多项式插值如果在原点附近用插值多项式来近似表示一个偶函数，我们可以用仅含偶数次项的多项式作为插值多项式。仅含偶数次项的多项式的一般形式为 pm(x)=a0+a1x2+a2x4+amx2m(4)所以，如果是偶函数在原点附近的多项式插值，我们可以利用m+1对函数值列表得到2m阶精度的插值公式。,3利用外推方法加速收敛,具体计算方法假设利用(hk,zk)、(hk+1,zk+1)、(hk+2,zk+2)作3个基点的多项式插值，利用hk+1=hk/2,hk+2=hk/4,不难写出把所求得的a0记为Ak+2即可得到,3利用外推方法加速收敛,案例计算只需把程序10.01中计算zk和Ak的地方稍作修改，即可得到利用中心差商外推的程序。利用程序10.02计算前面的案例所得到的结果如下，可以看出采用中心差商外推效果更好一些。,图2.1,4.几点具体说明,对于求数值微分来说，对于给定的x0以及充分小的步长h，计算差商时会遇到许多不利的运算：计算x0+h和x0-h时，可能出现大数吃小数的现象；计算f(x0+h)-f(x0-h)时,可能出现绝对值相近符号相同的两个数相减，从而产生较大的舍入误差。但是，采用变步长方法结合外推加速，基本上可以化解这个问题。尽管一般情况下利用中心差商计算效果更好一些，如果所给的问题只能求单边导数，也只能用一阶差商外推。在后面两节中，我们将给出更一般性的外推方法，如果确定了外推多项式的次数，那么在工程应用中就可以采用这两节介绍的基本方法，程序显得更简洁一些。,11.3 里查逊外推加速法,科学计算领域内有一大类类似的问题可以表述为：我们要计算的某个与x0有关的真值T(x0)可表为某个与步长h有关的连续函数F(x0,h)当h趋近于零时的极限值。为了得到T(x0)的尽可能精确的近似值，我们通常取一个单调趋近于零的正数速序列hk,从而得到T(x0)的近似值序列zk=F(x0,hk)。接下来对所有不小于m的k值，寻找zk,zk-1,zk-m的一个线性函数或线性组合，以得到一个比zk更为精确的Ak,从而又得到一个加速收敛的序列。里查逊外推加速法就是求解这一类问题的一种通用的算法。提示：虽然教材的后面两节中又给出一种更为通用的外推加速方法，但仍有必要研究里查逊外推加速方法。,1.精度的定义与问题的提出,假设F(x0,h)在h=0处可展为泰勒级数 F(x0,h)=T(x0)+a1h+a2h2+aphp+如果a10,则称F(x0,h)具有1阶精度；如果a1=A2=ap-1=0,ap0,则称F(t,h)具有p阶精度；结论：显然，一个近似计算公式F(x0,h)的精度愈高愈好。问题：我们能否利用利用低阶精度的公式构造高阶精度的公式？,2中点公式的优越性,如果最初给出的F(x0,h)是一个偶函数，那么它本身就具有2阶精度。然而，平时按数学或物理学中给出的定义，F(x0,h)通常只有1阶精度。如果F(x0,h)具有1阶或奇数p阶精度，而且允许h取负值，那么我们可以简单地令 从而得到G(x0,h)是偶函数，而且具有2阶或p+1阶精度。约定：不失一般性，我们总是假定得到的第一个近似公式是关于h的偶函数，具有2阶精度，并记为FT0(x0,h)。,3利用步长折半的方法得到4阶公式,假设FT0(x0,h)适当光滑，在h=0处泰勒展开，我们有,4如法炮制6阶公式,上面得到的FT1(x0,h)是一个4阶公式。不难验证，它还是偶函数。如果足够光滑，我们可以如法炮制公式。还是把FT1(x0,h)在h=0处泰勒展开，我们有：,5.获得高阶近似公式的一般方法,实际上我们已经得到了获得高阶近似公式的一般方法：假如已经得到了一个2k阶的近似公式FTk-1(x0,h)，如果他还是足够光滑的，那么 就是一个具有2k+2阶精度的近似公式,6 数值计算时采用的公式,为了便于数值计算，我们采用如下记号：对于给定的初始步长h0，记 H0=h0，Hk=Hk-1/2，k=1,2,由于x0是固定值，所以在记号中不予考虑。对任意 k=0,1,2,，以及m=0,1,k，记 Tmk=FTm(x0,Hk)利用上面的记号，前面的(*)式可表示为,6.三角形的计算表格,把k作为行号，m作为列号，即可把前面定义的诸Tmk 安排在下面的计算表格中。,计算方法过去人们利用手工的方式计算，就是利用前面的(*)式，在上面的表格中由上往下逐行计算，在每一行中，由左向右逐列计算。,7.编写实用程序的考虑,现在利用计算机进行外推，为了强调程序的通用性，也为了强调数值稳定性，建议只用34列计算，已经相当于采用68阶精度的近似公式计算，进一步加列数已无实际意义。采用固定列数计算时，可以不用指示列的变量，从而可以用4个数组T0,T1,T2,T3来存放表格各列中的数据。为了得到精确结果，可适当增加计算的行数，所以实际的计算表格像一个直角梯形。矩形表格右上角单元格中的数据没有实际意义，但可以设置得与同行左边相邻的数据相同，这样程序会特别简单。即使增加计算的行数，数组的长度一般不要超过15，无论如何不要超过31。所以对额外内存的需求可忽略不计。,8.实用里查逊外推方法,结合计算机数值计算的特点，我们可以把实现里查逊外推方法的程序设计得非常简单：只用一个简单的循环结构即可，简化后的递推格式如下，完整的代码讲程序11.03，推广应用时只需改写计算FT0(x,h)的代码即可。,11.4 涅维尔递推插值方法,大量科学计算问题都是利用迭代方法产生一对自变量和因变量的收敛序列来获得其中一个序列的理论上的极限值。数值计算面临的问题是一方面论分析存在截断误差，另一方面实际计算又存在舍入误差，所以按最初的数学公式获得的序列通常收敛较慢，甚至难以得到具有一定精度的近似值。我们已经感受到采用插值方法利用一个已知的序列来产生一个加速收敛序列，可以节省计算时间、提高计算精度。埃特金方法和涅维尔方法都是利用利用已知的数值序列，通过两个低阶多项式插值的结果的线性组合来产生一个较高阶的多项式插值的结果，从而能自动调整插值精度。埃特金方法和涅维尔方法本质上都是相同的，我们更看重的是把涅维尔插值方法改造为通用的外推方法。,1.术语和记号,问题的提法假如我们已经有了一对相关变量的长度为n+1的观测值的有序的序列(x0,y0),(x1,y1),(xn,yn)，且满足多项式插值的基本假设：诸xk两两互不相同。对适当给定的x，对m=0,1,n，我们希望充分利用所得到的数据的所有长度为m的连续的子序列进行m次多项式插值，由此可以产生最好的插值结果。关键记号对给定的k=0,1,n，对所有m=0,1,k，记pk,m(x)为利用xk,xk-1,xk-m这m+1个基点所确定的m次插值多项式。对任意k，当m=0时，pk,0(x)yk；对任意k，当m=k时，pk,k(x)就是以x0,x1,xk这k+1个基点所确立的插值多项式。,2.涅维尔插值的递推格式,定理：对于前面给出的问题及记号，我们有 证明：不难验证，对m=0及任意k值，(*)式均成立。对k-mjk+1，由于 所以有pk+1,m+1(xj)=pk+1,m(xj)=pk,m(xj)=yj。显然，pk+1,m+1(xk+1)=pk+1,m(xk+1)=yk+1;pk+1,m+1(xk-m)=pk,m(xk+m)=yk-m;所以，pk+1,m+1(x)满足m+2个插值条件。,3.几点具体说明,实际数值计算时，应采用教材第313页(4.4)式计算，可节省一点计算时间，程序稍简单一点。无论是手工列表计算还是编程由计算机计算，均可仿照上一节中的里查逊外推法的实际计算方法。事实上，如果我们把涅维尔插值方法应用于外推，里查逊外推方法便成为这里的特例。在工程应用中，可以先用涅维尔插值方法计算，并根数值结果得出一个比较好的插值多项式的阶数，在实际应用中再采用具体的方法编程。这样一来，程序显得更简洁一些。实际使用里查逊外推方法也可采用这种模式。在后面的数值微分领域，寻找比较好的外推方法是很有意义的，涅维尔插值方法为此提供了基础方法。,4.利用涅维尔插值方法构造外推格式,在涅维尔递推插值公式中，如果诸xk是一个单调，不妨假设是单调减，收敛于x的序列，这时相应的计算问题便是一个外推问题。不失一般性，可在涅维尔递推插值公式(1)中令xk-x=hk,此时，hk是单调收敛于零的序列。再 把公式中的诸pk,m,(x)改写为pk,m,则有在进一步讨论上面的公式的用法之前，我们进一步寻找偶函数在原点处的外推格式。,偶函数在原点处的外推格式,与数值微分方法类似，大量数值计算问题均可等价地变换为下面的“步长”收敛于零的极限问题：而且，在一般情况下还可以假定F(x0,h)是关于h的偶函数。偶函数在原点处的泰勒展式只含偶次幂的项，简单的处理方法是把前面的(2)式中所有hk改写为(hk)2,从而有 其中，pk,0=F(x0,hk)。,6.里查逊外推格式的推导,为了说明如何利用前面给出的外推格式，首先我们利用上面的(3)式来推导里查逊外推格式。仍假定h0为初始步长，令hk=h0/2k,k=0,1,n，即可得到：,11.5 实用外推加速方法讨论,对于寻找某个二元函数F(x0,h)当h趋近于+0市的极限问题，上面一节中已经给出了3个外推格式，首先熟悉一下他们各自的特点是很有好处的：(2)式是最一般性的公式，适用于所有场合；如果F(x0,h)是关于h的偶函数，虽然也能使用(2)式，但使用(3)式的效果更好一些，包括计算量更少一些，数值性能更好一些，所能得到的最好结果会更精确一些。(4)式是(3)式的特例，按每次步长减半的规则确定步长序列，所以程序最为简单。大家也可以在(2)式中按每次步长减半的规则确定步长序列，也可以得到类似于(4)式得简单结果。(2)式和(3)式的意义在于，我们对外推方法的研究，可以专注于寻找更好的步长序列，或者构造更好的步长序列的策略。,利用等比序列外推加速,我们最容易想到的构造步长序列的方法是产生一个单调收敛于零的等比数列作为步长序列，这样处理的确有很多好处：方法简单，取0q1，对适当给定的h0，取hk=h0qk,k=0,1,，即可搞定。由于几何序列具有无记忆性，与初始步长h0的选取没有关系，在数值计算领域，这是难得的优点。具体做法：取q=1/r，则有hk=h0/rk，一般情况下r取1.22之间的数为宜，权衡各种细节因素，取r=1.6是比较理想的。把诸hk分别代入到(2)式和(3)式中，化简，即可得到具体的外推格式。,2.利用单调增的整数序列外推加速,设an是一个单调增的整数序列，对于给定的h00，记hk=h0/ak,则hk是单调减收敛于零的序列。把他们分别代入到(2)式和(3)式中，也可以得到相应的外推格式。利用整数序列外推加速的好处是方法简单，可选择的范围更大，所以更容易寻找具体的外推格式。最简单的整数序列可采用等差数列，由于首项、公差、以及系列的长度均能对算法性能产生影响，再加之与初始步长有密切关系，不宜作为一般性的方法使用。利用斐波拉奇数列外推也许是个好主意：构造简单，从1、2、3开始，以后每一项都是前面相邻两项之和，相邻两项之比大体为1.6，与初始步长的选取基本没有关系。,11.6 二阶导数的计算方法,对于求f(x)在x0处的二阶导数来说，也可以采用基本的处理方法：构造一个与步长有关的近似公式F(x0,h)，最好是关于步长的偶函数；确定适当的初始步长；构造变步长系列；利用变步长系列外推得到加速收敛序列。对于求f(x)在x0处的二阶导数来说，可以采用 不难验证，它是具有二阶精度的偶函数。一般地，直接采用理查逊外推方法(利用程序11.03重写函数FT0(double x,double h)即可。)也能得到非常精确的结果。建议大家以这个例子体验一下各种外推方法的性能以及编程时的微小差别。,