第六章数值微分与数值积分本利用函数的插值理论，构造了相应的数值微分与数值积分计算公式，并通过Richardson （理查森）外推技巧提高了这些公式的计算精度；另外，本章还对所建立公式的收敛性及数值稳定性进行了分析；对于数值积分公式，还给出了精度的评价标准――代数精度，并由此出发，建立了精度最高的Gauss型求积公式．§6.1 引言尽管微积分是现代科学的重要基础与起点，并且已在物理、力学、化学、生物等自然科学领域以及经济、管理等社会科学领域中有了非常广泛的应用，然而在微积分计算中，还至少面临着如下问题．(1) 被积函数的原函数不易求出或者根本不能用初等函数表出．例如，概率积分2()(0)t xp t e dx t-=≤<∞和椭圆积分()(02)te t tπ=≤≤⎰(2) 函数的表达式形式过于复杂，对其直接进行积分、求导运算，计算量很大；甚至对于一些形式简单的函数进行积分，得到的原函数也可能非常复杂．例如Cxtgarctgxdx+⎪⎪⎭⎫⎝⎛=+⎰233332cos2更重要的是，尽管上式从形式上看是精确的，然而，用上式计算定积分时，因函数值不能做到精确计算，最终得到的结果仍然是近似的，并且所引入的正切值和反正切值的计算量也不小．(3) 函数本身就是离散函数，如实验数据等，用经典的微积分知识无法求其导数值和积分值．在计算机发展迅速的今天，上述困难可以用数值的方法予以解决．本章主要介绍常用的数值微积分公式及其相关理论．§6.2 数值微分公式在微分学中，函数的导数是通过导数定义或求导法则求得的．然而如果需要利用函数在相关离散节点处的函数值近似计算其导数，就要寻求其它的方法．一种方法是利用离散数据进行插值，然后用插值函数的导数作为被插值函数导数的近似；另一种方法是将不同点的函数值在某一点Taylor展开，然后用其线性组合建立函数的导数值近似表达；另外，还可以根9293据数值微分公式的截断误差，通过Richardson 外推技巧建立更高精度的数值微分公式．一、插值法考虑以表格形式给出的，定义于区间],[b a 上的离散函数：),,2,1,0( )(n i x f y i i ==，用第四章的插值方法建立该函数的插值多项式)(x p n ．由于多项式的导数容易求得，可以取)(x p n 的导数作为函数)(x f 导数的近似，即取)()(x p x f n'≈'．由 (1)1()()()() ()(1)!n n n f f x p x x a b n ξωξ++=+<<+ (6.1)得到)()!1()()()!1()()()()1(11)1(ξωωξ++++++'++'='n n n n n f dxd n x x n f x p x f (6.2)上式中的ξ是与x 有关的未知函数，因此对于)()1(ξ+n f dxd 很难做出估计．但是，在插值节点处的导数值)0()()!1()()()(1)1(n i x n f x p x f i nn i n i ≤≤'++'='++ωξ (6.3) 进而得到函数在插值节点处的数值微分公式的截断误差)()!1()()()()(1)1(i nn i n i i n x n f x p x f x R ++'+='-'=ωξ (6.4) 若函数)()1(ξ+n f在插值区间上有上界M ，即有M fn ≤+)()1(ξ则得到数值微分公式的误差估计)()!1()()()(1i ni ni i n x n Mx p x f x R +'+≤'-'=ω 可以看出，当插值节点之间的距离较为接近时，通过插值方法得到的数值微分公式有较高的精度．常用的数值微分公式是1=n 、2=n 时的插值型微分公式⑴ 一阶两点公式01()x x <94))()((1)()(0110x f x f h x f x f -≈'=' (6.5)100001111101()(),(,)2()(),(,)2hR x f x x hR x f x x ξξξξ''=-∈''=∈ (6.6)⑵ 一阶三点公式012()x x x <<()()()()()()()()()()()⎪⎪⎪⎩⎪⎪⎪⎨⎧+-≈'-≈'-+-≈'210202*********)(21)(4321)(x f x f x f h x f x f x f h x f x f x f x f hx f (6.7) ()()()2200002221110222222021()(,)31()(,)61()(,)3R x h f x x R x h f x x R x h f x x ξξξξξξ⎧'''=∈⎪⎪⎪'''=-∈⎨⎪⎪'''=∈⎪⎩(6.8)利用类似的思路，还可以建立高阶导数的数值微分公式．从上面的数值微分公式的余项上来看，理论上步长h 取的越小，精度就越高．但在实际计算时并不这样简单．例6.1 设xe xf =)(，分别采用不同步长)15,2,1( 100 =⨯=-k h h k，使用公式(6.5)计算)1(f '的近似值．解 显然8459052.71828182)1(≈='e f ．取初始步长1.00=h ，分别记)1()1( f h f -+=∆，h d /∆=，e d error -=，用Matlab 软件编程计算，计算结果见表6.1．从表6.1可以看出，当步长从1.0=h 到7101.0-⨯=h 逐步减小时，得到的导数近似值的误差逐步减少；然而，随着步长的进一步减小，误差逐步增大．这是因为实际计算结果的误差是由截断误差和舍入误差共同决定的，由截断误差表达式可知，h 不宜过大；但当h 小到95表6.1 例6.1计算结果一定程度时，公式的分子出现相近数相减，同时分母趋于零，舍入误差导致有效数字位数急剧减少，也同样导致求解精度的降低．二、Taylor 展开法设函数)(x f 充分光滑，对于给定的步长h ，利用 Taylor 公式有+'''+''+'+=+32!3)(!2)()()()(h x f h x f h x f x f h x f (6.9) +'''-''+'-=-32!3)(!2)()()()(h x f h x f h x f x f h x f (6.10) 从而可建立数值微分公式hx f h x f h O h x f h x f x f )()()()()()(-+≈+-+=' (6.11)其截断误差为)(h O ．类似地，可以建立另一截断误差为)(h O 数值微分公式hh x f x f h O h h x f x f x f )()()()()()(--≈+--=' (6.12)如将式(6.9)和式(6.10)相减，可建立截断误差为)(2h O 数值微分公式96 hh x f h x f h O h h x f h x f x f 2)()()(2)()()(2--+≈+--+=' (6.13)类似地，将式(6.9)和式(6.10)相加，可得截断误差为)(2h O 的计算)(x f ''的数值微分公式222)()(2)()()()(2)()(hh x f x f h x f h O h h x f x f h x f x f -+-+≈+-+-+='' (6.14) 若利用更多点处的函数值(如 ),3(),2(h x h x ±±)的Taylor 展开式的线性组合，将可建立具有更高阶的截断误差的数值微分公式．三、Richardson 外推法假设利用某种数值方法得到某一量S 与步长h 有关的近似值)(*h S ，截断误差为+++=-k p k p p h a h a h a h S S 2121*)( (6.15)式中 <<<<<k p p p 210，系数 ,,21a a 非零，且),2,1(, =i a p i i 均是与步长h 无关的常数．用2/h 代替上述公式中的步长h ，得+⎪⎭⎫⎝⎛++⎪⎭⎫ ⎝⎛+⎪⎭⎫ ⎝⎛=⎪⎭⎫ ⎝⎛-kp k p p h a h a h a h S S 22222121*(6.16)如将上述两式进行加权平均，有望消去误差级数中的第一项，得到精度更高的数值计算公式．如在式(6.15)中，取*()()(),()2f x h f x h S f x S h h+--'==并且根据(5)24()()()()()23!5!f x h f x h f x f x f x h h h '''+--'=--- (6.17)可得(5)*24()()()()()()23!5!f x h f x h f x f x S S h f x h h h '''+--'-=-=--- (6.18)若将式(6.18)中的h 用2/h 替代，有24(5)*()()()()22()23!25!2h hf x f x h f x h f x h S S f x h +--'''⎛⎫⎛⎫⎛⎫'-=-=--- ⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭ (6.19)97将以上两式线性组合，并消去2h 的系数得)(2)()(31)2()2(34)(4h O hh x f h x f h h x f h x f x f +--+---+=' )()(31)2(344**h O h S h S +-= (6.20) 这是Richardson 外推算法的第一步．当然若有必要，还可以对式(6.20)继续进行外推运算（将式中的)(4h O 写成按h 幂次展开的形式）．例6.2 设x x x f cos )(2=，试用Richardson 外推法建立的数值公式(6.20)计算)0.1(f '的近似值（外推3次），初始步长取1.0=h ．计算结果精确到小数点后7位．解 计算结果见表6.2．表6.2 例6.2计算结果§6.3 Newton -Cotes 求积公式由高等数学的知识知，定积分⎰=badx x f f I )()(定义为如下和式的极限)()(01i ni i i f x xξ∑=+-式中011, [,](0,1,,1)n i i i a x x x b x x i n ξ+=<<<=∈=- ，极限过程为最大的子区间长度趋于0．这个定义本身就给出了一种计算积分近似值的方法．数值积分实际上就是用一些点处函数值的线性组合来近似计算定积分，即)(][)()()(0i ni i i n i i bax f A f E x f A dx x f f I ∑∑⎰==≈+== (6.21)式中{}ni i A 0=称为求积系数，{}ni i x 0=称为求积节点，它们均与被积函数)(x f 无关，][f E 表示98 求积余项，也称为求积公式的截断误差．追求理想的求积公式自然是希望其精确程度最高．概括起来，数值积分需研究如下三个问题 (1) 求积公式的具体构造；(2) 求积公式的精确程度衡量标准； (3) 求积公式的误差估计；这里通过引入被积函数的插值函数以及求积公式代数精度的概念分析上述三个问题．一、插值法建立求积公式设给定的1+n 个互异节点为{}],[0b a x ni i ⊂=，关于这些节点做被积函数的Lagrange 插值函数)(x L n ，于是有())()!1()()()(11x n f x L x f n n n ++++=ωξ (6.22)对式(6.22)两端在积分区间],[b a 上积分，得())()!1()()()( 11 ⎰⎰⎰++++=ban n ban badxx n f dx x L dx x f ωξ () )()!1()()( )( 11 0⎰⎰∑++=++=ban n bak nk kdx x n f dx x l x f ωξ (6.23)其中)(x l k 为关于节点{}ni i x 0=的Lagrange 插值基函数, 取⎰ban dx x L )(作为积分的近似值，构造出插值型求积公式)()(0k nk k bax f A dx x f ∑⎰=≈ (6.24)其中求积系数⎰=bak k dx x l A )( (6.25)同时，求得该求积公式的求积余项为()⎰+++=ban n dx x n f f E 11)()!1()(][ωξ (6.26)二、Newton -Cotes 求积公式99将节点等距分布时建立的插值求积公式称为Newton -Cotes 求积公式．它的思想最早出现于1676年Newton 写给Leibnitz （莱布尼茨） 的信中，后来Cotes （柯特斯） 对Newton 的方法进行了系统研究，于1711年给出了点数10≤的所有公式的权．下面列出几个最简单、最著名的Newton -Cotes 求积公式(1) 用区间中点做零次插值多项式，得到一点Newton -Cotes 求积公式，称为中点求积公式，即()():()2baa b f x dx b a f Mf +⎛⎫≈-= ⎪⎝⎭⎰(6.27) (2) 用区间两个端点做一次插值多项式，得到两点Newton -Cotes 求积公式，称为梯形求积公式，即()(()()):()2bab af x dx f a f b T f -≈+=⎰(6.28) (3) 用区间两个端点及中点做二次插值多项式，得到三点Newton -Cotes 求积公式，称为Simpson （辛普森）求积公式，即()()4():()62bab a a b f x dx f a f f b S f ⎛⎫-⎛+⎫⎛⎫≈++= ⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭⎰(6.29)(4) 将区间四等分得到的五个节点（包括区间两个端点）做四次插值多项式，得到五点Newton -Cotes 求积公式，称为Cotes 求积公式，即()337()3212327() :()90424baf x dx b a a b a b a b f a f f f f b C f ≈-⎡+++⎤⎛⎫⎛⎫⎛⎫++++= ⎪ ⎪ ⎪⎢⎥⎝⎭⎝⎭⎝⎭⎣⎦⎰(6.30)关于更详细的Newton -Cotes 求积公式的求积系数，可以查表6.3得到． 例6.3 试分别用一点、二点、三点、四点以及五点Newton -Cotes 公式计算积分dxx ⎰+10 11的近似值（计算结果取5位小数）． 解 利用中点求积公式，得66667.021111111 0 =⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡+⨯≈+⎰dx x利用梯形求积公式，得75000.0211211110 =⎥⎦⎤⎢⎣⎡+⨯≈+⎰dx x 利用Simpson 求积公式，得100 表6.3 Newton -Cotes 求积公式系数69444.02121114161111 0 =⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡++⨯+⨯≈+⎰dx x利用四点Newton -Cotes 求积公式，得69375.0213211331113181111 0 =⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡++⨯++⨯+⨯≈+⎰dx x 利用Cotes 求积公式，得69317.0274311322111124111327901111 0 =⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡++⨯++⨯++⨯+⨯≈+⎰dx x 事实上原积分的准确值为12 01ln 0.693151dx x =≈+⎰可见这五个公式的近似积分的精确程度不同，求积节点越多，得到积分近似值越精确．101三、代数精度由Weierstrass （维尔斯特拉斯）多项式逼近定理[9]可知，对任意给定的精度要求，闭区间上的连续函数都可以用多项式近似．一个很自然的做法就是用精确积分的被积多项式的次数做为求积公式精确程度的度量．定义6.1 若求积公式对于任意的不高于m 次的代数多项式做为被积函数时都精确成立，而对某个m ＋1次多项式被积函数不精确成立，则称该求积公式具有m 次代数精度．由(6.26)式可知，n ＋1个求积节点的插值型求积公式的代数精度至少为n ．容易证明，求积公式具有m 次代数精度的充要条件是它对于m x x x x f ,,,1)(2=都准确成立，而对于1)(+=m x x f 不准确成立．通过代数精度的定义容易验证，中点公式、Simpson 公式、Cotes 公式的代数精度分别为1阶、3阶、5阶，但它们分别用到的是具有1个、3个以及5个求积节点的Newton -Cotes 公式．从中可以看到这样一个规律：当插值节点的个数)1(+n 是奇数时的Newton -Cotes 求积公式的代数精度至少为1+n 阶，其证明见参考文献[6]．利用代数精度的定义也可以确定数值积分公式中的相关参数，目标是使求积公式的代数精度尽可能高．只需要依次取被积函数)(x f 为 ,,,12x x ，并通过使求积公式准确成立来建立适当数目的方程，进而求解未知参数．例6.4 试确定参数C B A ,,，使得求积公式)2()1()0()(2Cf Bf Af dx x f ++≈⎰的代数精度尽可能高．解 根据代数精度的定义，令2,,1)(x x x f =时，求积公式精确成立，得⎪⎩⎪⎨⎧=+=+=++3/84222C B C B C B A 解之得31,34,31===C B A取()3x x f =，则102 423=⎰dx x()()423134)2(103=⨯+=++Cf Bf Af 此时，该求积公式至少具有3次代数精度． 取()4x x f =，有53224=⎰dx x 但是()()32023134)2(104=⨯+=++Cf Bf Af 因此，当32,34,0===C B A 时，求积公式达到最高的代数精度，且代数精度的次数为3次． 例6.4实际上是积分区间取]2,0[时的Simpson 求积公式．可以看到，在固定求积节点的情况下，用待定系数法得到的求积系数和用插值法建立的求积公式的求积系数是相同的．这是因为，用待定系数法建立的形如式(6.21)的求积公式至少具有n 次代数精度，而至少具有n 次代数精度的求积公式一定是插值型的[6]．四、Newton -Cotes 求积公式的截断误差定理6.1 设函数)(x f 在积分区间],[b a 上具有连续的二阶导数，则中点求积公式的截断误差为3()[]() ()48M b a E f f a b ηη-''=<< (6.31)证明 由Taylor 展开公式2()()()2222!2a b a b a b f a b f x f f x x a b ξξ''++++⎛⎫⎛⎫⎛⎫⎛⎫'=+-+-<< ⎪ ⎪⎪ ⎪⎝⎭⎝⎭⎝⎭⎝⎭对上式两端在区间],[b a 上积分，有⎰⎰⎪⎭⎫ ⎝⎛+-''+⎪⎭⎫ ⎝⎛+-=b a badx b a x f b a f a b x dx x f 22!2)(2)()(ξ由于函数22)(⎪⎭⎫ ⎝⎛+-=b a x x g 在积分区间],[b a 内不变号，而函数)(ξf ''在],[b a 上连103续，由积分第二中值定理知，在(,)a b 内存在一点η，使得)(48)(2!2)(][32ηηf a b dx b a x f f E b a M ''-=⎪⎭⎫ ⎝⎛+-''=⎰ 定理6.2 设函数)(x f 在积分区间],[b a 上具有连续的二阶导数，则梯形求积公式的截断误差为3()[]()()12T b a E f f a b ηη-''=-<< (6.32)证明 梯形求积公式的余项为()[]()() ()2!bT af E f x a x b dx a b ξξ''=--<<⎰由于))(()(2a x a x x --=ω在积分区间],[b a 内不变号，而函数)(ξf ''在],[b a 上连续，由积分第二中值定理知，在(,)a b 内存在一点η，使得)(12)())((!2)(][3ηηf a b dx b x a x f f E b a T ''--=--''=⎰ 定理6.3 设函数)(x f 在区间],[b a 上有连续的四阶导数，则Simpson 公式的截断误差为5(4)()[]()()2880S b a E f f a b ηη-=-<< (6.33)证明 对区间],[b a 上的函数)(x f ，构造次数不高于3次的插值多项式)(3x H ,使得)()( ),()(33b f b H a f a H ==22 ,2233⎪⎭⎫⎝⎛+'=⎪⎭⎫ ⎝⎛+'⎪⎭⎫ ⎝⎛+=⎪⎭⎫ ⎝⎛+b a f b a H b a f b a H不难得到()b a b x b a x a x f x H x f <<-⎪⎭⎫ ⎝⎛+--+=ξξ )(2)(!4)()()(243对上式两端在区间],[b a 上积分104 ()⎰⎰⎰-⎪⎭⎫ ⎝⎛+--+=bab abadx b x b a x a x f dx x H dx x f 24 3 )(2)(!4)()()(ξ 因Simpson 公式的代数精确度是3，所以()][)(24)(6 )(2)(!4)( )(24)(6)( 24333 f E b f b a f a f a b dx b x b a x a x f b H b a H a H a b dx x f S b aba+⎭⎬⎫⎩⎨⎧⎪⎪⎭⎫ ⎝⎛+⎪⎭⎫ ⎝⎛++-=-⎪⎭⎫ ⎝⎛+--+⎭⎬⎫⎩⎨⎧⎪⎪⎭⎫ ⎝⎛+⎪⎭⎫ ⎝⎛++-=⎰⎰ξ式中()()())(2880)()(2)(!4)()(2)(!4)(][4524 24ηηξf a b dx b x b a x a x fdx b x b a x a x f f E b a baS --=-⎪⎭⎫ ⎝⎛+--=-⎪⎭⎫ ⎝⎛+--=⎰⎰定理6.3所使用的证明方法也可以用来证明定理6.1．对于更为复杂的Cotes 公式，有 定理6.4 设函数)(x f 在区间],[b a 上有连续的6阶导数，则Cotes 公式的截断误差为()6642()[]() ()9454b a b a E f f a b ηη--⎛⎫=-<< ⎪⎝⎭(6.34)五、Newton -Cotes 求积公式的稳定性和收敛性本节简要分析Newton –Cotes 公式序列的收敛性及稳定性．定义6.2 Newton -Cotes 求积公式的收敛性是指，对于任意连续的被积函数，当求积节点的个数∞→n 时，截断误差序列满足0][lim =∞→f E n n ．Newton -Cotes 求积公式是基于等距插值建立起来的，但因当∞→n 时，插值函数不保证收敛到被插值函数，进而也不能保证0][lim =∞→f E n n ．Newton -Cotes 求积公式序列的稳定性主要考虑函数值)(i x f 计算的舍入误差对数值积分结果的影响．定义6.3 对于一个求积公式序列{}0()n n I f +∞=，假设在计算)(i x f 时，引入舍入误差i e ，即有i i i e x f x f +=)()(~，并记105)(~)~(),()(0i ni i n i ni i n x f A f I x f A f I ∑∑====如果对任意给定的小正数0>ε，只要误差i ie max =δ充分小，就有ε≤-=-∑=)](~)([)~()(0i i ni i n n x f x f A f I f I则称该求积公式序列{}0()n n I f +∞=是稳定的． 对于数值求积公式，有δ∑∑∑∑====≤≤=-ni in i ii in i i i i ni i A e A e A x f x f A 0)](~)([由于Newton -Cotes 求积公式对常数均能精确积分，当被积函数1)(≡x f 时，有在∞→n 的过程中，Newton -Cotes 求积公式序列的求积系数有正有负，nii A=∑的有界性不能保证，因而Newton -Cotes 求积公式得稳定性不能保证．基于上述稳定性、收敛性原因，在实际计算中，人们并不追求高阶的Newton -Cotes 求积公式，而是通过细化积分区间的方法，用复化求积公式来提高数值积分的精度．§6.4 复化求积法从积分的定义可以看出，增加求积节点数目可以提高数值积分的代数精度．但通过增加插值节点建立的高阶Newton -Cotes 求积公式的稳定性差．一种提高求积精度的方法是先将整个积分区间细分，然后在每个小的区间上使用低阶的Newton -Cotes 求积公式，这就是复化求积法．理论和数值结果都表明，这种方案可以获得理想的数值结果．一、复化梯形公式首先将积分区间],[b a n 等分，区间长度nab h -=，节点记为 (0,1,,)i x a ih i n =+= ，然后在每个小区间],[1+i i x x 上使用梯形公式，得106 ∑∑∑∑⎰⎰-=-=+-=-=''-⎪⎭⎫ ⎝⎛++=''-+==+131131101)(12)()(2)(2)](12))()((2[)()(1n i i n i i i i i n i n i x x baf h b f x f a f h f h x f x f h dxx f dx x f i iηη略去上式最后一项，得到复化梯形公式⎪⎭⎫⎝⎛++=∑-=11)()(2)(2n i i n b f x f a f h T (6.35)其截断误差为3110() ()12n n T i i i i i h E f x x ηη-+=''=-<<∑定理6.5 设函数)(x f 在积分区间],[b a 上具有连续的二阶导数，则复化梯形公式的截断误差为2() ()12n T b a E h f a b ηη-''=-≤≤ (6.36)证明 因为)(x f ''在],[b a 上连续，则其在区间一定有最大值M 和最小值m ，即1(),(,)[,]i i i i m f M x x a b ηη+''≤≤∈⊂进而有M f n m n i i ≤''≤∑-=1)(1η由连续函数的介值定理知，存在一点[]b a ,∈η，使得∑-=''=''1)(1)(n i i f n f ηη故)(12)(122103ηηf h a b f h E n i i T n ''--=''-=∑-=从截断误差可以看出，复化梯形公式的代数精度仍为一阶．107二、复化Simpson 公式若在每个小区间上使用Simpson 公式，有()()∑∑∑∑∑∑⎰⎰-=-=-=+-=++-=-=-⎪⎪⎭⎫ ⎝⎛+⎪⎪⎭⎫ ⎝⎛++=-⎪⎪⎭⎫ ⎝⎛+⎪⎪⎭⎫ ⎝⎛+==+10451110211045121101)(2880)(4)(2)(6)(2880)(4)(6)()(1n i i n i n i i i n i i i i i n i n i x x b af h b f x f x f a f h f h x f x f x f h dxx f dx x f i iηη上式略去最后一项，得到复化Simpson 公式⎪⎪⎭⎫ ⎝⎛+⎪⎪⎭⎫ ⎝⎛++=∑∑-=-=+111021)(4)(2)(6n i n i i i n b f x f x f a f h S (6.37)其截断误差为()51410() ()2880n n S i i i i i h E f x x ηη-+==-<<∑定理6.6 设函数()x f 在区间[]b a ,上有连续的四阶导数，则复化Simpson 公式的截断误差为()44() ()2880n S i b a E h f a b ηη-=-≤≤ (6.38)从截断误差可以看出，复化Simpson 公式的代数精度仍为三阶．对于上节的中点求积公式，使用同样的方法，也可以获得其对应的复化求积公式以及截断误差11022() ()48n n n i i M M h f x b a R h f a b ηη-+=⎛⎫= ⎪⎝⎭-''=≤≤∑可以证明这三个复化求积公式都是收敛的、稳定的[6]． 例6.5 使用复化梯形公式和复化Simpson 公式计算定积分⎰2/ 0cos πxdx 的近似值时，要求误差不超过510-，分别需要多少个求积节点？108 解 这里，(), , max 122a b b a b a h f n nηππη≤≤-''-====．从复化梯形公式的截断误差可知，要使得近似值的误差不超过510-，则需满足()521012max -≤≤≤''--ηηf h a b ba 解不等式得179.7170≥n取180=n ，相应的需要1811=+n 个求积节点；同理，用复化Simpson 公式计算时，需满足()5)4(4102880max -≤≤≤--ηηf h a b ba 由于()1max )4(=≤≤ηηf b a ，解541022880-≤⎪⎭⎫ ⎝⎛--n a b π得， 4.2688≥n ，因此，需要将区间进行5=n 等分，相应的需要1112=+n 个求积节点．三、区间逐次分半求积法使用复化求积公式计算积分值，要保证满足指定精度，依据截断误差，需要将中值定理使用时产生的被积函数的高阶导数在某一点η的绝对值放大为整个积分区间上的最大值，若最大值难求，还需进一步进行放大运算，势必导致过多地选取求积节点，影响计算效率．区间逐次分半求积法则避免了这种不便．它是根据规定的精度要求，在计算过程中把积分区间逐次分半，利用相邻两次求积结果之差来判别误差的大小，从而得到满足精度要求的积分近似值．下面以复化梯形公式为例来说明区间逐次分半求积法的计算过程．由复化梯形公式的误差估计式可以看到，)(2h O E n T ≈．当对n 较大时，可以认为2ch E n T ≈ (6.39)这里c 为常数，当将区间n 2等分时有222⎪⎭⎫⎝⎛≈h c E nT (6.40) 进而42≈--nnT I T I (6.41)即109)(3122n n n T T T I -+≈ (6.42)上式说明，若用n T 2作为准确值I 的近似值时，误差大约为)(312n n T T -．因此可以使用ε<-)(312n n T T 判断近似值n T 2的近似程度，其中ε为容许误差，这样就避免了被积函数高阶导数最值的估计．需要指出的是，在实际计算n T 2时，由于每次总是在前一次对分的基础上将区间再次对分，分点加密一倍，原分点上的函数值不需要重复计算，这样可以采用如下的递推公式减少计算量n a b h h i a f h T T n i n n -=⎪⎪⎭⎫⎝⎛⎪⎭⎫ ⎝⎛-++=∑=,212212 (6.43) 同理，也可以类似的给出基于复化Simpson 公式和复化Cotes 公式的区间逐次分半求积法，此时有)(15122n n n S S S I -+≈ (6.44) )(63122n n nC C C I -+≈ (6.45) 算法6.1 （复化梯形求积算法） 输入 区间端点b a ,及精度ε； 输出 积分近似值T ；Step 1 令1=n ，2/)(a b h -=，))()((0b f a f h T +=，/2h h =；Step 2 令0=F ，对n i ,,2,1 =，计算))12((h i a f F F -++=； //计算新节点函数值和 Step 3 计算hF T T +=2/0； //计算式(6.43)Step 4 若ε≤-0T T ，算法终止，输出T ；否则，赋值T T h h n n ===0,2/,2，转Step2． 算法6.2 （复化simpson 求积算法） 输入 区间端点b a ,及精度ε； 输出 积分近似值S ；110 Step 1 令2=n ，()4/a b h -=；Step 2 计算01()(),(()/2),F f a f b F f a b =+=+()6/)4(100F F a b S +-=，/2h h =； Step 3 令02=F ，对n i ,,2,1 =，计算()h i a f F F )12(22-++=；//计算新节点函数值和 Step 4 计算3/)42(210F F F h S ++=；Step 5 若ε≤-0S S ，算法终止，输出S ；否则，令S S h h n n ===0,2/,2，322F F F +=转Step3．例6.6 使用复化梯形公式和复化Simpson 公式计算定积分⎰10 sin dx x x的近似值，要求误差不超过510-=ε．解 使用区间逐次分半求积法的复化梯形公式和复化Simpson 公式进行求解．计算结果如表6.4所示，其中等分区间间距为：k a b h 2/)(-=，算例使用 ε<-n n T T 2和ε<-n n S S 2作为误差终止判断条件．表6.4 例6.6的计算结果§6.5 Romberg 求积法一、Romberg 求积法区间逐次分半求积法将)(312n n T T -或)(1512n n S S -处理成误差，将n T 2或n S 2作为积分111的近似值．分析近似表达式(6.45)和(6.47)后可以发现，既然)(n n S T 和)(22n n S T 都已经算出，自然地，误差)(312n n T T -或)(1512n n S S -也可以算出，那么将)(3122n n n T T T -+或)(15122n n n S S S -+整体作为积分的近似值，显然可进一步改善积分计算精度． 事实上，可以验证n n n n S T T T =-+)(3122 (6.46)上式表明对基于区间逐次分半求积的复化梯形公式进行误差校正得到复化Simpson 求积公式的计算值，从而把误差从)(2h O 提高到)(4h O ．类似地，有n n n n C S S S =-+)(15122(6.47) 即对基于区间逐次分半求积的复化Simpson 公式进行误差校正得到复化Cotes 求积公式的计算值，误差从)(4h O 提高到)(6h O ．我们称基于区间逐次分半求积的复化Cotes 公式进行误差校正得到的公式为Romberg （龙贝格）求积公式n n n n R C C C =-+)(63122 (6.48)它的误差阶是)(8h O ．二、Richardson 外推法Romberg 求积算法还可以通过Richardson 外推技巧得到．文献[4]给出了复化梯形公式的渐近展开形式的误差表达246246 ()bn af x dx T a h a h a h =++++⎰(6.49)其中， 642,,a a a 是与步长h 无关的常数．使用基于复化梯形公式的区间逐次分半求积法，可得梯形序列{}1n n T +∞=． 现将Richardson 外推技巧用到式(6.49)，经过一次外推得到112 (1)4(1)6246 41()33bn n af x dx T T a h a h =-+++⎰(6.50)而由式(6.46)知24133n n n S T T =-，于是有(1)4(1)646()bn af x dx S a h a h =+++⎰(6.51)这与定理6.6一致．类似地使用Richardson 外推技巧，还可以依次得到式(6.47)和式(6.48) ． 算法6.3（Romberg 算法）输入 积分区间端点参数b a ,及精度ε； 输出 积分近似值；Step 1计算初值()()1()()/2T b a f a f b =-+； Step2 令(1,2,;2)i b ah i n n-=== ，按式 (6.43) 计算梯形序列{}n T ； Step 3 分别按式 (6.46) 、式 (6.47)、 式 (6.48) 计算序列{}n S 、{}n C 、{}n R ； Step 4若(2)()R n R n ε-≤，算法终止，输出(2)R n ；否则，区间分半，转Step2．例6.7 使用Romberg 求积法计算定积分⎰+10 214dx x 的近似值，要求误差不超过710-=ε. 解 计算结果如表 6.5所示，其中等分区间间距为：n a b h 2/)(-=，算例使用ε<-n n C R 作为误差终止判断条件．表6.5 例6.7的计算结果容易计算358973.141592651410 2==+⎰πdx x可以看到，当4=k 时，使用基于复化梯形公式的区间分半求积法，误差是1130006510.03.1409416=-π但是通过外推之后，得到的龙贝格序列的误差减小为00000001.03.14159264=-π计算精度得到了很大提高．§6.6 Gauss 型求积公式前面研究的数值求积公式都是在积分区间内事先确定了1+n 个求积节点，然后通过待定系数法，按照使求积公式的代数精度达到最高的原则选取最佳求积系数（和通过插值的方法得到的求积系数是一致的）．已经知道，对于具有1+n 个求积节点的求积公式，其代数精度至少可以达到n 次．那么能否适当地选择求积节点的位置和相应的求积系数，使得求积公式具有更高甚至最高的代数精度？答案是肯定的．例如下面的求积公式()⎪⎪⎭⎫⎝⎛+⎪⎪⎭⎫ ⎝⎛-≈⎰-33331 1 f f dx x f 可以验证，该公式的代数精确度是3次，而前述的两点Newton -Cotes 公式即梯形求积公式的代数精度只能达到1次．实际上，对任意一组1+n 个互异求积节点n i i x 0}{=的插值型求积公式，当被积函数)(x f 取为22+n 次的多项式()()()2212021)(nn x x x x x x x ---=+ ω 时，总有0][][)()()(21021>=+==+=+∑⎰f E f E x A dx x f I i n ni i ban ωω(6.52)即对任意的1+n 个互异求积节点构造的求积公式的代数精度不能达到22+n 次，那么是否能够达到12+n 次，则是本节考虑的问题．一、Gauss 型求积公式一般地，如果选取1+n 个求积节点，要使求积公式对任意的12+n 次多项式精确成立，根据待定系数法，有114 01 22 0011 2222212121210011 1222b n a b n n a n n b n n n n n n a A A A dx b ab a A x A x A x xdx ba A x A x A x x dx n ++++++⎧++==-⎪⎪-⎪++==⎪⎨⎪⎪-⎪++==⎪+⎩⎰⎰⎰ (6.53) 理论上可以证明，上述非线性方程组的解是存在惟一的．也就是说，通过合理地选取求积节点及求积系数，可以使得求积公式的代数精度达到12+n 次．将这种具有1+n 个求积节点的具有12+n 次代数精度的求积公式称为Gauss 型求积公式，对应的一组求积节点称为一组Gauss 点．Gauss 型求积公式的求积节点和求积系数可以通过求解非线性方程组(6.53)得到，然而，直接求解该非线性方程组是复杂的．通常构造Gauss 型求积公式的方法是，首先通过正交多项式确定一组Gauss 点，然后再使用待定系数法或插值的方法求出Gauss 型求积公式的求积系数．定理6.7 对于插值型求积公式，其互异节点n i i x 0}{=是一组Gauss 点的充分必要条件是以这些点为零点的多项式函数101()()()()n n x x x x x x x ω+=--- 与任意的次数不超过n 次的多项式函数)(x p 在积分区间[]b a ,上正交，即满足0)()(1=⎰+ban dx x p x ω (6.54)证明 (充分性） 对任意不超过2n +1次的多项式)(x g ，存在不超过n 次的多项式)(x q 和)(x r ，使得)()()()(1x r x x q x g n +=+ω (6.55)成立．对上式两端积分，得⎰⎰⎰⎰=+=+bab ab an badx x r dx x r dx x x q dx x g )()()()()(1ω (6.56)对互异节点组n i i x 0}{=，构造插值型求积公式，即取求积系数115),,1,0(0n k dx x xx x A b an kj j jkjk =--=⎰∏≠=(6.57)这样求积公式至少具有n 次代数精度． 利用代数精度概念有∑∑⎰====ni i i n i i i b ax g A x r A dx x r 0)()()((6.58)比较式(6.56)和式(6.58)，有∑⎰==ni i i b ax g A dx x g 0)()((6.59)这样求积公式具有2n +1次代数精度，积分节点是一组Gauss 点．(必要性) 节点n i i x 0}{=是一组Gauss 点，n i i A 0}{=是相应的Gauss 求积系数，这样的求积公式具有2n +1次代数精度．设)(x p 是任意不超过n 次的多项式，)()(1x x p n +ω不超过2n +1次．0)()()()(011==∑⎰=++ni i n i i b an x x p A dx x x p ωω(6.60)这说明式(6.54)成立．由定理6.7可知，只要能够找到与任意的次数不超过n 次的多项式函数)(x p 在区间],[b a 上正交的n +1次多项式函数，那么，这个函数的所有零点即为Gauss 点．下面以Legendre (勒让德) 多项式为例，介绍Gauss 型求积公式的构造过程．二、Gauss-Legendre 求积公式称形如∑⎰=-≈ni i i x f A dx x f 01 1)()((6.61)的Gauss 型求积公式为Gauss-Legendre 求积公式． 可以证明，该求积公式对应的n +1个Gauss 点正好是n +1次勒让德多项式121111)1(2)!1(1)(+++++-+=n n n n n x dxd n x p (6.62)116 的一组零点．在求得一组Gauss 求积节点后，求积系数可通过待定系数法或插值型求积公式的求积系数公式求得．表6.6给出了部分Gauss-Legendre 求积公式的求积节点和求积系数．表6.6 Gauss-Legendre 求积公式中的数据和上一节描述的求积公式不同的是，Gauss-Legendre 求积公式的积分区间是]1,1[-，对于定积分⎰badx x f )(，可以通过变量代换t ab b a x 22-++=(6.63) 将区间],[b a 上的积分转化为]1,1[-上的积分⎰⎰-⎪⎭⎫⎝⎛-++-=1 1 222)(dt t a b b a f a b dx x f ba(6.64) 然后再使用Gauss-Legendre 求积公式进行计算．例6.8 分别使用两点Gauss-Legendre 求积公式和四点Gauss-Legendre 求积公式计算定 积分⎰2/ 0cos πxdx 的近似值．解 作变量代换)1(4t x +=π，则⎰⎰-⎪⎭⎫ ⎝⎛+==11 2/ 0)1(4cos 4cos dt t xdx I πππ117记⎪⎭⎫⎝⎛+=)1(4cos )(t t πφ，因为节点5773503.00=t ，5773503.01-=t ，110==A A ，由两点Gauss-Legendre 求积公式得160.99847260))()((410≈+≈t t I φφπ同样地，取160.861136310=t ，10.8611363116t =-，360.339981042=t ，360.339981043-=t ，50.3478548410==A A ，490.6521451532==A A ，由四点Gauss-Legendre 求积公式得720.99999997))()()()((433221100≈+++≈t A t A t A t A I φφφφπ易知，1cos 2/ 0==⎰πxdx I ．四点Gauss-Legendre 求积公式的误差为710720.999999971-≤-而由例6.3可知，要求误差不超过510-，复化梯形公式和复化Simpson 公式分别需要181个求积节点和11个求积节点．可见，Gauss 型求积公式具有较高的计算精度．但是，当求积节点数目增加时，Gauss 点发生了变化，前面计算的函数值不能在后面使用．三、Gauss 型求积公式的截断误差及稳定性定理6.8 设被积函数)(x f 的2n +2阶导函数在区间],[b a 上连续，则Gauss 求积公式的截断误差为),())(()!22()(][21)22(b a dx x n f f R b a n n ∈+=⎰++ξωξ (6.65)证明 设)(12x H n +是满足如下插值条件的不超过2n +1次的插值多项式),,1,0()()()()(1212n k x f x H x f x H k kn k k n =⎩⎨⎧'='=++(6.66)则其插值余项可表示为118 ),())(()!22()()()(21)22(12b a x n f x H x f n n n ∈+=-+++ηωη(6.67)由于Gauss 型求积公式具有12+n 次代数精度，因此∑∑⎰==++==ni i i n i i n i b an x f A x H A dx x H 01212)()()((6.68)故，求积公式截断误差可以表示为 ∑⎰--=ni i i b ax f A dx x f f R 0)()(][⎰⎰+-=ban b adxx H dx x f )()(12dxx n f b an n ⎰+++=21)22())(()!22()(ωη (6.69)因为21))((x n +ω在求积区间上不变号，)(x f 的2n +2阶导函数在求积区间上连续，对式(6.69)使用积分第二中值定理便得到式(6.65) ．对于Gauss 型求积公式，其求积系数220(())(())0nbk i k k ai A A l x l x dx ===>∑⎰(6.70)其中， n i i x l 0)}({=是以Gauss 点n i i x 0}{=为插值节点的Lagrange 插值基函数．上式表明Gauss求积系数都大于零，进而得到如下定理定理6.9 Gauss 型求积公式序列是稳定的、收敛的．。
