1、14.1 4.1 引言引言 4.1.1 4.1.1 数值求积的基本思想数值求积的基本思想 依据微积分基本定理，对于积分,)(badxxfI只要找到被积函数 的原函数 ，便有下列牛顿-莱布尼兹(Newton-Leibniz)公式：)(xf)(xF).()()(aFbFdxxfba但对于下列情形：2 （1）被积函数，诸如 等等，找不到用初等函数表示的原函数；2sin,sinxxx （2）当 是由测量或数值计算给出的一张数据表时，牛顿-莱布尼兹公式也不能直接运用.)(xf 因此有必要研究积分的数值计算问题.由积分中值定理知，在积分区间 内存在一点，成立,ba),()()(fabdxxfba3就是说，
2、底为 而高为 的矩形面积恰等于所求 ab)(f曲边梯形的面积 (图4-1).I图4-14 问题在于点的具体位置一般是不知道的，因而难以 准确算出 的值.)(f 将 称为区间 上的平均高度.)(f,ba 这样，只要对平均高度 提供一种算法，相应地便)(f获得一种数值求积方法.用两端点“高度“与 的算术平均作为平均高度)(af)(bf)(f的近似值，这样导出的求积公式)()(2bfafabT（1.1）是梯形公式(几何意义参看图4-2).5图4-2 用区间中点 的“高度”近似地取代平均高度 ，则又可导出所谓中矩形公式(简称矩形公式)2bac)(cf)(f).2()(bafabR（1.2）6 一般地，
3、可以在区间 上适当选取某些节点 ，,bakx然后用 加权平均得到平均高度 的近似值，)(kxf)(f,)()(0nkkkbaxfAdxxf（1.3）式中 称为求积节点求积节点；kx权 仅仅与节点 的选取有关，kAkx的具体形式.)(xf这样构造出的求积公式具有下列形式：kA称为求积系数求积系数，亦称伴随节点kx的权权.而不依赖于被积函数7 这类数值积分方法通常称为机械求积，其特点是将积分求值问题归结为函数值的计算，这就避开了牛顿-莱布尼兹公式需要寻求原函数的困难.8 4.1.2 4.1.2 代数精度的概念代数精度的概念 定义定义1 1 如果某个求积公式对于次数不超过 的多项式m均能准确地成立，
4、但对于 次多项式就不准确成立，1m则称该求积公式具有 次代数精度次代数精度.m 梯形公式(1.1)和矩形公式(1.2)均具有一次代数精度.数值求积是近似方法，为保证精度，自然希望求积公式对尽可能多的函数准确成立.9nkkabA0,欲使求积公式(1.3)具有 次代数精度，则只要令它m对 都准确成立，就得到mxxf,2,1)(（1.4）nkmmmkkabmxA011).(11nkkkabxA022),(2110 如果事先选定求积节点 ，譬如，以区间 的等距分点作为节点，这时取 ，求解方程组(1.4)即可确定求积系数 ，而使求积公式(1.3)至少具有 次代数精度.kx,banm kAn 构造形如(1
5、.3)的求积公式，原则上是一个确定参数 和 的代数问题.kxkA11 4.1.3 4.1.3 插值型的求积公式插值型的求积公式 设给定一组节点,210bxxxxan且已知函数 在这些节点上的值，)(xf作插值函数 .)(xLn取 banndxxLI)(作为积分 的近似值，badxxfI)(nkkknxfAI0)(（1.5）这样构造出的求积公式12称为是插值型插值型的，式中求积系数 通过插值基函数 积分得出 kA)(xlk.)(bakkdxxlA（1.6）由插值余项定理(第2章的定理2)即知，对于插值型的求积公式(1.5)，其余项 nIIfR式中与变量 有关，x).()()(10nxxxxxxx
6、（1.7）,)()!1()()1(bandxxnf13 当 是次数不超过 的多项式时，插值多项式就是n)(xf函数本身，余项 为零，fR至少具有 次代数精度.n 反之，如果求积公式(1.5)至少具有 次代数精度，则n.)()(0banjjkjkxlAdxxl 事实上，这时公式(1.5)对于插值基函数 应准确)(xlk它必定是插值型的.成立，即有所以这时插值型求积公式14 定理定理1 1形如(1.5)的求积公式至少有 次代数精度的n注意到,)(kjjkxl上式右端实际上即等于 ，kA因而式(1.6).)(bakkdxxlA成立.这样，有充分必要条件是，它是插值型的.15 4.1.4 4.1.4
7、求积公式的收敛性与稳定性求积公式的收敛性与稳定性 定义定义2 2.)()(lim00bankkkhndxxfxfA其中),(max11iinixxh 在求积公式(1.3)中，由于计算 可能产生误差 ，)(kxfk实际得到将是 ，kf即.)(kkkfxf,)()(0nkkknxfAfI在求积公式(1.3)中，若则称求积公式(1.3)是收敛的.记nkkknfAfI0.)(16如果对任给小正数,0只要误差 充分小就有 knkkkknnfxfAfIfI0)()()(（1.8）,则表明求积公式(1.3)计算是稳定的，由此给出:定义定义3 3),1,0()(nkfxfkk就有(1.8)成立，则称求积公式(
8、1.3)是稳定的.,0对任给,0若只要17 定理定理2 2 证明证明取,ab,)(kkfxf),1,0(0nkAk若求积公式(1.3)中系数 则此求积公式是稳定的.,0对任给都有nk,1,0若对则当 时有0kAnkkkknnfxfAfIfI0)()()(nkkkkfxfA0)(18由定义3，知求积公式(1.3)是稳定的.nkkA0)(ab.194.2 4.2 牛顿牛顿-柯特斯公式柯特斯公式 4.2.1 4.2.1 柯特斯系数柯特斯系数 设将积分区间 划分为 等分，,ban选取等距节点 构造出的插值型求积公式khaxknkknknxfCabI0)()()(（2.1）称为牛顿牛顿-柯特斯公式柯特斯
9、公式，式中 称为柯特斯系数柯特斯系数.)(nkC 按(1.6)式，引进变换,thax,nabh步长则利用等距节点的插值公式，有20 nnkjjnkdtjkjtabhC00)(（2.2）.)()!(!)1(00 nnkjjkndtjtknnk 当 时，1n,21)1(1)1(0 CC这时的求积公式就是梯形公式(1.1)()(2bfafabT21 当 时，按(2.2)式，2n,61)2)(1(4120)2(0dtttC相应的求积公式是辛普森辛普森(Simpson)公式公式),()2(4)(6bfbafafabS（2.3）,64)2(2120)2(1dtttC.61)1(4120)2(2dtttC柯
10、特斯系数为 22 的牛顿-柯特斯公式称为柯特斯公式，4n),(7)(32)(12)(32)(79043210 xfxfxfxfxfabC（2.4）这里.4,abhkhaxk 按（2.2）式，可构造柯特斯系数表.其形式是 2328350989228350588828350928283501049628350454028350104962835092828350588828350989817280751172803577172801323172802989172802989172801323172803577172807517840413592809105342809359840416288199
11、62514425144259625288195907451615245169074818383813613261221211)(nkCn1表424 从柯特斯系数表看到 时，柯特斯系数 出现负值，8n)(nkC,10)(0)(nknknknkCC特别地，假定,0)()(kknkfxfCnkkknknnfxfCfIfI0)()()()(于是有,)(kkfxf且则有 nkkknkfxfC0)()(nkkknkfxfC0)()(25它表明初始数据误差将会引起计算结果误差增大，即计算不稳定，故 的牛顿-柯特斯公式是不用的.8n.0)(nknkC26 4.2.2 4.2.2 偶阶求积公式的代数精度偶阶求积
12、公式的代数精度 由定理1，阶的牛顿-柯特斯公式至少具有 次的代数精度.nn 先看辛普森公式(2.3)，它是二阶牛顿-柯特斯公式，因此至少具有二次代数精度.用 进行检验，3)(xxf.)2(46333bbaaabS本节讨论代数精度的进一步提高问题.按辛普森公式计算得 27均能准确成立，.4443abdxxIba另一方面，直接求积得 这时有 ，IS 而它对 通常是不准确的，4)(xxf辛普森公式实际上具有三次代数精度.即辛普森公式对次数不超过三次的多项式因此，定理定理3 3当阶 为偶数时，牛顿-柯特斯公式(2.1)至少n有 次代数精度.1n28 证明证明我们只要验证，当 为偶数时，牛顿-柯特斯n公
13、式对 的余项为零.1)(nxxf 由于这里,)!1()()1(nxfn.)(0 banjjdxxxfR引进变换 并注意到 有,thax,jhaxj按余项公式(1.7)nIIfR,)()!1()()1(bandxxnf有 29,)2(2202 nnnjndujnuhfR因为被积函数.0fR若 为偶数，则 为整数，n2n,)(002 nnjndtjthfRnjjnuuH0)2()(为奇函数，所以,2nut再令进一步有2/2/)(nnjju30 4.2.3 4.2.3 几种低阶求积公式的余项几种低阶求积公式的余项 按余项公式(1.7)，梯形公式(1.1)的余项,)(2)(baTdxbxaxfTIR这
14、里积分的核函数 在区间 上保号(非正)，)(bxax,ba应用积分中值定理，在 内存在一点 使,ba,baTdxbxaxfR)(2)(（2.5）.,)(12)(3baabf 31 为研究辛普森公式(2.3)的余项 构造次数,SIRs不超过3的多项式 满足 ),(xH),()(),()(bfbHafaH（2.6）).()(),()(cfcHcfcH其中.2bac 辛普森公式具有三次代数精度，对于这样构造出的三次式 应是准确的，即 )(xH),()(4)(6)(bHcHaHabdxxHba32basdxxHxfSIR)()(对于多项式 ，其插值余项由第2章(5.11)得)(xH),()(!4)()
15、()(2)4(bxcxaxfxHxf由插值条件(2.6)，上式右端实际上等于按辛普森公式(2.3)求得的积分值 ，S因此积分余项 故有.)()(!4)(2)4(basdxbxcxaxfR33这时积分的核函数 在 上保号)()(2bxcxax,babasdxbxcxaxfR)()(!4)(2)4(类似的，对于柯特斯公式(2.4)，结果如下：).(4945)(2)6(6fababCIRC（2.8）(非正)，再用积分中值定理有（2.7）).(2180)4(4fabab344.3 4.3 复化求积公式复化求积公式 复化求积的基本思想是把积分区间分成若干子区间(通常是等分)，再在每个子区间上用低阶求积公
16、式，目的是提高精度.4.3.1 4.3.1 复化梯形公式复化梯形公式 将区间 划分为 等分，,ban,1,0nk在每个子区间 上)1,1,0(,1nkxxkk,nabhkhxk分点采用梯形公式(1.1)，则得 35badxxfI)(101)(nkxxkkdxxf（3.1）).()()(2101fRxfxfhnnkkk记 101)()(2nkkknxfxfhT称为复化梯形公式复化梯形公式.（3.2）,)()(2)(211nkkbfxfafh36 由(2.5)，其余项nnTIfR)(由于 ,且,)(2baCxf 10)(1nkkfn所以 使),(ba.)(1)(10 nkkfnf于是复化梯形公式余
17、项为 )(min10knkf).(max10knkf 103,)(12nkkfh).,(1kkkxx37).(12)(2fhabfRn（3.3）误差是 阶，2h且当 时有,)(2baCxf,)(limbanndxxfT即复化梯形公式是收敛的.将 改写为 nT.)()(21110nkknkknxfnabxfnabT38 此外，的求积系数为正，由定理2知复化梯形公式是稳定的.nT只要 则当 时，上式右端括号内的两个和式均收敛到积分 所以复化梯形公式(3.2)收敛.,)(baCxfn,)(badxxf39 4.3.2 4.3.2 复化辛普森求积公式复化辛普森求积公式 将区间 分为 等分，,ban在每
18、个子区间 上,1kkxx若记,212/1hxxkkbadxxfI)(记 1012/1)()(4)(6nkkkknxfxfxfhS（3.5）采用辛普森公式(2.3)，则得 101)(nkxxkkdxxf（3.4）).()()(4)(61012/1fRxfxfxfhnnkkkk40),()(2)(4)(611102/1bfxfxfafhnkknkk称为复化辛普森求积公式复化辛普森求积公式.由(2.7)，其余项nnSIfR)(于是当 时，,)(4baCxf).()2(180)()4(4fhabSIfRnn（3.6）误差阶为 ，显然是收敛的.4h10)4(4),()2(180nkkfhh).,(1kk
19、kxx与复化梯形公式相似有 41 实际上，只要,)(baCxf.)(limbanndxxfS则可得到收敛性，即 此外，由于 中求积系数均为正数，故知复化辛普森公式计算稳定.nS42 例例1 1对于函数 ，xxxfsin)(给出 的函数表8n,sin10dxxxI并估计误差.解解将积分区间 划分为8等分，1,08414709.018771925.08/79088516.04/39361556.08/59588510.02/19767267.08/39896158.04/19973978.08/110)(xfx2-表4;9456909.08T(见表4-2)，计算积分 应用复化梯形法求得 试用复化梯
20、形公式(3.2)及复化辛普森公式(3.5)43而如果将 分为4等分，应用复化辛普森法有 1,0.9460832.04S 以上得到的两个结果 与 ，都需要提供9个点上的4S8T 同积分的准确值 比较，9460831.0I 接下来看误差估计，由于,)cos(sin)(10dtxtxxxf所以有 函数值，计算量基本相同，然而精度却差别很大.只有两位有效数字.复化梯形法的结果4410)()(cos)(dtxtdxdxfkkk，10)2cos(dtkxttk于是 10)2cos(dttkxtk由(3.3)得复化梯形公式误差 88)(TIfR)(max)(10 xfkx10dttk.11k)(max121
21、02xfhx 31811212.000434.045对复化辛普森公式，由(3.6)得 44)(SIfR5141288014.10271.06464.4 4.4 龙贝格求积公式龙贝格求积公式 4.4.1 4.4.1 梯形法的递推化梯形法的递推化 复化求积方法可提高求积精度，实际计算时若精度不够可将步长逐次分半.设将区间 分为 等分，共有 个分点，,ban1n如果将求积区间再二分一次，则分点增至 个，12n我们将二分前后两个积分值联系起来加以考察.47).()(2)(4121kkkxfxfxfh用复化梯形公式求得该子区间上的积分值为 每个子区间 经过二分只增加了一个分点,1kkxx),(21121
22、kkkxxx这里 代表二分前的步长.nabh 将每个子区间上的积分值相加得,)(2)()(410211012nkknkkknxfhxfxfhT48 例例2 2.sin10dxxxI 解解先对整个区间 使用梯形公式.1,0从而利用式(3.2)可导出下列递推公式.)(22110212nkknnxfhTT（4.1）计算积分值 对于函数,sin)(xxxf定义它在 的值0 x,1)0(f,8414709.0)1(f而由梯形公式49.9207355.0)1()0(211ffT将区间二等分，求出中点的函数值,9588510.0)21(f利用递推公式(4.1)，有.9397933.0)21(212112fT
23、T进一步二分求积区间，并计算新分点上的函数值.9088516.0)43(,9896158.0)41(ff再利用式(4.1)，有 50.9445135.0)43()41(412124ffTT这样不断二分下去，计算结果见下表.9460831.09460830.09460827.09460815.09460769.01098769460596.09459850.09456909.09445135.09397933.054321nnTkTk3-表4 它表明用复化梯形公式计算积分 要达到7位有效数字的精度需要二分区间10次，即要有分点1025个，计算量很大.I51 4.4.2 4.4.2 龙贝格算法龙贝
24、格算法 梯形法计算简单但收敛慢，本节讨论如何提高收敛速度以节省计算量.根据复化梯形公式的余项表达式(3.3);,(),(122bafhabTIn).,(),(21222bafhabTIn 假定),()(ff .412nnTITI则有52移项整理，得).(3122nnnTTTI（4.2）由此可见，只有二分前后的两个积分值 与 相当接近，就可以保证计算结果 的误差很小.nTnT2nT2 这样直接用计算结果来估计误差的方法通常称作误差的事后估计法事后估计法.按式(4.2)，积分近似值 的误差大致等于 nT2),(312nnTT若用这个误差值作为 的一种补偿，可以期望所得到的nT253)(3122nn
25、nTTTT可能是更好的结果.直接验证知，按公式(4.3)组合得到的近似值恰为 ,nS即.31342nnnTTS（4.4）也就是说，利用梯形法二分前后的两个积分值 与 ,nTnT2按式(4.3)做线性组合，结果得到辛普森法的积分值.nS（4.3）nnTT3134254 再考察辛普森法，按误差公式(3.6)，其截断误差大致与 成正比，4h,1612nnSISI由此得到.15115162nnSSI因此，若将步长折半则误差将减至原有误差的1/16，即有不难直接验证，上式右端的值等于 ,为复化柯特斯公式，nC).(6hO它的精度为 55 就是说，用辛普森法二分前后两个积分值 与 ，nSnS2可得到 的近
26、似误差估计为nS2),(1512nnSS并且按上式做线性组合可得复化柯特斯的积分值 ，nC.15115162nnnSSC（4.5）重复同样的手续，依据柯特斯法的误差阶为 ，6h.63163642nnnCCR（4.6）即有可进一步导出下列龙贝格龙贝格(Romberg)公式公式：56 在变步长的过程中运用公式(4.4)，(4.5)和(4.6)，就能将粗糙的梯形值 逐步加工成精度较高的辛普森值nTnS柯特斯值 和龙贝格值 .nCnR 例例3 3计算结果见下表(代表二分次数)：k9460831.09460831.09460833.09456909.039460830.09460869.09445135
27、.029461459.09397933.019207355.003212222kkkkRCSTk4-表4用加速公式(4.4)，(4.5)和(4.6)加工例2得到的梯形值，57 可以看到，这里利用二分3次的数据(它们的精度都很差，只有两三位有效数字)，通过三次加速求得 9460831.01R这个结果的每一位数字都是有效数字，可见加速的效果是十分显著的.58 4.4.3 4.4.3 理查森外推加速法理查森外推加速法 由梯形公式出发，将区间 逐次二分可提高求积公式的精度，上述加速过程还可继续下去，设,ba),(122fhabTIn 若记),(hTTn当区间 分为 等分时，,ban2),(12)(2f
28、habIhT 梯形公式余项可展成级数形式，即,ba.nabh并且有),2(2hTTn则有 ,)0()(lim0IThTh59 定理定理4 4 设 则有,)(baCxf,)(24221llhhhIhT（4.7）其中系数 与 无关.),2,1(llh 定理4表明 是 阶，IhT)()(2hO代替 ，有 h,2164)2(24221llhhhIhT（4.8）若用4乘(4.8)式，减去(4.7)式再除3记之为),(1hT2h在(4.7)中，若用则得 603)(24)(1hThThT这里 以及后面将出现的 均为与 无关的系数，,21kk,h这样构造的 与积分值 近似的阶为 .)(1hTI)(4hO（4.
29、9）,6241hhI 比较(4.9)与(4.4)可知，这样构造的序列),2(),(11hThT就是辛普森公式序列.,2nnSS 又根据(4.9)，有,16)2(8362411hhhIhT61若令),(161)2(1516)(112hThThT则又可进一步从余项展开式中消去 项，而有 4h.)(82612hhIhT这样构造出的 ，其实就是柯特斯公式序列.)(2hT它与积分值 的逼近阶为 I).(6hO 如此继续下去，每加速一次，误差的量级便提高2阶.62 一般地，若记 则有),()(0hThT).(1412144)(11hThThTmmmmmm（4.10）经过 次加速后，),2,1(mm.)()
30、2(22)1(21mmmhhIhT（4.11）上述处理方法通常称为理查森外推加速方法理查森外推加速方法.设以 表示二分 次后求得的梯形值，且以 表示)(0kTk)(kmT序列 的 次加速值，则依递推公式(4.10)可得)(0kTm余项便取下列形式：63计算过程：(1)取,0abhk 令k1(记区间 的二分次数).k,ba (2)求梯形值,20kabT即按递推公式(4.1)计算.)(0kT (3)求加速值，按公式(4.12)逐个求出T表(见表4-5)的第 行其余各元素 k).,2,1()(kjTjkj公式(4.12)也称为龙贝格求积算法龙贝格求积算法.（4.12）).,2,1(141144)(1
31、)1(1)(kTTTkmmkmmmkm).()(2)0(0bfafhT求64)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表T5-4表 (4)若 (预先给定的精度)，则终止计算，)0(1)0(kkTT否则令 转(2)继续计算.kk1;)0(ITk并取65 可以证明，如果 充分光滑，那么T表每一列的元素及对角线元素均收敛到所求的积分值 ，即)(xfI，）（固定mITkmk)(li
32、m 对于 不充分光滑的函数也可用龙贝格算法计算，)(xf.lim)0(ITmm只是收敛慢一些，这时也可以直接使用复化辛普森公式计算.66 例例4 4 解解在 上仅是一次连续可微，2/3)(xxf1,0.102/3dxxI用龙贝格算法计算积分 用龙贝格算法计算结果见表4-6.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)(0kkkkkkTTTTTTk6-表467 从表中看到用龙贝格算到 的精度与辛普森求积精度相当.这里 的精确值为 5kI.4.0
