第4章数值分析

1、1第第4 4章章 数值积分与数值微分数值积分与数值微分4.1 数值积分概论4.2 牛顿-柯特斯公式 4.3 复合求积公式 4.4 龙贝格求积公式 4.5 自适应积分方法4.6 高斯求积公式4.7 多重积分4.8 数值微分 24.1 数值积分概论数值积分概论 数值积分的基本思想数值积分的基本思想 依据微积分基本定理，对于积分,)(badxxfI只要找到被积函数 的原函数 ，便有下列牛顿-莱布尼兹(Newton-Leibniz)公式： )( xf)( xF).()()(aFbFdxxfba但对于下列情形：3 （1）被积函数，诸如 等，找不到用初等函数表示的原函数，或者即使能求得原函数但原函数的表达

2、式非常复杂，计算困难； 2e),0(sinxxxx （2）当 是由测量或数值计算给出的一张数据表时，牛顿-莱布尼兹公式也不能直接运用. )( xf 因此有必要研究积分的数值计算问题. 由积分中值定理知，在积分区间 内存在一点，成立 ,ba),()()(fabdxxfba4就是说，底为 而高为 的矩形面积恰等于所求 ab )(f曲边梯形的面积 (图4-1).I图4-15 问题在于点的具体位置一般是不知道的，因而难以 准确算出 的值.)(f 将 称为区间 上的平均高度. )(f,ba 这样，只要对平均高度 提供一种算法，相应地便)(f获得一种数值求积方法. 用两端点“高度“ 与 的算术平均作为平均

3、高度 的近似值，这样导出的求积公式)(af)(bf)(f)()(2)(bfafabdxxfba（1.1）是梯形公式梯形公式(几何意义参看图4-2). 6图4-2 用区间中点 的“高度” 近似地取代平均高度 ，则又可导出所谓中矩形公式中矩形公式(简称矩形公式矩形公式) 2bac)(cf)(f).2()()(bafabdxxfba（1.2）7 一般地，可以在区间 上适当选取某些节点 ，,bakx然后用 加权平均得到平均高度 的近似值，)(kxf)(f, )()(0nkkkbaxfAdxxf（1.3）式中 称为求积节点求积节点； 称为求积系数求积系数，亦称伴随节点 的权权. kx 权 仅仅与节点 的

4、选取有关，而不依赖于被积函数 的具体形式. kAkx这样构造出的求积公式具有下列形式：kAkx)( xf8 这类数值积分方法通常称为机械求积，其特点是将积分求值问题归结为函数值的计算，这就避开了牛顿-莱布尼兹公式需要寻求原函数的困难. 9 代数精度的概念代数精度的概念 定义定义1 1 如果某个求积公式对于次数不超过 的多项式m均能准确地成立，但对于 次多项式就不准确成立，1m则称该求积公式具有 次代数精度次代数精度. m 梯形公式(1.1)和矩形公式(1.2)均具有一次代数精度. 数值求积是近似方法，为保证精度，自然希望求积公式对尽可能多的函数准确成立.10nkkabA0, 欲使求积公式(1.

5、3)具有 次代数精度，则只要令它m对 都准确成立，就得到mxxf,2, 1)(（1.4）nkmmmkkabmxA011).(11nkkkabxA022),(21, )()(0nkkkbaxfAdxxf（1.3）11 如果事先选定求积节点 ，譬如，以区间 的等距分点作为节点，这时取 ，求解方程组(1.4)即可确定求积系数 ，而使求积公式(1.3)至少具有 次代数精度. kx,banm kAn 构造形如(1.3)的求积公式，原则上是一个确定参数 和 的代数问题. kxkA, )()(0nkkkbaxfAdxxf（1.3）12 例如 时，取 ，求积公式为bxax10,1n).()()()(10bfA

6、afAdxxffIba在线性方程组（1.4）中令 ，则得1m),(21,221010abbAaAabAA解得 于是得).(2110bfAA).()(2)()(bfafabdxxffIba这就是梯形公式，表明利用线性方程组（1.4）推出的求积公式，与用通过两点 与 的直线近似曲线 得到的结果是一致的.)(,(afa)(,(bfb)( xfy 13 当 时（1.4）式的第3个式子不成立，因为2)(xxf).(31)(233222abdxxbaabba所以梯形公式（1.1）的代数精度为1. 在（1.4）中如果节点和系数都不确定，那么（1.4）就是关于 及 的 个参数的非线性方程组，该方程组在 时求解

7、是很困难的.ix), 1 ,0(niAi22n1n 但在 和 时还是可以通过求解（1.4）得到相应的求积公式的.0n1n140n 如 ，此时求积公式为),()()(00 xfAdxxffIba其中， 及 为待定参数.0A0 x 根据代数精度的定义可令 ，由（1.4）知xxf,1)(),(21,22000abxAabA于是 ).(210bax 所得到的就是（1.2）式的中矩形公式. 15 再令 ，代入（1.4）的第3式有2)(xxf),(31)(4)2)(332222200abdxxbaabbaabxAba说明公式（1.2）对 不精确成立，故它的代数精度为1.2)(xxf 方程组（1.4）是根据

8、形如（1.3）式的求积公式得到的，按照代数精度的定义，如果求积公式中除了 还有 在某些节点上的值，也同样可得到相应的求积公式. )(ixf)(xf 16 例例1 1 给定形如 的求积公式，试确定系数 ，使公式具有尽可能高的代数精度.)0() 1 ()0()(01010fBfAfAdxxf010,BAA 解解 根据题意可令 分别代入求积公式使它精确成立2, 1)(xxxf 当 时，得1)(xf;111010dxAA 当 时，得xxf)(;211001dxxBA17 当 时，得2)(xxf.311021dxxA解得 ,于是得61,32,31000BAA).0(61)1(31)0(32)(10fff

9、dxxf 当 时， 而上式右端为 ，故公式对 不精确成立，其代数精度为2.3)(xxf.41103dxx313)(xxf18 插值型的求积公式插值型的求积公式 设给定一组节点 ,210bxxxxan且已知函数 在这些节点上的值，)(xf作插值函数 .)(xLn取 banndxxLI)(作为积分 的近似值，badxxfI)(nkkknxfAI0)(（1.5）这样构造出的求积公式称为是插值型插值型的，式中求积系数 通过插值基函数 积分得出 kA)( xlk19., 1 ,0,)(nkdxxlAbakk（1.6）求积公式的值余项式中依赖于 ， x).()()(101nnxxxxxxx（1.7）,)(

10、)()(dxxRdxxLxffRbanban),()!1()()(1)1(xnfxRnnn其中20 当 是次数不超过 的多项式时，插值多项式就是n)( xf函数本身，余项 为零， fR. )()(0banjjkjkxlAdxxl 反之，如果求积公式(1.5)至少具有 次代数精度，则n它必定是插值型的. 事实上，这时公式(1.5)对于插值基函数 应准确)(xlk成立，即有至少具有 次代数精度.n所以这时插值型求积公式nkkknxfAI0)(（1.5）21注意到 上式右端实际上即等于 ，因而,)(kjjkxlkA.)(bakkdxxlA成立. 这样，有 定理定理1 1 形如(1.5)的求积公式至少

11、有 次代数精度的充分必要条件是，它是插值型的. nnkkknxfAI0)(（1.5）22 若求积公式(1.3)的代数精度为 ，则有求积公式余项的表达式(1.7)可以证明余项形如m（1.8）),()()()1(0mnkkkbafKxfAdxxffR其中 为不依赖于 的待定参数,K)(xf).,(ba 结果表明当 是次数小于等于 的多项式时，由于 ，故此时 ，即求积公式(1.3)精确成立.)(xfm0)()1(xfm0fR 而当 时， (1.8)的右端故可求得1)(mxxf,)!1()()1(mxfm, 0 xRn 求积公式的余项求积公式的余项23（1.9）.)()2(1)!1(1)!1(1012

12、2011nkmkkmmnkmkkbamxAabmmxAdxxmK代入余项(1.8)中可以得到更细致的余项表达式. 梯形公式(1.1)的代数精度为1，可以证明它的余项表达式为),(),(bafKfR 其中.)(121)(6121)(2)(3121332233ababbaababK于是得到梯形公式(1.1)的余项为).,(),(12)(3bafabfR （1.10）24 对中矩形公式(1.2)，其代数精度为1，可以证明),(),(bafKfR 其中.24)()2)()(31213233abbaababK于是得到梯形公式(1.1)的余项为).,(),(24)(3bafabfR （1.11）25 例例

13、2 2 求例1中求积公式)0(61) 1 (31)0(32)(10fffdxxf的余项 解解 由于此求积公式的代数精度为2，故余项表达式为 . 令 ，得 ，于是有)(fKfR 3)(xxf! 3)( f.721)3141(! 31)0(61) 1 (31)0(32(! 31103fffdxxK故得).1 , 0(),(721 ffR26 求积公式的收敛性与稳定性求积公式的收敛性与稳定性 定义定义2 2 在求积公式(1.3)中，若.)()(lim00bankkkhndxxfxfA其中),(max11iinixxh 在求积公式(1.3)中，由于计算 可能产生误差 ，)(kxfk实际得到将是 ，kf

14、即.)(kkkfxf, )()(0nkkknxfAfI则称求积公式(1.3)是收敛收敛的. 记nkkknfAfI0.)(, )()(0nkkkbaxfAdxxf（1.3）27如果对任给小正数,0只要误差 充分小就有 knkkkknnfxfAfIfI0)()()(（1.12）,则表明求积公式(1.3)计算是稳定的，由此给出: 定义定义3 3 对任给 若 只要),1 ,0()(nkfxfkk就有(1.12)成立，则称求积公式(1.3)是稳定稳定的. ,0,0, )()(0nkkkbaxfAdxxf（1.3）28 定理定理2 2 若求积公式(1.3)中系数 证明证明取,ab ,)(kkfxf),1

15、,0(0nkAk则此求积公式是稳定的. ,0对任给都有nk,1 ,0若对则当 时有0kAnkkkknnfxfAfIfI0)()()(nkkkkfxfA0)(, )()(0nkkkbaxfAdxxf（1.3）29 定理2表明，只要求积系数 ，就能保证计算的稳定性. 由定义3 ，知求积公式(1.3)是稳定的. nkkA0)(ab .0kA, )()(0nkkkbaxfAdxxf（1.3）304.2 牛顿牛顿- -柯特斯公式柯特斯公式 柯特斯系数与辛普森公式柯特斯系数与辛普森公式 设将积分区间 划分为 等分，,ban选取等距节点 构造出的插值型求积公式khaxknkknknxfCabI0)()()(

16、（2.1）称为牛顿牛顿- -柯特斯公式柯特斯公式，式中 称为柯特斯系数柯特斯系数. )( nkC 按(1.6)式，引进变换,thax,nabh步长则利用等距节点的插值公式，有.)(bakkdxxlA（1.6）31 nnkjjnkdtjkjtabhC00)(（2.2）.)()!(!)1(00 nnkjjkndtjtknnk 当 时，1n,21)1(1)1(0 CC这时的求积公式就是梯形公式)()(2bfafabT32 当 时，按(2.2)式，2n,61)2)(1(4120)2(0dtttC相应的求积公式是辛普森辛普森(Simpson)公式公式 ),()2(4)(6bfbafafabS（2.3）,

17、64)2(2120)2(1dtttC.61)1(4120)2(2dtttC柯特斯系数为 nnkjjnkdtjkjtabhC00)(（2.2）.)()!(!)1(00 nnkjjkndtjtknnk33 的牛顿-柯特斯公式称为柯特斯公式，其形式是 4n),(7)(32)(12)(32)(79043210 xfxfxfxfxfabC（2.4）这里 .4,abhkhaxk 按（2.2）式，可构造柯特斯系数表.3428350989228350588828350928283501049628350454028350104962835092828350588828350989817280751172803

18、57717280132317280298917280298917280132317280357717280751784041359280910534280935984041628819962514425144259625288195907451615245169074818383813613261221211)(nkCn1表435 从柯特斯系数表看到 时，柯特斯系数 出现负值，8n)( nkC,10)(0)(nknknknkCC特别地，假定,0)()(kknkfxfCnkkknknnfxfCfIfI0)()()()(于是有,)(kkfxf且则有 nkkknkfxfC0)()(nkkknkfxf

19、C0)()(36它表明初始数据误差将会引起计算结果误差增大，即计算不稳定，故 的牛顿-柯特斯公式是不用的. 8n.0)(nknkC37 偶阶求积公式的代数精度偶阶求积公式的代数精度 由定理1， 阶的牛顿-柯特斯公式至少具有 次的代数精度. nn 先看辛普森公式(2.3)，它是二阶牛顿-柯特斯公式，因此至少具有二次代数精度. 用 进行检验，3)(xxf.)2(46333bbaaabS本节讨论代数精度的进一步提高问题. 按辛普森公式计算得 ),()2(4)(6bfbafafabS（2.3）38这时有 ，即辛普森公式对次数不超过三次的多项式均能准确成立，而它对 通常是不准确的，因此，辛普森公式实际上

20、具有三次代数精度.4443abdxxIba另一方面，直接求积得 IS 4)(xxf 定理定理3 3 当阶 为偶数时，牛顿-柯特斯公式(2.1)至少有 次代数精度. n1nnkknknxfCabI0)()()(（2.1）39 证明证明 我们只要验证，当 为偶数时，牛顿-柯特斯公式对 的余项为零. n1)(nxxf 由于这里,)!1()()1(nxfn.)(0 banjjdxxxfR引进变换 并注意到 有 ,thax,jhaxj所以按余项公式有,)(002 nnjndtjthfR若 为偶数，则 为整数，n2n,2nut再令进一步有40,)2(2202 nnnjndujnuhfR因为被积函数.0fR

21、njjnuuH0)2()(为奇函数，所以2/2/)(nnjju41 辛普森公式的余项辛普森公式的余项 对牛顿-柯特斯求积公式，通常只使用 时的三个公式， 时为梯形公式，余项为 为辛普森公式4,2,1n1n).,(),(12)(3bafabfR 2n),()2(4)(6bfbafafabS代数精度为3，可以证明余项表达式为),(),()4(bafKfR其中 由(1.9)及(2.3)可得K42,)2(180120)(! 4)2(46)(51! 414544455abababbbaaababK).,(),()2(180)4(4bafababfR从而可得辛普森公式的余项为（2.5） 为柯特斯公式4n代

22、数精度为5，可以证明余项),(7)(32)(12)(32)(79043210 xfxfxfxfxfabC).,(),()4(945)(2)6(6bafababfR（2.6）434.3 复合求积公式复合求积公式 复合求积的基本思想是把积分区间分成若干子区间(通常是等分)，再在每个子区间上用低阶求积公式，目的是提高精度. 复合梯形公式复合梯形公式 将区间 划分为 等分，分点 在每个子区间 上采用梯形公式(1.1)，则得n, 1 , 0nk,nabhkhaxk)()(2bfafabT（1.1）,ba)1, 1 ,0(,1nkxxkk44badxxfI)( 101)(nkxxkkdxxf（3.1）).

23、()()(2101fRxfxfhnnkkk记 101)()(2nkkknxfxfhT称为复合梯形公式复合梯形公式. . （3.2）,)()(2)(211nkkbfxfafh45 由(1.10) ，余项由于 , 且 ,)(2baCxf 10)(1nkkfn所以 使 ),(ba于是复合梯形公式余项为 )(min10knkf).(max10knkf ).,(,)(121103 kkknkknxxfhTIfR. )(1)(10 nkkfnf（2.5）.,)(12)(3baabfRT 46).(12)(2fhabfRn （3.3）误差是 阶，且当 时有 ,)(2baCxf,)(limbanndxxfT即

24、复合梯形公式是收敛的. 事实上只要设 ，就可以得到收敛性，因为只要将 改写为 .)()(21110nkknkknxfnabxfnabT2hnT,)(baCxf47 此外， 的求积系数为正，由定理2知复合梯形公式是稳定的. nT当 时，上式右端括号内的两个和式均收敛到积分 所以复化梯形公式(3.2)收敛.n.)(badxxf48 复合辛普森求积公式复合辛普森求积公式 将区间 分为 等分，在每个子区间 上采用辛普森公式(2.3)，若记 ，则得nbadxxfI)(记 1012/1)()(4)(6nkkkknxfxfxfhS（3.5） 101)(nkxxkkdxxf（3.4）).()()(4)(610

25、12/1fRxfxfxfhnnkkkk),()2(4)(6bfbafafabS（2.3）hxxkk212/1,ba,1kkxx),()(2)(4)(611102/1bfxfxfafhnkknkk称为复合辛普森求积公式复合辛普森求积公式. 49 由(2.5)，其余项于是当 时，,)(4baCxf).,()()2(180)()4(4bafhabSIfRnn（3.6）误差阶为 ，显然是收敛的. 与复合梯形公式相似有 ).,(, )()2(1801110)4(4kkknkknnxxfhSIfR4h).,(),()2(180)4(4bafababfR 实际上，只要 就有,)(baCxf.)(limban

26、ndxxfS 此外，由于 中求积系数均为正数，故知复合辛普森公式计算稳定. nS50 例例3 3 对于函数 ，给出 的函数表(见表4-2)，试用复合梯形公式(3.2)及复合辛普森公式(3.5)计算积分 xxxfsin)(8n,sin10dxxxI并估计误差. 解解将积分区间 划分为8等分，8414709.018771925.08/79088516.04/39361556.08/59588510.02/19767267.08/39896158.04/19973978.08/110)(xfx2-表4应用复化梯形法求得 1 , 0;9456909.08T,)()(2)(211nkknbfxfafhT

27、),()(2)(4)(611102/1bfxfxfafhSnkknkkn51 同积分的准确值 比较，复合梯形法的结果 只有两位有效数字，而复合辛普森法的结果 却有6位有效数字.而如果将 分为4等分，应用复化辛普森法有 接下来看误差估计 ，由于,)cos(sin)(10dtxtxxxf所以有 以上得到的两个结果 与 ，都需要提供9个点上的函数值，计算量基本相同，然而精度却差别很大. 1 , 0.9460832.04S8T4S9460831.0I;9456909.08T.9460832.04S5210)()(cos)(dtxtdxdxfkkk，10)2cos(dtkxttk于是 10)2cos(d

28、ttkxtk由(3.3)得复合梯形公式误差 )(max)(10 xfkx10dttk.11k.000434.031)81(121)(max12)(210288 xfhTIfRx对复合辛普森公式，由(3.6)得 .10271.051)41(28801)(6444SIfR).(12)(2fhabfRn （3.3）).()2(180)()4(4fhabSIfRnn（3.6）53 例例4 4 计算积分 ，若用复合梯形公，问区间 应分多少等份才能使误差不超过 ，若改用复合辛普森公式，要达到同样的精度，区间 应分多少等份？,e10dxIx1 ,01 ,051021 解解 本题只要根据 及 的余项公式即可求

29、得其截断误差应满足的精度.nTnS 由于 ，由复合梯形公式的余项公式得误差的上界为1,e)(,e)(,e)()4( abxfxfxfxxx.1021e)1(121)(12)(522 nfhabfR54因此有 ，可取 ，即将区间 213等份，即可使误差不超过85.212,106e52nn213n1 ,0.10215 若采用复合辛普森公式计算积分，则由余项公式，要满足精度要求，必须使,1021e)1(28801)(2880)(54)4(4nfhabfRn由此得.707.3,10144e44nn可取 ，即用 的复合辛普森公式计算即可达到精度要求，此时区间 实际上应分为8等份.4n4n1 ,055 从

30、这个例子可以看出，为达到同样的精度，复合辛普森公式只需计算9个函数值，而复合梯形公式则需214个函数值，工作量相差近24倍.564.4 龙贝格求积公式龙贝格求积公式 梯形法的递推化梯形法的递推化 复化求积方法可提高求积精度，实际计算时若精度不够可将步长逐次分半. 设将区间 分为 等分，共有 个分点，,ban1n如果将求积区间再二分一次，则分点增至 个，12n我们将二分前后两个积分值联系起来加以考察. 57用复合梯形公式求得该子区间上的积分值为 每个子区间 经过二分只增加了一个分点这里 代表二分前的步长. nabh 将每个子区间上的积分值相加得 , )(2)()(410211012nkknkkk

31、nxfhxfxfhT,1kkxx)(21121kkkxxx).()(2)(4121kkkxfxfxfh58从而利用式(3.2)可导出下列递推公式 . )(22110212nkknnxfhTT（4.1）（3.2）,)()(2)(211nkknbfxfafhT59 例例5 5.sin10dxxxI 解解先对整个区间 使用梯形公式.对于函数计算积分值 ,sin)(xxxf定义它在 的值0 x, 1)0(f,8414709.0)1 (f而由梯形公式 1 , 0.9207355.0)1()0(211ffT将区间二等分，求出中点的函数值 ,9588510.0)21(f60利用递推公式(4.1)，有 .93

32、97933.0)21(212112fTT进一步二分求积区间，并计算新分点上的函数值 .9088516.0)43(,9896158.0)41(ff再利用式(4.1)，有 . )(22110212nkknnxfhTT.9445135.0)43()41(412124ffTT这样不断二分下去，计算结果见下表. 619460831.09460830.09460827.09460815.09460769.01098769460596.09459850.09456909.09445135.09397933.054321nnTkTk3-表4 它表明用复化梯形公式计算积分 要达到7位有效数字的精度需要二分区间1

33、0次，即要有分点1025个，计算量很大. I62 外推技巧外推技巧 由梯形公式出发，将区间 逐次二分可提高求积公式的精度，当 分为 等份时，有,ba),(122fhabTIn 若记 当区间 分为 等份时，则有并且有),(hTTn,ban2),(12)(2fhabIhT 梯形公式余项可展成级数形式，即 ,ba.nabh),2(2hTTn,)0()(lim0IThTh,ban63 定理定理4 4 设 则有 ,)(baCxf,)(24221llhhhIhT（4.2）其中系数 与 无关. ),2,1(llh 定理4表明 是 阶，在(4.2)中，若用代替 ，有IhT)()(2hOh,2164)2(242

34、21llhhhIhT（4.3）若用4乘(4.3)式，减去(4.2)式再除3并记之为 则有 ),(hS2h64 这种将计算 的近似值的误差阶由 提高到 的方法称为外推算法，也称为(Richardson)外推算法.3)(2/4)(hThThS这里 是与 无关的系数.,21hI)(4hO（4.4）,6241hhI 用 近似积分值 ，其误差阶为 ，这比复合梯形公式的误差阶 提高了，容易看到 ，即将 分为 等份得到的复合辛普森公式. )(hSI)(4hO)(2hOnShS)(,ban)(2hO 只要真值与近似值的误差能表示成 的幂级数，如(4.2)，都可以使用外推算法提高精度.h65与上述做法类似，从(

35、4.4)出发，当 再增加一倍，即 减少一半时，有hn.)2()2()2(6221hhIhS（4.5）用16乘(4.5)式再减去(4.4)式后除以15，将所得的式子记为 ，则有 )(hC15)(2/16)(hShShC（4.6）.8261hhI它就是把区间 分为 个子区间的复合柯特斯公式， ，它的精度为,ban.)()(6hOIhCnChC)(66 这个公式相当于由辛普森法二分前后的两个积分值 与 组合得到的，即nS)2(2hSSn).16(1512nnnSSC（4.7）).()2(64631)(hChChR（4.8） 从(4.6)出发，利用外推技巧还可得到逼近阶为 的算法公式)(8hO如此继续

36、下去就可得到龙贝格(Romberg)算法.67).(1412144)(11hThThTmmmmmm（4.9）经过 次加速后，余项便取下列形式： ),2, 1(mm.)()2(22)1(21mmmhhIhT（4.10） 上述处理方法通常称为理查森外推加速方法理查森外推加速方法. 龙贝格算法龙贝格算法 将上述外推技巧得到的公式(4.4)、(4.6)、(4.8)重新引入记号)()(),()(),()(),()(3210hRhThChThShThThT从而可将上述公式写成统一形式68 设以 表示二分 次后求得的梯形值，且以 表示)(0kTk)(kmT序列 的 次加速值，则依递推公式(4.9)可得 )(

37、0kTm公式(4.11)也称为龙贝格求积算法龙贝格求积算法. . （4.11）).,2, 1(141144)(1)1(1)(kTTTkmmkmmmkm69计算过程： (1) 取,0abhk 令k1( 记区间 的二分次数). k,ba (2) 求梯形值,20kabT即按递推公式(4.1)计算.)(0kT (3) 求加速值，按公式(4.12)逐个求出T表(见表4-4)的第 行其余各元素 k).,2, 1()(kjTjkj).()(2)0(0bfafhT求 (4) 若 (预先给定的精度)，则终止计算，)0(1)0(kkTT;)0(ITk并取否则令 转(2)继续计算. kk1. )(22110212n

38、kknnxfhTT（4.1）70)0(410) 1 (39)2(28)3(17)4(0)0(36) 1 (25)2(14)3(0)0(23) 1 (12)2(0)0(11) 1 (0)0(0)(4)(3)(2)(1)(01648342210TTTTTabTTTTabTTTabTTabTabTTTTThkkkkkk表T4-4表71 可以证明，如果 充分光滑，那么T表每一列的元素及对角线元素均收敛到所求的积分值 ，即 )(xfI，）（固定mITkmk)(lim 对于 不充分光滑的函数也可用龙贝格算法计算，只是收敛慢一些，这时也可以直接使用复合辛普森公式计算. )( xf.lim)0(ITmm72

39、例例6 6 用龙贝格算法计算积分 解解在 上仅是一次连续可微，2/3)(xxf1 ,0.102/3dxxI用龙贝格算法计算结果见表4-5. 400002. 0400002. 0400002. 0400002. 0400002. 0400118. 05400009. 0400009. 0400009. 0400014. 0400463. 04400050. 0400054. 0400077. 0401812. 03400302. 0400432. 0407018. 02402369. 0426777. 01500000. 00)(5)(4)(3)(2)(1)(0kkkkkkTTTTTTk5-表4

40、73 从表中看到用龙贝格算到 的精度与辛普森求积精度相当. 这里 的精确值为 5kI.4.0744.5 自适应积分方法自适应积分方法 复合求积方法是用于被积函数变化不太大的积分. 如果在求积区间中被积函数变化很大，有的部分函数值变化剧烈，另一部分变化平缓，这时统一将区间等份用复合求积公式计算工作量就会很大. 要达到误差要求对变化剧烈部分必须将区间细分，而平缓部分则可用大步长，即针对被积函数在区间上不同情形采用不同的步长，使得在满足精度前提下积分计算的工作量尽可能小.75 针对这类问题的算法技巧是在不同区间上预测被积函数变化的剧烈程度确定相应的步长. 这种方法称为自适应积分方法. 以常用的复合辛

41、普森公式为例说明方法的基本思想.76 设给定精度要求 ，计算积分的近似值.先取步长 ，应用辛普森公式有其中若把区间 对分，步长 ，在每个小区间上用辛普森公式，则得0badxxfI)(abh),(),()2(180),()()4(4bafhabbaSdxxfIba（5.1）).()2(4)(6),(bfbafafhbaS,ba222abhh),(),()2(180),()()4(422bafhabbaSfI（5.2）77实际上(5.2)即为与(5.1)比较，若 在 上变化不大，可假定其中).()43(4)2(6),2(),2()4(4)(6)2,(),2()2,(),(222bfhafhafhb

42、baShafhafafhbaaSbbaSbaaSbaS).,(),()4(180),()()4(42bafhabbaSfI（5.2）)()4(xf),(ba)()()4()4(ff从而可得).()2(180),(),(1516)4(42fhabbaSbaS78 若不等式(5.3)不成立，则应分别对子区间 及 再用辛普森公式，此时步长 ,得到这里 .如果有则可期望得到与(5.2)比较，则得,151),(),(151),()(2122SSbaSbaSbaSfI),(),(221baSSbaSS,1521 SS,),()(2baSfI),(2baS此时可取 作为 的近似，则可达到给定的误差精度 .b

43、adxxfI)(2,baa,2bba 2321hh79 及 . 只要分别考察 及是否成立. 对满足要求的区间不再细分，对不满足要求的还要继续上述过程，直到满足要求为止，最后还要应用龙贝格法则求出相应区间积分的近似值.)2,(3baaS),2(3bbaS2)2,(3baaSI2),2(3bbaSI80 例例7 7 计算积分 若用复合辛普森法(3.5)，计算结果见表4-6.（此处 即为公式中的 ，积分精确值为4） 计算到 为止，此时 的近似值 ，若再用龙贝格法则得到整个计算是将 做32等分，即需要计算33个 的值.112.02dxxInhh02.01nnSS.1)(12.02dxxfI000154

44、.41 ,2 .05S.00002.4151 ,2.0455SSSRS1, 2 .0)(xf81 现在若用自适应积分法，当 时有4 .02h,66851852.01,6 .0,51851852.36 .0,2 .022SS由于 大于允许误差，故要对 及 两区间再用 做积分.,187037.41,6 .06 .0,2 .01,2 .02222SSSS761111.021 SS6 .0, 2 .01, 6 .02321hh 先计算 的积分1, 6 .0.25002572.01, 8 .0,41678477.08 .0,6 .033SS 由于001708.066681049.066851852.0)

45、1, 8 .08 .0,6 .0(1,6 .0332SSS82小于允许误差0.01，故在 区间的积分值为下面再计算子区间 的积分，其中 由于1, 6 .0.66669662.0)66851852.066681049.001,6 .0RS6 .0, 2 .0而对 可求得,51851852.36 .0,2 .02S2321hh.83425926.06 .0,4 .0,52314815.24 .0,2 .033SS161111.0)6 .0,4 .04 .0,2 .0(6 .0,2 .0332SSS83大于允许误差0.01，因此还要分别计算 及的积分. 当 时可求得而小于允

46、许误差0.005，故可得 的积分近似而对区间 ，其误差 不小于0.005，故还要分别计算 及 的积分，4 .0, 2 .06 .0, 4 .0234hh,33334864.06 .0, 5 .0,50005144.05 .0,4 .044SS000859.0)6 .0, 5 .05 .0,4 .0(6 .0,4 .0443SSS6 .0, 4 .0.8333428.06 .0,4 .0RS4 .0, 2 .04 .0,2 .04 .0,2 .043SS3 .0, 2 .04 .0, 3 .084且小于允许误差0.0025，故有最后子区间 的积分可检验出它的误差小于0.0025，且可得其中 ，当

47、 可求得83356954.04 .0, 3 .04S245hh,35714758.04 .0,35.0,47620166.035.0, 3 .055SS000220.0)4 .0,35.035.0, 3 .0(4 .0, 3 .0554SSS,83333492.04 .0, 3 .0RS3 .0, 2 .0.666686.13 .0,2 .0RS85将以上各区间的积分近似值相加可得它一共只需计算17个 的值.,00005957.41,6 .06 .0,4 .04 .0, 3 .03 .0,2 .0)(RSRSRSRSfI)(xf86 一般理论一般理论 求积公式 nkkkbaxfAdxxf0)(

48、)(含有 个待定参数22n).,1 ,0(,nkAxkk 当 为等距节点时得到的插值求积公式其代数精度至少为 次. kxn 如果适当选取 有可能使求积公式具有 次代数精度.), 1 ,0(nkxk12n4.6 高斯求积公式高斯求积公式87试确定节点 及 和系数 ，使其具有近可能高的代数精度.（6.2）.0131030AxAx 例例8 8 求积公式),()()(110011xfAxfAdxxf（6.1） 解解 令公式(6.1)对于 准确成立，得32, 1)(xxxxf,210 AA,00000AxAx,32121020AxAx0 x1x10, AA88 用(6.2)式中的第3式减去 乘(6.2)

49、中的第2式有0 x用第4式减去第2式乘 ，得 ,0)(202111 xxxA由此得于是可取 .20 x.01xx.32)(0111 xxxA用前一式代入则得,3110 xx由此得出 与 异号，即 ，从而有0 x1x01xx.31, 1211xA33,3310 xx89 再由(6.2)式的第1式得 ，于是有110 AA).33()33()(11ffdxxf（6.3）当 时，(6.3)式两端分别为 及 ，(6.3)式对 不精确成立，故公式(6.3)的代数精度为3.4)(xxf52924)(xxf 实际上，形如(6.1)的求积公式其代数精度不可能超过3，因为当 时，设这是4次多项式，代入(6.12)

50、左端有 ，而右端为0.表明两个节点的求积公式的代数精度为3.1 , 1,10 xx2120)()()(xxxxxf0)(11dxxf 一般 节点的求积公式的代数精度最高为 次.1n12n90kx为求积节点，可适当选取 及 使5.1）具有 次代数精度. 下面研究带权积分 这里 为权函数，类似(1.3)，求积公式为,)()(badxxxfI)( x, )()()(0nkkkbaxfAdxxxf（6.4）为不依赖于 的求积系数.), 1 , 0(nkAk)( xf12n), 1 ,0(nkxk, ),1 ,0(nkkA, )()(0nkkkbaxfAdxxf（1.3）91 根据定义要使(6.4)具有

51、 次代数精度，只要对12n),12,1 ,0(,)(nmxxfm.12, 1 ,0)(0nmdxxxxAnkbammkk（6.5）当给定权函数 ，求出右端积分，则可由(6.5)解得 )( x).,1 ,0(nkAxkk及令(6.4)精确成立，即, )()()(0nkkkbaxfAdxxxf（5.1） 定义定义4 4 如果求积公式(6.4)具有 次代数精度，则称其节点 为高斯点高斯点，相应公式(6.4)称为高斯求积公式高斯求积公式.12n),1 ,0(nkxk92 (6.4)是关于 及 的非线性方程组，当 时求解非常困难.), 1 ,0(nkxkkA1n 如果事先确定了节点 ，则可以利用(6.5

52、)求解 .此时(6.5)是关于 的线性方程组.), 1 ,0(nkxk), 1 ,0(nkAkkA 下面讨论如何选取节点 才能使求积公式(6.4)具有 次代数精度.), 1 ,0(nkxk12n93 设 上的 个节点的拉格朗日插值多项式为)(,10 xfbxxxan1n, )()()(0nkkknxlxfxL,ba其中,)()()(0nkjjjkjkxxxxxl则).,(),()()!1(1)()()(1)1(0baxxfnxlxfxfnnnkkk94用 乘上式并从 到 积分，则得,)()()()!1(1)()()(1)1(0dxxxxfnxfAdxxfxbannnkkkba)(xab其中,)

