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

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