《ch11：数值计算方法之数值微分与外推方法》由会员分享，可在线阅读，更多相关《ch11：数值计算方法之数值微分与外推方法（39页珍藏版）》请在金锄头文库上搜索。
1、第11章： 数值微分与外推方法n设f(x)是x0-h, x0+h上连续可微的实函数，数值微分就是直 接利用计算f(x)的程序来计算f(x)在x0处的导数值。n对于工程应用来说，数值微分还是非常重要的：l如果函数是用列表的方式给出的，如果函数不是初等函 数，这时我们只能用数值微分方法求导数值。l为了争夺市场，现在的软件商更喜欢开发所谓傻瓜软件 ，这种软件不能要求用户给出一个函数的导函数。l直接采用数值微分方法可使相关的软件开发简单一些。n现在的计算工具也为我们求数值微分提供了极大的方便，主 要是计算时间和数值稳定性都得到明显改善，从而提升了的 数值微分的实用价值。11.1 利用差商外推加速n求数
2、值微分的方法还是比较多的，当我们强烈推荐的方法是 外推加速方法，基本步骤是：l构造一个与步长有关的近似公式；l构造一个步长收敛于零的变步长序列；l利用变步长序列外推，得到一个加速收敛序列。n利用外推方法加速的优点是：l尽管理论上有一定的高度，但学习起来并不困难，而编 程方法特别简单；l方法具有通用性，在后面求数值积分和常微分方程数值 解都用到了这种方法，也收到显著效果；l外推方法的效果特别好，在几乎不增加计算量的同时， 可以大幅度提升计算结果的精度。1.利用差商替代微商n计算公式假设f(x)是x0-h, x0+h上连续可微的实函数，当h的值充分小时，我们用f(x)在x0, x0+h处的一阶差商
3、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,即可简单地得到一个步长序列，相应地得到一个导数的近 似
4、值序列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可以取零，所 以是一种外推式的插值方法，亦即插值
5、点在所有插值基点所 在的最小区间的外边。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.111.2 利
6、用中心差商外推加速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 利用差商替代
7、微商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
8、偶函数在原点附近的多项式插值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中计算
9、zk和Ak的地方稍作修改，即可得 到利用中心差商外推的程序。l利用程序10.02计算前面的案例所得到的结果如下，可以 看出采用中心差商外推效果更好一些。图2.14.几点具体说明n对于求数值微分来说，对于给定的x0以及充分小的步长h， 计算差商时会遇到许多不利的运算：l计算x0+h和x0-h时，可能出现大数吃小数的现象；l计算f(x0+h)-f(x0-h)时, 可能出现绝对值相近符号相同的两 个数相减，从而产生较大的舍入误差。但是，采用变步长方法结合外推加速，基本上可以化解这个 问题。n尽管一般情况下利用中心差商计算效果更好一些，如果所给 的问题只能求单边导数，也只能用一阶差商外推。n在后面两节
10、中，我们将给出更一般性的外推方法，如果确定 了外推多项式的次数，那么在工程应用中就可以采用这两节 介绍的基本方法，程序显得更简洁一些。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,从 而又得到一个
11、加速收敛的序列。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)是一个偶函数，那么它本身就
12、具有 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处
13、泰勒展开，我们有：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 安排在下
14、面的计算表格中。 m=0m=1m=2m=3k=0T00k=1T01T11k=2T02T12T22k=3T03T13T23T3 3n n计算方法计算方法过去人们利用手工的方式计算，就是利用前面的过去人们利用手工的方式计算，就是利用前面的(*)(*)式，式， 在上面的表格中由上往下逐行计算，在每一行中，由左向在上面的表格中由上往下逐行计算，在每一行中，由左向 右逐列计算。右逐列计算。7.编写实用程序的考虑n现在利用计算机进行外推，为了强调程序的通用性，也为了 强调数值稳定性，建议只用34列计算，已经相当于采用68 阶精度的近似公式计算，进一步加列数已无实际意义。n采用固定列数计算时，可以不用指示列
15、的变量，从而可以用 4个数组T0,T1,T2,T3来存放表格各列中的数据。n为了得到精确结果，可适当增加计算的行数，所以实际的计 算表格像一个直角梯形。n矩形表格右上角单元格中的数据没有实际意义，但可以设置 得与同行左边相邻的数据相同，这样程序会特别简单。n即使增加计算的行数，数组的长度一般不要超过15，无论如 何不要超过31。所以对额外内存的需求可忽略不计。8.实用里查逊外推方法n结合计算机数值计算的特点，我们可以把实现里查逊外推方 法的程序设计得非常简单：只用一个简单的循环结构即可， 简化后的递推格式如下，完整的代码讲程序11.03，推广应用 时只需改写计算FT0(x,h)的代码即可。11
16、.4 涅维尔递推插值方法n大量科学计算问题都是利用迭代方法产生一对自变量和因变 量的收敛序列来获得其中一个序列的理论上的极限值。n数值计算面临的问题是一方面论分析存在截断误差，另一方 面实际计算又存在舍入误差，所以按最初的数学公式获得的 序列通常收敛较慢，甚至难以得到具有一定精度的近似值。n我们已经感受到采用插值方法利用一个已知的序列来产生一 个加速收敛序列，可以节省计算时间、提高计算精度。n埃特金方法和涅维尔方法都是利用利用已知的数值序列，通 过两个低阶多项式插值的结果的线性组合来产生一个较高阶 的多项式插值的结果，从而能自动调整插值精度。n埃特金方法和涅维尔方法本质上都是相同的，我们更看重的 是把涅维尔插值方法改造为通用的外推方法。1. 术语和记号n