53、()(dxxxlAbakk余项.)()()()!1(11)1(dxxxxfnfRbann95显然当 取为 时有 ，此时有)(xfnxx, 10fR, )()()(0nkkkbaxfAdxxfx即求积公式至少具有 次代数精度.n 现在考察如何选取节点 才能使求积公式精度提高到 次.), 1 ,0(nkxk12n 此时要求 为 次多项式时 ，而当 时， 为 次多项式.)(xf12n0fR12)(nHxf)()1(xfnn96 若要求对 ，积分nHxp)(,0)()()(1dxxxxpban即相当与要求 与每个 带权 在 上 正交.)(1xnnHxp)()(x,ba 也就是以节点 为零点的 次多项式

54、 是 上带权 正交的多项式，故有以下定理.), 1 ,0(nkxk1n)(1xn,ba)(x97 定理定理5 5 插值型求积公式(6.4)的节点bxxxan10是高斯点的充分必要条件是以这些节点为零点的多项式)()()(101nnxxxxxxx与任何次数不超过 的多项式 带权 正交，n)( xP)( x.0)()()(1bandxxxxp（6.7） 证明证明 必要性. 即,H)(nxp设,H)()(121nnxxp则, )()()(0nkkkbaxfAdxxxf（6.4）98.)()()()(babadxxxqdxxxf（6.8）nxxx,10是高斯点，因此，如果. )()()()()(011

55、nkknkkbanxxpAdxxxxp因),1 ,0(0)(1nkxkn即有故(6.7)成立. 精确成立，)()()(1xxpxfn 则求积公式(6.4)对于 充分性. 用 除 ，)(1xn )( xf记商为 ，余式为 ，即 ，)( xp)(xq)()()()(1xqxxpxfn其中 . nxqxpH)(),(,H)(12nxf对于由(6.7)可得 , )()()(0nkkkbaxfAdxxxf（6.4）.0)()()(1bandxxxxp（6.7）99由于求积公式(6.4)是插值型的，它对于 是精确的，nxqH)(. )()()(0nkkkbaxqAdxxxq即 再注意到),1 ,0(0)(

56、1nkxkn), 1 ,0()()(nkxfxqkk知从而由(6.8)有babadxxxqdxxxf)()()()(. )(0nkkkxfA.)()()()(babadxxxqdxxxf（6.8）, )()()(0nkkkbaxfAdxxxf（6.4）100可见求积公式(6.4)对一切次数不超过 的多项式均精确成立. 因此， 为高斯点. 12n),1 ,0(nkxk 定理表明在 上带权 的 次正交多项式的零点就是求积公式(6.4)的高斯点. ,ba)( x1n 有了求积节点 ，再利用),1 ,0(nkxknkbammkkdxxxxA0)(对 成立，nm,1 ,0).,1 ,0(nkAk解此方程

57、则得 的线性方程.nAAA,10则得到一组关于求积系数101 也可直接由 的插值多项式求出求积系数 nxxx,10).,1 ,0(nkAk102 例例9 9 确定求积公式 ).()()(110010 xfAxfAdxxfx 解解 具有最高代数精度的求积公式是高斯型求积公式，其节点为关于权函数 的正交多项式零点 及 ,的系数 及节点 和 ，使它具有最高的代数精度.0 x1x10, AAxx )(0 x1x 设 由正交性知与1及 带权正交，即得,)()(210cbxxxxxxx)(xx.0)(,0)(1010dxxxxdxxx于是得,0527292,0325272cbcb103由此解得 即 ,21

58、5,910cb.215910)(2xxx令 ，则得 0)(x.821162.0,289949.010 xx 由于两个节点的高斯求积公式具有3次代数精度，故公式对 精确成立，即xxf, 1)( 当 时1)(xf;321010dxxAA.52100000dxxxAxAx 当 时xxf)(由此解出.389111.0,277556.010AA104 下面讨论高斯求积公式(6.4)的余项. 利用 在节点 的埃尔米特插值)( xf),1 ,0(nkxk,12nH., 1 , 0),()(),()(1212nkxfxHxfxHkknkkn于是 )()!22()()(21)22(12xnfHxfnnn即 ,

59、)()()(0nkkkbaxfAdxxxf（6.4）105两端乘 ，并由 到 积分，则得 )( xab.)()()()(12fRdxxxHdxxxfInbanba（6.9）其中右端第一项积分对 次多项式精确成立，故 12nnkkknxfAIfR0)(由于,0)()(21xxn.)()()!22()(21)22(bannndxxxnffR（6.10）.)()()!22()(21)22(banndxxxnf由积分中值定理得(6.4)的余项为 关于高斯求积公式的稳定性与收敛性，有： , )()()(0nkkkbaxfAdxxxf（6.4）106 定理定理6 6 高斯求积公式(6.4)的求积系数全是正

60、的. 证明证明 考察,)(0nkjjjkjkxxxxxl它是 次多项式，n因而 是 次多项式，)(2xlkn2. )()()(0022niikibakxlAdxxxl注意到,)(kiikxl),1 ,0(nkAk故高斯求积公式(6.4)对于它能准确成立，即有 ,kA上式右端实际上即等于从而有 , )()()(0nkkkbaxfAdxxxf（6.4）107由本定理及定理2，则得 推论推论 高斯求积公式(6.4)是稳定的. 定理定理7 7 设 则高斯求积公式(6.4)收敛，即.)()()(lim0bankkkndxxxfxfA.0)()(2bakkdxxxlA定理得证. ,)(baCxf, )()