数值分析知识内容 (35).pdf
6.6 数值微分 数值微分(Numerical Differentiation)就是用函数值的线性组合近似函数在某点的导数值.设函数)(xfy 在节点kx上的函数值为),2,1)(nkxfk，则我们可以用插值多项式的导数近似函数的导数，也可以按导数定义用差商近似导数，这些方法都可以得到数值微分的计算方法.6.6.1 差商求导方法 1、一阶求导公式 设)(xfy 在ix附近可以作泰勒（Taylor）展开，即 )(!3)(!2)()()(32iiiiixfhxfhxf hxfhxf.（7.43）由公式（7.43）解出)(ixf，得到 )(6)(2)()()(2iiiiixfhxfhhxfhxfxf)()()(iiixEhxfhxf，其中h为步长，)(ixE称为截断误差，若略去截断误差，得到一阶向前差商求导公式 hxfhxfxfiii)()()(.（7.44）以hxi代替ix，又可得到一阶向后差商求导公式 hhxfxfxfiii)()()(.（7.45）若将公式（7.44）与公式（7.45）相加，则有中点公式)(2)()()(hGhhxfhxfxfiii.（7.46）截断误差为 )(6)(2iixfhxE.【注】公式（7.44）、（7.45）、（7.46）均可用于计算一阶导数，尤其是中点公式更为常用.由导数定义知，当0h时，按公式去计算导数，精度会高些；但此时公式右端的分子是两相近的数相减，有可能造成有效数字的严重损失.因此需要选择适当的h.2、二阶求导公式 以hxi代替hxi，利用泰勒展开公式（7.43）得到 )(!3)(!2)()()(32iiiiixfhxfhxf hxfhxf.（7.47）将公式（7.47）与（7.43）相加，约去)(ixf 项，可得 )(12)()(2)()()4(22iiiiixfhhhxfxfhxfxf)()()(2)(2iiiixEhhxfxfhxf，其中，)(ixE为截断误差，若略去截断误差，得到二阶求导公式 2)()(2)()(hhxfxfhxfxfiiii.（7.48）例 11 利用中点公式计算xexf)(在0 x的一阶导数，取三位数字计算.解 利用一阶求导中点公式 heehGfhh2)()0(.取不同的步长h，计算结果见表 7-10（准确值为1)0(f）.表 7-10 例 10 计算值 h 1 0.5 0.1 0.05 0.01 0.005 0.001 0.0005 )(hG 1.18 1.04 1.03 0.99 1.00 1.50 0.50 0.00 由表 7-10 可知，01.0h的逼近效果最好，步长越增加（或越减小）的逼近效果都越差.因此，步长的选择不能太大也不能太小，最优步长的选择可根据误差分析得到.3、误差分析 对于一阶求导公式hhxfhxfhGii2)()()(.由截断误差 )(6)(3iixfhxE可得，MhxEhGxfii6)()()(2.（7.49）其中，)(maxxfMhxxi .显然，步长h越小，截断误差越小.再根据舍入误差，设)(hxfi和)(hxfi分别有舍入误差1和2，则计算)(ixf 的舍入误差上限为 hhhGxfxfii2)()()(21.（7.50）其中，,max21.显然，步长h越小，舍入误差)(ixf 越大，所以是病态的.由公式（7.49）和（7.50）可得计算)(ixf 的误差上限为hMhhE6)(2.要使)(hE最小，最优步长opth应为3/3Mhopt.同理，对于二阶求导)(12)()(2)()()4(22iiiiixfhhhxfxfhxfxf，由截断误差和舍入误差可得计算)(ixf 的误差上限为22412)(hMhhE.要使)(hE最小，最优步长opth应为4/48Mhopt.其中，)(max)4(xfMhxxi.6.6.2 插值型求导方法 插值型数值微分的基本思想：用插值多项式近似函数)(xf，从而用插值多项式在节点上的导数近似计算)(kxf.设函数)(xfy 定义在区间,ba上，在节点),2,1,0(nkxk处的函数值为ky，利用 Lagrange 插值方法，则有)()()(xRxLxfnn)()()!1()()(100)1(0nnknknkiiikixxxxxxnfyxxxx，从而)()()(xRxLxfnn.因此，得到插值型求导公式)()(xLxfn，)(xRn为数值求导的截断误差.当讨论节点处的导数时，求导公式和截断误差公式都比较简单，例如截断误差)()()!1()()(0)1(nnnxxxxnfxR )()()!1()()()()()!1(10)1(0)1(nnnnxxxxnfxxxxdxdfn.当kxx 时，上式右端第一项为零，因此 knnknxxxxxxnfxR)()()!1()()(0)1(.几个常用的等距节点的求导公式如下.（1）两点公式.设01xxh，对经过两点的插值多项式)(1xL求导，得到)()(1xLxf=1001yhxxyhxx=)(101yyh.因此，10,xxx 时均有)(1)()(0110yyhxfxf.（7.51）截断误差为 kkxxxkxxxxfxxxxfxR x01101 2)()(2)()(.当0k时，)(2)(01fhxR；当1k时，)(2)(11fhxR，其中，在0 x与1x之间.（2）三点公式.设hxxxx0112，二次插值多项式为 2120210121012002010212)()()()()()()(yxxxxxxxxyxxxxxxxxyxxxxxxxxxL，则有)(R)()(22xxLxf.（7.52）其中)(22)(0122022112202xxxxhyxxxxhyxxxxhyxL，)()()(!31)(210)3(2xxxxxxfdxdxR)()()(!3)(102021)3(xxxxxxxxxxxxf.当210,xxxx 时，分别得到)(322423)(22100fhhyhyhyxf ，（7.53）)(622)(2201fhhyhyxf ，（7.54）)(323242)(22102fhhyhyhyxf .（7.55）其中，公式（7.54）为一阶求导的中点公式，它只用两点函数值的组合，达到三点精度.由于对)(2xL求二阶导数，得到22212022)(hyhyhyxL，因此二阶求导的中点公式为)(12)2(1)()4(221021fhyyyhxf.（7.56）利用n次插值多项式)(xLn作为)(xf的近似函数，则可以得到高阶数值微分公式,2,1),()()()()()(kxRxLxfknknk，高阶数值微分公式需作误差分析，相应的数值计算精度更差.类似的可利用样条函数建立数值微分公式.6.6.3 数值微分的外推算法 利用中点公式计算一阶导数公式)()(21)()(hxfhxfhhGxf.对)(xf在点x做泰勒级数展开有 4221)()(hchchGxf，其中),2,1(ici与h无关.将h分半，4221)2/()2/()2/()(hchchGxf.上面二式消去2h项，则得 624114)()2/(4)(hdhdhGhGxf.其中),2,1(idi与h无关.若记)()(0hGhG，则有14)()2/(4)(001hGhGhG.同理，将h逐次分半，则有理查逊(Richardson)外推公式,2,1,14)()2/(4)(11mhGhGhGmmmm （7.57）截断误差为)()()()1(2mhOhGxf.由此可知，当m越大时，计算越精确.考虑到舍入误差，一般m不能取太大，在例 12 中，2m.例 12 设xexxf2)(，当h分别取025.0,05.0,1.0时，用中点公式求出5.0 x的一阶导数，进行外推，并与精确值454897994.0)5.0(f进行比较.解 先分别取025.0,05.0,1.0h，用中点公式)5.0()5.0(21)()5.0(2)5.0(20hhehehhhG 求出5.0 x的一阶导数值，见表 7-11.再由公式（7.57）外推两次，结果列于表 7-11 中.表 7-11 例 12 外推计算值 h)(0hG)(1hG)(2hG 0.1 0.05 0.025 0.4516049081 0.4540761693 0.4546926288 0.4548999231 0.4548981152 0.454897994 从表 7-11 可以看出，025.0h时只有 3 位有效数字，外推一次达到 5 位有效数字，外推两次达到 9 效有效数字，通常外推一次就能达到满意的结果.6.6.4 MATLAB 的微分函数 1、MATLAB 的符号微分 diff（）函数用于计算导数，其命令格式为 diff(S,v,n)表示求函数 S 关于变量 v 的 n 阶导数.例如，对于函数2sin)(xxf，求)(xf.利用 MATLAB 的符号微分 diff，输入 syms x,diff(sin(x2),x,1)或 syms x,diff(sin(x2)，得到 ans=2*cos(x2)*x.2、MATLAB 的数值微分 diff（）函数用于计算数值微分，其命令格式为 diff(S,n)表示求函数 S 的 n 阶差分.例如，对于函数2sin)(xxf，求)(xf.利用 MATLAB 的差分函数 diff，输入 h=0.001;x=0:h:pi;yy=diff(sin(x.2)/h;显然，取 h=.001 时,差分 diff(sin(x.2)/h 是导函数 2*cos(x.2).*x 的近似函数.
