数值计算方法第六章

92第六章数值微分与数值积分本利用函数的插值理论，构造了相应的数值微分与数值积分计算公式，并通过Richardson（理查森）外推技巧提高了这些公式的计算精度；另外，本章还对所建立公式的收敛性及数值稳定性进行了分析；对于数值积分公式，还给出了精度的评价标准――代数精度，并由此出发，建立了精度最高的Gauss型求积公式．§6.1引言尽管微积分是现代科学的重要基础与起点，并且已在物理、力学、化学、生物等自然科学领域以及经济、管理等社会科学领域中有了非常广泛的应用，然而在微积分计算中，还至少面临着如下问题．(1)被积函数的原函数不易求出或者根本不能用初等函数表出．例如，概率积分202()(0)txptedxt和椭圆积分220()1sin(02)tetkxdxt(2)函数的表达式形式过于复杂，对其直接进行积分、求导运算，计算量很大；甚至对于一些形式简单的函数进行积分，得到的原函数也可能非常复杂．例如Cxtgarctgxdx233332cos2更重要的是，尽管上式从形式上看是精确的，然而，用上式计算定积分时，因函数值不能做到精确计算，最终得到的结果仍然是近似的，并且所引入的正切值和反正切值的计算量也不小．(3)函数本身就是离散函数，如实验数据等，用经典的微积分知识无法求其导数值和积分值．在计算机发展迅速的今天，上述困难可以用数值的方法予以解决．本章主要介绍常用的数值微积分公式及其相关理论．§6.2数值微分公式在微分学中，函数的导数是通过导数定义或求导法则求得的．然而如果需要利用函数在相关离散节点处的函数值近似计算其导数，就要寻求其它的方法．一种方法是利用离散数据进行插值，然后用插值函数的导数作为被插值函数导数的93近似；另一种方法是将不同点的函数值在某一点Taylor展开，然后用其线性组合建立函数的导数值近似表达；另外，还可以根据数值微分公式的截断误差，通过Richardson外推技巧建立更高精度的数值微分公式．一、插值法考虑以表格形式给出的，定义于区间],[ba上的离散函数：),,2,1,0()(nixfyii，用第四章的插值方法建立该函数的插值多项式)(xpn．由于多项式的导数容易求得，可以取)(xpn的导数作为函数)(xf导数的近似，即取)()(xpxfn．由(1)1()()()()()(1)!nnnffxpxxabn(6.1)得到)()!1()()()!1()()()()1(11)1(nnnnnfdxdnxxnfxpxf(6.2)上式中的是与x有关的未知函数，因此对于)()1(nfdxd很难做出估计．但是，在插值节点处的导数值)0()()!1()()()(1)1(nixnfxpxfinnini(6.3)进而得到函数在插值节点处的数值微分公式的截断误差)()!1()()()()(1)1(inniniinxnfxpxfxR(6.4)若函数)()1(nf在插值区间上有上界M，即有Mfn)()1(则得到数值微分公式的误差估计)()!1()()()(1ininiinxnMxpxfxR可以看出，当插值节点之间的距离较为接近时，通过插值方法得到的数值微分公94式有较高的精度．常用的数值微分公式是1n、2n时的插值型微分公式⑴一阶两点公式01()xx))()((1)()(0110xfxfhxfxf(6.5)100001111101()(),(,)2()(),(,)2hRxfxxhRxfxx(6.6)⑵一阶三点公式012()xxx210202121003421)(21)(4321)(xfxfxfhxfxfxfhxfxfxfxfhxf(6.7)2200002221110222222021()(,)31()(,)61()(,)3RxhfxxRxhfxxRxhfxx(6.8)利用类似的思路，还可以建立高阶导数的数值微分公式．从上面的数值微分公式的余项上来看，理论上步长h取的越小，精度就越高．但在实际计算时并不这样简单．例6.1设xexf)(，分别采用不同步长)15,2,1(100khhk，使用公式(6.5)计算)1(f的近似值．解显然8459052.71828182)1(ef．取初始步长1.00h，分别记)1()1(fhf，hd/，ederror，用Matlab软件编程计算，计算结果见表6.1．95从表6.1可以看出，当步长从1.0h到7101.0h逐步减小时，得到的导数近似值的误差逐步减少；然而，随着步长的进一步减小，误差逐步增大．这是因为实际计算结果的误差是由截断误差和舍入误差共同决定的，由截断误差表达式可知，h不宜过大；但当h小到表6.1例6.1计算结果kderror10.027319186557872.731918655787080.0136368273280320.002719641422532.719641422532780.0013595940737430.000271841774712.718417747078480.0001359186194440.000027182954202.718295419912310.0000135914532650.000002718283192.718283186986530.0000013585274860.000000271828202.718281963964840.0000001355058070.000000027182822.71828177744737-0.0000000510116780.000000002718282.71828159981169-0.0000002286473690.000000000271832.71827893527643-0.00000289318262100.000000000027182.71827005349223-0.00001177496681110.000000000002722.71827005349223-0.00001177496681120.000000000000272.71338507218388-0.00489675627516130.000000000000032.66453525910038-0.05374656935867140.000000000000002.66453525910038-0.053746569358671500-2.71828182845905一定程度时，公式的分子出现相近数相减，同时分母趋于零，舍入误差导致有效数字位数急剧减少，也同样导致求解精度的降低．二、Taylor展开法设函数)(xf充分光滑，对于给定的步长h，利用Taylor公式有32!3)(!2)()()()(hxfhxfhxfxfhxf(6.9)32!3)(!2)()()()(hxfhxfhxfxfhxf(6.10)从而可建立数值微分公式hxfhxfhOhxfhxfxf)()()()()()((6.11)其截断误差为)(hO．类似地，可以建立另一截断误差为)(hO数值微分公式96hhxfxfhOhhxfxfxf)()()()()()((6.12)如将式(6.9)和式(6.10)相减，可建立截断误差为)(2hO数值微分公式hhxfhxfhOhhxfhxfxf2)()()(2)()()(2(6.13)类似地，将式(6.9)和式(6.10)相加，可得截断误差为)(2hO的计算)(xf的数值微分公式222)()(2)()()()(2)()(hhxfxfhxfhOhhxfxfhxfxf(6.14)若利用更多点处的函数值(如),3(),2(hxhx)的Taylor展开式的线性组合，将可建立具有更高阶的截断误差的数值微分公式．三、Richardson外推法假设利用某种数值方法得到某一量S与步长h有关的近似值)(*hS，截断误差为kpkpphahahahSS2121*)((6.15)式中kppp210，系数,,21aa非零，且),2,1(,iapii均是与步长h无关的常数．用2/h代替上述公式中的步长h，得kpkpphahahahSS22222121*(6.16)如将上述两式进行加权平均，有望消去误差级数中的第一项，得到精度更高的数值计算公式．如在式(6.15)中，取*()()(),()2fxhfxhSfxShh并且根据(5)24()()()()()23!5!fxhfxhfxfxfxhhh(6.17)可得(5)*24()()()()()()23!5!fxhfxhfxfxSShfxhhh(6.18)97若将式(6.18)中的h用2/h替代，有24(5)*()()()()22()23!25!2hhfxfxhfxhfxhSSfxh(6.19)将以上两式线性组合，并消去2h的系数得)(2)()(31)2()2(34)(4hOhhxfhxfhhxfhxfxf)()(31)2(344**hOhShS(6.20)这是Richardson外推算法的第一步．当然若有必要，还可以对式(6.20)继续进行外推运算（将式中的)(4hO写成按h幂次展开的形式）．例6.2设xxxfcos)(2，试用Richardson外推法建立的数值公式(6.20)计算)0.1(f的近似值（外推3次），初始步长取1.0h．计算结果精确到小数点后7位．解计算结果见表6.2．表6.2例6.2计算结果h(6.12)式计算结果第一次外推第二次外推第三次外推0.10.226736160.050.236030920.239129170.0250.238357740.239133350.239133630.01250.238939640.239133610.239133630.23913363§6.3Newton-Cotes求积公式由高等数学的知识知，定积分badxxffI)()(定义为如下和式的极限)()(01iniiifxx式中011,[,](0,1,,1)niiiaxxxbxxin，极限过程为最大的子区98间长度趋于0．这个定义本身就给出了一种计算积分近似值的方法．数值积分实际上就是用一些点处函数值的线性组合来近似计算定积分，即)(][)()()(00iniiiniibaxfAfExfAdxxffI(6.21)式中niiA0称为求积系数，niix0称为求积节点，它们均与被积函数)(xf无关，][fE表示求积余项，也称为求积公式的截断误差．追求理想的求积公式自然是希望其精确程度最高．概括起来，数值积分需研究如下三个问题(1)求积公式的具体构造；(2)求积公式的精确程度衡量标准；(3)求积公式的误差估计；这里通过引入被积函数的插值函数以及求积公式代数精度的概念分析上述三个问题．一、插值法建立求积公式设给定的1n个互异节点为],[0baxnii，关于这些节点做被积函数的Lagrange插值函数)(xLn，于是有)()!1()()()(11xnfxLxfnnn(6.22)对式(6.22)两端在积分区间],[ba上积分，得)()!1()()()(11bannbanbad